<p><b>ringler@lanl.gov</b> 2010-03-25 10:59:54 -0600 (Thu, 25 Mar 2010)</p><p><br>
added grid generator based on hexagons that included necessary<br>
changes for ocean simulations.<br>
</p><hr noshade><pre><font color="gray">Copied: branches/ocean_projects/grid_gen_ocean/hex_periodic (from rev 144, trunk/grid_gen/periodic_hex)
Modified: branches/ocean_projects/grid_gen_ocean/hex_periodic/Makefile
===================================================================
--- trunk/grid_gen/periodic_hex/Makefile 2010-03-17 20:33:00 UTC (rev 144)
+++ branches/ocean_projects/grid_gen_ocean/hex_periodic/Makefile 2010-03-25 16:59:54 UTC (rev 160)
@@ -6,19 +6,19 @@
#LDFLAGS = -g -C
# pgf90
-FC = pgf90
-CC = pgcc
-FFLAGS = -r8 -O3
-CFLAGS = -O3
-LDFLAGS = -O3
-
-# ifort
-#FC = ifort
-#CC = icc
-#FFLAGS = -real-size 64 -O3
+#FC = pgf90
+#CC = pgcc
+#FFLAGS = -r8 -O3
#CFLAGS = -O3
#LDFLAGS = -O3
+# ifort
+FC = ifort
+CC = icc
+FFLAGS = -real-size 64 -O1 -g -traceback -check all
+CFLAGS = -O3
+LDFLAGS = -O1 -g -traceback -check all
+
# absoft
#FC = f90
#CC = gcc
Modified: branches/ocean_projects/grid_gen_ocean/hex_periodic/module_write_netcdf.F
===================================================================
--- trunk/grid_gen/periodic_hex/module_write_netcdf.F 2010-03-17 20:33:00 UTC (rev 144)
+++ branches/ocean_projects/grid_gen_ocean/hex_periodic/module_write_netcdf.F 2010-03-25 16:59:54 UTC (rev 160)
@@ -51,6 +51,7 @@
integer :: wrVarIDh_s
integer :: wrVarIDu
integer :: wrVarIDuBC
+ integer :: wrVarIDu_src
integer :: wrVarIDv
integer :: wrVarIDh
integer :: wrVarIDvh
@@ -215,6 +216,9 @@
nferr = nf_def_var(wr_ncid, 'uBC', NF_INT, 2, dimlist, wrVarIDuBC)
dimlist( 1) = wrDimIDnVertLevels
dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'u_src', NF_DOUBLE, 2, dimlist, wrVarIDu_src)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
dimlist( 3) = wrDimIDTime
nferr = nf_def_var(wr_ncid, 'v', NF_DOUBLE, 3, dimlist, wrVarIDv)
dimlist( 1) = wrDimIDnVertLevels
@@ -289,6 +293,7 @@
fVertex, &
h_s, &
uBC, &
+ u_src, &
u, &
v, &
h, &
@@ -343,6 +348,7 @@
real (kind=8), dimension(:), intent(in) :: fVertex
real (kind=8), dimension(:), intent(in) :: h_s
integer, dimension(:,:), intent(in) :: uBC
+ real (kind=8), dimension(:,:), intent(in) :: u_src
real (kind=8), dimension(:,:,:), intent(in) :: u
real (kind=8), dimension(:,:,:), intent(in) :: v
real (kind=8), dimension(:,:,:), intent(in) :: h
@@ -537,7 +543,12 @@
start2(2) = 1
count2( 1) = wrLocalnVertLevels
count2( 2) = wrLocalnEdges
- nferr = nf_put_vara_int(wr_ncid, wrVarIDuBC, start2, count2, u)
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDuBC, start2, count2, uBC)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDu_src, start2, count2, u_src)
start3(3) = time
count3( 1) = wrLocalnVertLevels
Modified: branches/ocean_projects/grid_gen_ocean/hex_periodic/namelist.input
===================================================================
--- trunk/grid_gen/periodic_hex/namelist.input 2010-03-17 20:33:00 UTC (rev 144)
+++ branches/ocean_projects/grid_gen_ocean/hex_periodic/namelist.input 2010-03-25 16:59:54 UTC (rev 160)
@@ -1,8 +1,8 @@
&periodic_grid
- nx = 500,
- ny = 500,
- dc = 10000.,
- nVertLevels = 1,
+ nx = 125,
+ ny = 250,
+ dc = 20000.,
+ nVertLevels = 4,
nTracers = 1,
nproc = 2, 4, 8,
/
Modified: branches/ocean_projects/grid_gen_ocean/hex_periodic/periodic_grid.F
===================================================================
--- trunk/grid_gen/periodic_hex/periodic_grid.F 2010-03-17 20:33:00 UTC (rev 144)
+++ branches/ocean_projects/grid_gen_ocean/hex_periodic/periodic_grid.F 2010-03-25 16:59:54 UTC (rev 160)
@@ -8,6 +8,10 @@
! specify grid resolution with variable dc (meters)
! real (kind=8), parameter :: dc = 10000.0 ! Distance between cell centers
+ real (kind=8) :: ymid, f0, beta, ytmp, ymax, stress
+ real (kind=8), allocatable, dimension(:,:,:) :: utmp
+ integer :: iVert1, iEdge1, iEdge2, iEdge3
+
real (kind=8), parameter :: pi = 3.141592653589793
real (kind=8), parameter :: ONE = 1.0_8
real (kind=8), parameter :: TWO = 2.0_8
@@ -25,7 +29,7 @@
real (kind=8), allocatable, dimension(:) :: latCell, lonCell, xCell, yCell, zCell
real (kind=8), allocatable, dimension(:) :: latEdge, lonEdge, xEdge, yEdge, zEdge
real (kind=8), allocatable, dimension(:) :: latVertex, lonVertex, xVertex, yVertex, zVertex
- real (kind=8), allocatable, dimension(:,:) :: weightsOnEdge, kiteAreasOnVertex
+ real (kind=8), allocatable, dimension(:,:) :: weightsOnEdge, kiteAreasOnVertex, u_src
real (kind=8), allocatable, dimension(:) :: fEdge, fVertex, h_s
real (kind=8), allocatable, dimension(:,:,:) :: u, v, h, vh, circulation, vorticity, ke
real (kind=8), allocatable, dimension(:,:,:,:) :: tracers
@@ -89,6 +93,7 @@
allocate(fVertex(nVertices))
allocate(h_s(nCells))
allocate(uBC(nVertLevels, nEdges))
+ allocate(u_src(nVertLevels, nEdges))
allocate(u(nVertLevels,nEdges,1))
allocate(v(nVertLevels,nEdges,1))
@@ -99,6 +104,8 @@
allocate(ke(nVertLevels,nCells,1))
allocate(tracers(nTracers,nVertLevels,nCells,1))
+ allocate(utmp(0:nx+1, 0:ny+1, 3))
+
do iRow = 1, ny
do iCol = 1, nx
iCell = cellIdx(iCol,iRow)
@@ -293,17 +300,74 @@
tracers(:,:,:,:) = 0.0
h(:,:,:) = 1.0
- do i=1,nCells
- r = sqrt((xCell(i) - (nx/2)*(10.0*dc))**2.0 + (yCell(i) - (ny/2)*(10.0*dc))**2.0)
- if (r < 10.0*10.0*dc) then
- tracers(1,1,i,1) = (20.0 / 2.0) * (1.0 + cos(pi*r/(10.0*10.0*dc))) + 0.0
- h(1,i,1) = 1.0 + 0.1*cos(pi*r/(20.0*10.0*dc))
- else
- tracers(1,1,i,1) = 0.0
- h(1,i,1) = 1.0
- end if
- end do
+ h = 100.0/nVertLevels
+ stress = 0.1
+ beta = 5.0e-11
+ f0 = 1.4e-4
+ ymid = dc*real(ny/2)*sqrt(THREE) / TWO
+
+ fVertex = 0.0
+ do i = 1,nVertices
+ fVertex(i) = f0 + (yVertex(i) - ymid) * beta
+ enddo
+
+ fEdge = 0.0
+ do i = 1,nEdges
+ fEdge(i) = f0 + (yEdge(i) - ymid) * beta
+ enddo
+
+ do iRow=1,ny
+ do iCol=1,nx
+ ymax = dc*real(ny)*sqrt(THREE) / TWO
+ ytmp = dc*real(iRow)*sqrt(THREE) / TWO - 0.0
+ utmp(iCol, iRow, 1) = - stress * cos ( TWO * pi * ytmp / ymax)
+
+ ytmp = dc*real(iRow)*sqrt(THREE) / TWO - 0.5 * dc * sin(pi/THREE)
+ utmp(iCol, iRow, 2) = - 0.5*stress * cos ( TWO * pi * ytmp / ymax)
+
+ ytmp = dc*real(iRow)*sqrt(THREE) / TWO - 0.5 * dc * sin(pi/THREE)
+ utmp(iCol, iRow, 3) = + 0.5*stress * cos ( TWO * pi * ytmp / ymax)
+
+ if(iCol.eq.1) write(6,*) utmp(iCol,iRow,1)
+
+ enddo
+ enddo
+
+ u_src = 0.0
+ do iRow=1,ny
+ do iCol=1,nx
+ iCell=cellIdx(iCol,iRow)
+ iEdge1=3*(iCell-1)+1
+ iEdge2=3*(iCell-1)+2
+ iEdge3=3*(iCell-1)+3
+ u_src(1,iEdge1) = utmp(iCol,iRow,1)
+ u_src(1,iEdge2) = utmp(iCol,iRow,2)
+ u_src(1,iEdge3) = utmp(iCol,iRow,3)
+ enddo
+ enddo
+
+ call enforce_uBC(u(:,:,:), uBC(:,:), xCell(:), yCell(:), zCell(:), nCells, nEdges, nVertLevels, dc, cellsOnEdge(:,:))
+
+ j = 0
+ do i=1,nEdges
+ if(uBC(1,i).gt.0) j=j+1
+ enddo
+
+ write(6,*) ' summary of data '
+ write(6,10) ' grid res ', dc
+ write(6,11) ' domain ', nx, ny
+ write(6,10) ' thickness ', minval(h), maxval(h)
+ write(6,10) ' fVertex ', minval(fVertex), maxval(fVertex)
+ write(6,10) ' fEdge ', minval(fEdge), maxval(fEdge)
+ write(6,10) ' u_src ', minval(u_src), maxval(u_src)
+ write(6,10) ' u ', minval(u), maxval(u)
+ write(6,10) ' xCell ', minval(xCell), maxval(xCell)
+ write(6,10) ' yCell ', minval(yCell), maxval(yCell)
+ write(6,10 ) ' uBC ',float(j), float(j)/float(nx*ny)
+ 10 format(a12,2e15.5)
+ 11 format(a12,2i5)
+
!
! Write grid to grid.nc file
!
@@ -334,6 +398,7 @@
fVertex, &
h_s, &
uBC, &
+ u_src, &
u, &
v, &
h, &
@@ -381,21 +446,6 @@
end program hexagonal_periodic_grid
-subroutine enforce_uBC(u, uBC, xCell, yCell, zCell, nCells, nEdges, nVertLevels, dc)
-! this suboutine provides a hook into uBC. the uBC field is read into the ocean
-! model and used to enforce boundary conditions on the velocity field.
-! uBC is written to the grid.nc file, even if the forward model does not use it.
-
-real (kind=8), intent(in) :: dc
-real (kind=8), intent(inout), dimension(nVertLevels, nEdges, 1) :: u
-real (kind=8), intent(in), dimension(nCells) :: xCell, yCell, zCell
-integer, intent(inout), dimension(nVertLevels, nEdges) :: uBC
-
-uBC = -10
-
-end subroutine enforce_uBC
-
-
subroutine decompose_nproc(nproc, nprocx, nprocy)
implicit none
@@ -409,3 +459,33 @@
end do
end subroutine decompose_nproc
+
+
+subroutine enforce_uBC(u, uBC, xCell, yCell, zCell, nCells, nEdges, nVertLevels, dc, cellsOnEdge)
+
+implicit none
+
+real (kind=8), intent(in) :: dc
+real (kind=8), intent(inout), dimension(nVertLevels, nEdges, 1) :: u
+real (kind=8), intent(in), dimension(nCells) :: xCell, yCell, zCell
+integer, intent(inout), dimension(nVertLevels, nEdges) :: uBC
+integer, intent(in) :: cellsOnEdge(2,nEdges), nCells, nEdges, nVertLevels
+
+real (kind=8) :: d, b(3), t(3)
+integer :: i1, i2, iEdge
+
+uBC = -1
+do iEdge=1,nEdges
+ i1 = cellsOnEdge(1,iEdge)
+ i2 = cellsOnEdge(2,iEdge)
+ b(1) = xCell(i1); b(2) = yCell(i1); b(3) = zCell(i1)
+ t(1) = xCell(i2); t(2) = yCell(i2); t(3) = zCell(i2)
+ d = sqrt( (b(1)-t(1))**2 + (b(2)-t(2))**2 + (b(3)-t(3))**2 )
+ if(d.gt.1.1*dc) then
+ uBC(:,iEdge) = 1
+ u(:,iEdge,:) = 0.0
+ endif
+ 10 format(2f10.2)
+enddo
+
+end subroutine enforce_uBC
</font>
</pre>