<p><b>mpetersen@lanl.gov</b> 2010-03-23 13:17:10 -0600 (Tue, 23 Mar 2010)</p><p>Move grid generation code to branch project space.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/grid_gen_ocean/quad_periodic/Makefile
===================================================================
--- branches/ocean_projects/grid_gen_ocean/quad_periodic/Makefile (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/quad_periodic/Makefile 2010-03-23 19:17:10 UTC (rev 155)
@@ -0,0 +1,62 @@
+# IBM with Xlf compilers
+#FC = xlf90
+#CC = xlc
+#FFLAGS = -qrealsize=8 -g -C
+#CFLAGS = -g
+#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
+CFLAGS = -O3
+LDFLAGS = -O3
+NETCDF=/usr/projects/climate/bzhao/netcdf-3.6.1
+
+# absoft
+#FC = f90
+#CC = gcc
+#FFLAGS = -dp -O3
+#CFLAGS = -O3
+#LDFLAGS = -O3
+
+
+CPP = cpp -C -P -traditional
+CPPFLAGS =
+CPPINCLUDES =
+INCLUDES = -I$(NETCDF)/include
+LIBS = -L$(NETCDF)/lib -lnetcdf
+
+RM = rm -f
+
+##########################
+
+.SUFFIXES: .F .o
+
+
+OBJS = periodic_grid.o \
+ module_cell_indexing.o \
+ module_write_netcdf.o
+
+all: periodic_grid
+
+periodic_grid.o: module_cell_indexing.o module_write_netcdf.o
+
+periodic_grid: $(OBJS)
+ $(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
+
+clean:
+ $(RM) *.o *.mod periodic_grid
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(INCLUDES)
+ $(RM) $*.f90
Added: branches/ocean_projects/grid_gen_ocean/quad_periodic/module_cell_indexing.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/quad_periodic/module_cell_indexing.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/quad_periodic/module_cell_indexing.F 2010-03-23 19:17:10 UTC (rev 155)
@@ -0,0 +1,142 @@
+module cell_indexing
+
+! this subroutine provide index mapping for quad meshes dimensioned (nx, ny)
+
+ integer, parameter :: maxEdges = 4
+
+ integer :: nx, ny, nVertLevels, nTracers, nVertexDegree
+ integer, dimension(20) :: nproc
+
+
+ contains
+
+
+ subroutine cell_indexing_read_nl()
+
+ implicit none
+
+ namelist /periodic_grid/ nx, ny, nVertLevels, nTracers, nproc, nVertexDegree
+
+ nx = 200
+ ny = 200
+ nVertLevels = 4
+ nTracers = 2
+ nproc(:) = -1
+
+ open(20,file='namelist.input',status='old')
+ read(20,periodic_grid)
+ close(20)
+
+ end subroutine cell_indexing_read_nl
+
+
+ subroutine cellColRow(idx, iCol, iRow)
+
+ implicit none
+
+ integer, intent(in) :: idx
+ integer, intent(out) :: iCol, iRow
+
+ iRow = ((idx-1) / nx) + 1
+ iCol = mod((idx-1), nx) + 1
+
+ end subroutine cellColRow
+
+
+ integer function cellIdx(iCol, iRow)
+
+ implicit none
+
+ integer, intent(in) :: iCol, iRow
+
+ cellIdx = (iRow-1)*nx + iCol
+
+ end function cellIdx
+
+
+ integer function cellOnCell(iCol, iRow, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCol, iRow, neighborNumber
+
+ integer :: mx, px, my, py
+
+ mx = iCol - 1
+ if (mx == 0) mx = nx
+ my = iRow - 1
+ if (my == 0) my = ny
+ px = iCol + 1
+ if (px == nx + 1) px = 1
+ py = iRow + 1
+ if (py == ny + 1) py = 1
+
+ if (neighborNumber == 1) then
+ cellOnCell = cellIdx(mx, iRow)
+ else if (neighborNumber == 2) then
+ cellOnCell = cellIdx(iCol, my)
+ else if (neighborNumber == 3) then
+ cellOnCell = cellIdx(px, iRow)
+ else if (neighborNumber == 4) then
+ cellOnCell = cellIdx(iCol, py)
+ endif
+
+ end function cellOnCell
+
+
+ integer function edgeOnCell(iCell, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCell, neighborNumber
+
+ integer :: myRow, myCol
+
+ call cellColRow(iCell, myCol, myRow)
+
+ if (neighborNumber == 1) then
+ edgeOnCell = 2*(iCell - 1) + 1
+ else if (neighborNumber == 2) then
+ edgeOnCell = 2*(iCell - 1) + 2
+ else if (neighborNumber == 3) then
+ edgeOnCell = 2*(cellOnCell(myCol, myRow, 3) - 1) + 1
+ else if (neighborNumber == 4) then
+ edgeOnCell = 2*(cellOnCell(myCol, myRow, 4) - 1) + 2
+ end if
+
+ end function edgeOnCell
+
+
+ integer function vertexOnCell(iCell, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCell, neighborNumber
+
+ integer :: myRow, myCol, mx, my, px, py
+
+ call cellColRow(iCell, myCol, myRow)
+
+ mx = myCol - 1
+ if (mx == 0) mx = nx
+ my = myRow - 1
+ if (my == 0) my = ny
+ px = myCol + 1
+ if (px == nx + 1) px = 1
+ py = myRow + 1
+ if (py == ny + 1) py = 1
+
+ if (neighborNumber == 1) then
+ vertexOnCell = cellIdx(myCol,myRow)
+ else if (neighborNumber == 2) then
+ vertexOnCell = cellIdx(myCol,my)
+ else if (neighborNumber == 3) then
+ vertexOnCell = cellIdx(px,my)
+ else if (neighborNumber == 4) then
+ vertexOnCell = cellIdx(px,myRow)
+ end if
+
+ end function vertexOnCell
+
+
+end module cell_indexing
Added: branches/ocean_projects/grid_gen_ocean/quad_periodic/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/quad_periodic/module_write_netcdf.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/quad_periodic/module_write_netcdf.F 2010-03-23 19:17:10 UTC (rev 155)
@@ -0,0 +1,621 @@
+module write_netcdf
+
+ integer :: wr_ncid
+ integer :: wrDimIDTime
+ integer :: wrDimIDnCells
+ integer :: wrDimIDnEdges
+ integer :: wrDimIDnVertices
+ integer :: wrDimIDmaxEdges
+ integer :: wrDimIDmaxEdges2
+ integer :: wrDimIDTWO
+ integer :: wrDimIDvertexDegree
+ integer :: wrDimIDnVertLevels
+ integer :: wrDimIDnTracers
+ integer :: wrVarIDlatCell
+ integer :: wrVarIDlonCell
+ integer :: wrVarIDxCell
+ integer :: wrVarIDyCell
+ integer :: wrVarIDzCell
+ integer :: wrVarIDindexToCellID
+ integer :: wrVarIDlatEdge
+ integer :: wrVarIDlonEdge
+ integer :: wrVarIDxEdge
+ integer :: wrVarIDyEdge
+ integer :: wrVarIDzEdge
+ integer :: wrVarIDindexToEdgeID
+ integer :: wrVarIDlatVertex
+ integer :: wrVarIDlonVertex
+ integer :: wrVarIDxVertex
+ integer :: wrVarIDyVertex
+ integer :: wrVarIDzVertex
+ integer :: wrVarIDindexToVertexID
+ integer :: wrVarIDcellsOnEdge
+ integer :: wrVarIDnEdgesOnCell
+ integer :: wrVarIDnEdgesOnEdge
+ integer :: wrVarIDedgesOnCell
+ integer :: wrVarIDedgesOnEdge
+ integer :: wrVarIDweightsOnEdge
+ integer :: wrVarIDdvEdge
+ integer :: wrVarIDdcEdge
+ integer :: wrVarIDangleEdge
+ integer :: wrVarIDareaCell
+ integer :: wrVarIDareaTriangle
+ integer :: wrVarIDcellsOnCell
+ integer :: wrVarIDverticesOnCell
+ integer :: wrVarIDverticesOnEdge
+ integer :: wrVarIDedgesOnVertex
+ integer :: wrVarIDcellsOnVertex
+ integer :: wrVarIDkiteAreasOnVertex
+ integer :: wrVarIDfEdge
+ integer :: wrVarIDfVertex
+ integer :: wrVarIDh_s
+ integer :: wrVarIDu
+ integer :: wrVarIDuBC
+ integer :: wrVarIDu_src
+ integer :: wrVarIDv
+ integer :: wrVarIDh
+ integer :: wrVarIDvh
+ integer :: wrVarIDcirculation
+ integer :: wrVarIDvorticity
+ integer :: wrVarIDke
+ integer :: wrVarIDtracers
+
+ integer :: wrLocalnCells
+ integer :: wrLocalnEdges
+ integer :: wrLocalnVertices
+ integer :: wrLocalmaxEdges
+ integer :: wrLocalnVertLevels
+ integer :: wrLocalnTracers
+ integer :: wrLocalnVertexDegree
+
+ contains
+
+ subroutine write_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ nVertLevels, &
+ nTracers, &
+ nVertexDegree &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertices
+ integer, intent(in) :: maxEdges
+ integer, intent(in) :: nVertLevels
+ integer, intent(in) :: nTracers
+ integer, intent(in) :: nVertexDegree
+
+ integer :: nferr
+ integer, dimension(10) :: dimlist
+
+
+ wrLocalnCells = nCells
+ wrLocalnEdges = nEdges
+ wrLocalnVertices = nVertices
+ wrLocalmaxEdges = maxEdges
+ wrLocalnVertLevels = nVertLevels
+ wrLocalnTracers = nTracers
+ wrLocalnVertexDegree = nVertexDegree
+
+ nferr = nf_create('grid.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
+
+ !
+ ! Define dimensions
+ !
+ nferr = nf_def_dim(wr_ncid, 'nCells', nCells, wrDimIDnCells)
+ nferr = nf_def_dim(wr_ncid, 'nEdges', nEdges, wrDimIDnEdges)
+ nferr = nf_def_dim(wr_ncid, 'nVertices', nVertices, wrDimIDnVertices)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges', maxEdges, wrDimIDmaxEdges)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges2', 2*maxEdges, wrDimIDmaxEdges2)
+ nferr = nf_def_dim(wr_ncid, 'TWO', 2, wrDimIDTWO)
+! nferr = nf_def_dim(wr_ncid, 'THREE', 3, wrDimIDTHREE)
+ nferr = nf_def_dim(wr_ncid, 'vertexDegree', nVertexDegree, wrDimIDvertexDegree)
+ nferr = nf_def_dim(wr_ncid, 'nVertLevels', nVertLevels, wrDimIDnVertLevels)
+ nferr = nf_def_dim(wr_ncid, 'nTracers', nTracers, wrDimIDnTracers)
+ nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
+
+ !
+ ! Define variables
+ !
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'latCell', NF_DOUBLE, 1, dimlist, wrVarIDlatCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'lonCell', NF_DOUBLE, 1, dimlist, wrVarIDlonCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'xCell', NF_DOUBLE, 1, dimlist, wrVarIDxCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'yCell', NF_DOUBLE, 1, dimlist, wrVarIDyCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'zCell', NF_DOUBLE, 1, dimlist, wrVarIDzCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'indexToCellID', NF_INT, 1, dimlist, wrVarIDindexToCellID)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'latEdge', NF_DOUBLE, 1, dimlist, wrVarIDlatEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'lonEdge', NF_DOUBLE, 1, dimlist, wrVarIDlonEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'xEdge', NF_DOUBLE, 1, dimlist, wrVarIDxEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'yEdge', NF_DOUBLE, 1, dimlist, wrVarIDyEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'zEdge', NF_DOUBLE, 1, dimlist, wrVarIDzEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'indexToEdgeID', NF_INT, 1, dimlist, wrVarIDindexToEdgeID)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'latVertex', NF_DOUBLE, 1, dimlist, wrVarIDlatVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'lonVertex', NF_DOUBLE, 1, dimlist, wrVarIDlonVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'xVertex', NF_DOUBLE, 1, dimlist, wrVarIDxVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'yVertex', NF_DOUBLE, 1, dimlist, wrVarIDyVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'zVertex', NF_DOUBLE, 1, dimlist, wrVarIDzVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'indexToVertexID', NF_INT, 1, dimlist, wrVarIDindexToVertexID)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'cellsOnEdge', NF_INT, 2, dimlist, wrVarIDcellsOnEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnCell', NF_INT, 1, dimlist, wrVarIDnEdgesOnCell)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnEdge', NF_INT, 1, dimlist, wrVarIDnEdgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'edgesOnCell', NF_INT, 2, dimlist, wrVarIDedgesOnCell)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'edgesOnEdge', NF_INT, 2, dimlist, wrVarIDedgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'weightsOnEdge', NF_DOUBLE, 2, dimlist, wrVarIDweightsOnEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dvEdge', NF_DOUBLE, 1, dimlist, wrVarIDdvEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dcEdge', NF_DOUBLE, 1, dimlist, wrVarIDdcEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'angleEdge', NF_DOUBLE, 1, dimlist, wrVarIDangleEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'areaCell', NF_DOUBLE, 1, dimlist, wrVarIDareaCell)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'areaTriangle', NF_DOUBLE, 1, dimlist, wrVarIDareaTriangle)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'cellsOnCell', NF_INT, 2, dimlist, wrVarIDcellsOnCell)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'verticesOnCell', NF_INT, 2, dimlist, wrVarIDverticesOnCell)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'verticesOnEdge', NF_INT, 2, dimlist, wrVarIDverticesOnEdge)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'edgesOnVertex', NF_INT, 2, dimlist, wrVarIDedgesOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'cellsOnVertex', NF_INT, 2, dimlist, wrVarIDcellsOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'kiteAreasOnVertex', NF_DOUBLE, 2, dimlist, wrVarIDkiteAreasOnVertex)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'fEdge', NF_DOUBLE, 1, dimlist, wrVarIDfEdge)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'fVertex', NF_DOUBLE, 1, dimlist, wrVarIDfVertex)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'h_s', NF_DOUBLE, 1, dimlist, wrVarIDh_s)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'u', NF_DOUBLE, 3, dimlist, wrVarIDu)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ 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
+ dimlist( 2) = wrDimIDnCells
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'h', NF_DOUBLE, 3, dimlist, wrVarIDh)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'vh', NF_DOUBLE, 3, dimlist, wrVarIDvh)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnVertices
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'circulation', NF_DOUBLE, 3, dimlist, wrVarIDcirculation)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnVertices
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'vorticity', NF_DOUBLE, 3, dimlist, wrVarIDvorticity)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnCells
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'ke', NF_DOUBLE, 3, dimlist, wrVarIDke)
+ dimlist( 1) = wrDimIDnTracers
+ dimlist( 2) = wrDimIDnVertLevels
+ dimlist( 3) = wrDimIDnCells
+ dimlist( 4) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'tracers', NF_DOUBLE, 4, dimlist, wrVarIDtracers)
+
+ nferr = nf_enddef(wr_ncid)
+
+ end subroutine write_netcdf_init
+
+
+ subroutine write_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ xCell, &
+ yCell, &
+ zCell, &
+ indexToCellID, &
+ latEdge, &
+ lonEdge, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ indexToEdgeID, &
+ latVertex, &
+ lonVertex, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex, &
+ h_s, &
+ uBC, &
+ u_src, &
+ u, &
+ v, &
+ h, &
+ vh, &
+ circulation, &
+ vorticity, &
+ ke, &
+ tracers &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: time
+ real (kind=8), dimension(:), intent(in) :: latCell
+ real (kind=8), dimension(:), intent(in) :: lonCell
+ real (kind=8), dimension(:), intent(in) :: xCell
+ real (kind=8), dimension(:), intent(in) :: yCell
+ real (kind=8), dimension(:), intent(in) :: zCell
+ integer, dimension(:), intent(in) :: indexToCellID
+ real (kind=8), dimension(:), intent(in) :: latEdge
+ real (kind=8), dimension(:), intent(in) :: lonEdge
+ real (kind=8), dimension(:), intent(in) :: xEdge
+ real (kind=8), dimension(:), intent(in) :: yEdge
+ real (kind=8), dimension(:), intent(in) :: zEdge
+ integer, dimension(:), intent(in) :: indexToEdgeID
+ real (kind=8), dimension(:), intent(in) :: latVertex
+ real (kind=8), dimension(:), intent(in) :: lonVertex
+ real (kind=8), dimension(:), intent(in) :: xVertex
+ real (kind=8), dimension(:), intent(in) :: yVertex
+ real (kind=8), dimension(:), intent(in) :: zVertex
+ integer, dimension(:), intent(in) :: indexToVertexID
+ integer, dimension(:,:), intent(in) :: cellsOnEdge
+ integer, dimension(:), intent(in) :: nEdgesOnCell
+ integer, dimension(:), intent(in) :: nEdgesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnCell
+ integer, dimension(:,:), intent(in) :: edgesOnEdge
+ real (kind=8), dimension(:,:), intent(in) :: weightsOnEdge
+ real (kind=8), dimension(:), intent(in) :: dvEdge
+ real (kind=8), dimension(:), intent(in) :: dcEdge
+ real (kind=8), dimension(:), intent(in) :: angleEdge
+ real (kind=8), dimension(:), intent(in) :: areaCell
+ real (kind=8), dimension(:), intent(in) :: areaTriangle
+ integer, dimension(:,:), intent(in) :: cellsOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnVertex
+ integer, dimension(:,:), intent(in) :: cellsOnVertex
+ real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
+ real (kind=8), dimension(:), intent(in) :: fEdge
+ 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
+ real (kind=8), dimension(:,:,:), intent(in) :: vh
+ real (kind=8), dimension(:,:,:), intent(in) :: circulation
+ real (kind=8), dimension(:,:,:), intent(in) :: vorticity
+ real (kind=8), dimension(:,:,:), intent(in) :: ke
+ real (kind=8), dimension(:,:,:,:), intent(in) :: tracers
+
+ integer :: nferr
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ start1(1) = 1
+
+ start2(1) = 1
+ start2(2) = 1
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+
+ start4(1) = 1
+ start4(2) = 1
+ start4(3) = 1
+ start4(4) = 1
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatCell, start1, count1, latCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonCell, start1, count1, lonCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxCell, start1, count1, xCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyCell, start1, count1, yCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzCell, start1, count1, zCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToCellID, start1, count1, indexToCellID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatEdge, start1, count1, latEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonEdge, start1, count1, lonEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxEdge, start1, count1, xEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyEdge, start1, count1, yEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzEdge, start1, count1, zEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToEdgeID, start1, count1, indexToEdgeID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatVertex, start1, count1, latVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonVertex, start1, count1, lonVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxVertex, start1, count1, xVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyVertex, start1, count1, yVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzVertex, start1, count1, zVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToVertexID, start1, count1, indexToVertexID)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnEdge, start2, count2, cellsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnCell, start2, count2, edgesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnEdge, start2, count2, edgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDweightsOnEdge, start2, count2, weightsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdvEdge, start1, count1, dvEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdcEdge, start1, count1, dcEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDangleEdge, start1, count1, angleEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaCell, start1, count1, areaCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaTriangle, start1, count1, areaTriangle)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnCell, start2, count2, cellsOnCell)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnCell, start2, count2, verticesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnEdge, start2, count2, verticesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfEdge, start1, count1, fEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfVertex, start1, count1, fVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDh_s, start1, count1, h_s)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnEdges
+ 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
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDu, start3, count3, u)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDv, start3, count3, v)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnCells
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDh, start3, count3, h)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDvh, start3, count3, vh)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnVertices
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDcirculation, start3, count3, circulation)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnVertices
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDvorticity, start3, count3, vorticity)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnCells
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDke, start3, count3, ke)
+
+ start4(4) = time
+ count4( 1) = wrLocalnTracers
+ count4( 2) = wrLocalnVertLevels
+ count4( 3) = wrLocalnCells
+ count4( 4) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDtracers, start4, count4, tracers)
+
+
+ end subroutine write_netcdf_fields
+
+
+ subroutine write_netcdf_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(wr_ncid)
+
+ end subroutine write_netcdf_finalize
+
+end module write_netcdf
Added: branches/ocean_projects/grid_gen_ocean/quad_periodic/namelist.input
===================================================================
--- branches/ocean_projects/grid_gen_ocean/quad_periodic/namelist.input (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/quad_periodic/namelist.input 2010-03-23 19:17:10 UTC (rev 155)
@@ -0,0 +1,6 @@
+&periodic_grid
+ nVertLevels = 1,
+ nTracers = 2,
+ nVertexDegree = 4,
+ nproc = 2, 4, 8, 16, 32, 64,
+/
Added: branches/ocean_projects/grid_gen_ocean/quad_periodic/periodic_grid.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/quad_periodic/periodic_grid.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/quad_periodic/periodic_grid.F 2010-03-23 19:17:10 UTC (rev 155)
@@ -0,0 +1,453 @@
+program square_periodic_grid
+
+ use cell_indexing
+ use write_netcdf
+
+ implicit none
+
+! specify grid resolution with variable dc (meters)
+ real (kind=8), parameter :: dc = 50*1000.0 ! 50km distance between cell centers
+! real (kind=8), parameter :: dc = 10*1000.0 ! 10km distance between cell centers
+! real (kind=8), parameter :: dc = 5*1000.0 ! 5km distance between cell centers
+
+ real (kind=8), parameter :: &
+ lx = 2500*1000.0, & ! domain length in x, m
+ ly = 4300*1000.0, & ! domain length in x, m
+ h_total_max = 100.0, &
+ u_max = 0.0, &
+ u_src_max = 0.1, & ! max wind stress, N/m2
+ beta = 5.0e-11, &
+ f0 = 1.4e-4
+
+ real (kind=8), parameter :: pi = 3.141592653589793
+ real (kind=8), parameter :: ONE = 1.0_8
+ real (kind=8), parameter :: TWO = 2.0_8
+ real (kind=8), parameter :: FOUR = 4.0_8
+
+ integer, allocatable, dimension(:) :: indexToCellID, indexToEdgeID, indexToVertexID
+ integer, allocatable, dimension(:) :: nEdgesOnCell, nEdgesOnEdge
+ integer, allocatable, dimension(:,:) :: cellsOnCell, edgesOnCell, verticesOnCell
+ integer, allocatable, dimension(:,:) :: cellsOnEdge, edgesOnEdge, verticesOnEdge
+ integer, allocatable, dimension(:,:) :: edgesOnVertex, cellsOnVertex, uBC
+ real (kind=8), allocatable, dimension(:) :: areaTriangle, areaCell, angleEdge
+ real (kind=8), allocatable, dimension(:) :: dcEdge, dvEdge
+ 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, 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
+
+ integer :: i, j, np, iCell, iVertex, iEdge
+ integer :: nCells, nEdges, nVertices
+ integer :: iRow, iCol, ii, jj
+ integer :: nprocx, nprocy
+ real (kind=8) :: r
+ character (len=32) :: decomp_fname
+
+ real (kind=8), allocatable, dimension(:,:,:) :: utmp
+ integer :: iVert1, iEdge1, iEdge2
+ real (kind=8) :: ymid, ytmp, ymax
+
+ call cell_indexing_read_nl()
+
+ ! mrp 100305: replace nx and ny, specify using domain size.
+ nx = lx/dc;
+ ny = ly/dc;
+ ! mrp 100305: end
+
+ nCells = nx*ny
+ nEdges = 2*nCells
+ nVertices = nCells
+
+ allocate(indexToCellID(nCells))
+ allocate(indexToEdgeID(nEdges))
+ allocate(indexToVertexID(nVertices))
+
+ allocate(nEdgesOnCell(nCells))
+ allocate(cellsOnCell(maxEdges, nCells))
+ allocate(edgesOnCell(maxEdges, nCells))
+ allocate(verticesOnCell(maxEdges, nCells))
+
+ allocate(nEdgesOnEdge(nEdges))
+ allocate(cellsOnEdge(2,nEdges))
+ allocate(verticesOnEdge(2,nEdges))
+ allocate(edgesOnEdge(2*maxEdges,nEdges))
+ allocate(weightsOnEdge(2*maxEdges,nEdges))
+
+ allocate(edgesOnVertex(nVertexDegree,nVertices))
+ allocate(cellsOnVertex(nVertexDegree,nVertices))
+ allocate(kiteAreasOnVertex(nVertexDegree,nVertices))
+
+ allocate(areaTriangle(nVertices))
+ allocate(areaCell(nCells))
+
+ allocate(dcEdge(nEdges))
+ allocate(dvEdge(nEdges))
+ allocate(angleEdge(nEdges))
+
+ allocate(latCell(nCells))
+ allocate(lonCell(nCells))
+ allocate(xCell(nCells))
+ allocate(yCell(nCells))
+ allocate(zCell(nCells))
+ allocate(latEdge(nEdges))
+ allocate(lonEdge(nEdges))
+ allocate(xEdge(nEdges))
+ allocate(yEdge(nEdges))
+ allocate(zEdge(nEdges))
+ allocate(latVertex(nVertices))
+ allocate(lonVertex(nVertices))
+ allocate(xVertex(nVertices))
+ allocate(yVertex(nVertices))
+ allocate(zVertex(nVertices))
+
+ allocate(fEdge(nEdges))
+ 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))
+ allocate(vh(nVertLevels,nEdges,1))
+ allocate(h(nVertLevels,nCells,1))
+ allocate(circulation(nVertLevels,nVertices,1))
+ allocate(vorticity(nVertLevels,nVertices,1))
+ allocate(ke(nVertLevels,nCells,1))
+ allocate(tracers(nTracers,nVertLevels,nCells,1))
+
+ allocate(utmp(0:nx+1, 0:ny+1, 2))
+
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellIdx(iCol,iRow)
+ nEdgesOnCell(iCell) = nVertexDegree
+ do j=1,maxEdges
+ cellsOnCell(j,iCell) = cellOnCell(iCol,iRow,j)
+ edgesOnCell(j,iCell) = edgeOnCell(iCell,j)
+ verticesOnCell(j,iCell) = vertexOnCell(iCell,j)
+ end do
+ do j=1,2
+ cellsOnEdge(2,edgesOnCell(j,iCell)) = iCell ! Edges owned by this cell
+ end do
+ do j=3,4
+ cellsOnEdge(1,edgesOnCell(j,iCell)) = iCell
+ end do
+
+ verticesOnEdge(1,edgesOnCell(1,iCell)) = verticesOnCell(2,iCell) ! For edges owned by this cell
+ verticesOnEdge(2,edgesOnCell(1,iCell)) = verticesOnCell(1,iCell)
+ verticesOnEdge(1,edgesOnCell(2,iCell)) = verticesOnCell(3,iCell)
+ verticesOnEdge(2,edgesOnCell(2,iCell)) = verticesOnCell(2,iCell)
+
+ edgesOnEdge(1,edgesOnCell(3,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(2,edgesOnCell(3,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(3,edgesOnCell(3,iCell)) = edgesOnCell(2,iCell)
+
+ edgesOnEdge(1,edgesOnCell(4,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(2,edgesOnCell(4,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(3,edgesOnCell(4,iCell)) = edgesOnCell(3,iCell)
+
+ edgesOnEdge(4,edgesOnCell(1,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(5,edgesOnCell(1,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(6,edgesOnCell(1,iCell)) = edgesOnCell(4,iCell)
+
+ edgesOnEdge(4,edgesOnCell(2,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(5,edgesOnCell(2,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(6,edgesOnCell(2,iCell)) = edgesOnCell(1,iCell)
+
+ weightsOnEdge(1,edgesOnCell(3,iCell)) = ONE / FOUR
+ weightsOnEdge(2,edgesOnCell(3,iCell)) = 0.0
+ weightsOnEdge(3,edgesOnCell(3,iCell)) = ONE / FOUR
+
+ weightsOnEdge(1,edgesOnCell(4,iCell)) = -ONE / FOUR
+ weightsOnEdge(2,edgesOnCell(4,iCell)) = 0.0
+ weightsOnEdge(3,edgesOnCell(4,iCell)) = -ONE / FOUR
+
+ weightsOnEdge(4,edgesOnCell(1,iCell)) = ONE / FOUR
+ weightsOnEdge(5,edgesOnCell(1,iCell)) = 0.0
+ weightsOnEdge(6,edgesOnCell(1,iCell)) = ONE / FOUR
+
+ weightsOnEdge(4,edgesOnCell(2,iCell)) = -ONE / FOUR
+ weightsOnEdge(5,edgesOnCell(2,iCell)) = 0.0
+ weightsOnEdge(6,edgesOnCell(2,iCell)) = -ONE / FOUR
+
+ cellsOnVertex(1,verticesOnCell(1,iCell)) = iCell
+ cellsOnVertex(2,verticesOnCell(2,iCell)) = iCell
+ cellsOnVertex(3,verticesOnCell(3,iCell)) = iCell
+ cellsOnVertex(4,verticesOnCell(4,iCell)) = iCell
+
+ edgesOnVertex(3,verticesOnCell(1,iCell)) = edgesOnCell(1,iCell)
+ edgesOnVertex(1,verticesOnCell(2,iCell)) = edgesOnCell(1,iCell)
+ edgesOnVertex(4,verticesOnCell(2,iCell)) = edgesOnCell(2,iCell)
+ edgesOnVertex(2,verticesOnCell(3,iCell)) = edgesOnCell(2,iCell)
+
+ end do
+ end do
+
+ latCell = 0.0
+ lonCell = 0.0
+ latEdge = 0.0
+ lonEdge = 0.0
+ latVertex = 0.0
+ lonVertex = 0.0
+
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellIdx(iCol, iRow)
+ indexToCellID(iCell) = iCell
+ areaCell = dc*dc
+
+ xCell(iCell) = dc*real(iCol) - 0.5*dc
+ yCell(iCell) = dc*real(iRow) - 0.5*dc
+ zCell(iCell) = 0.0
+
+ xEdge(edgesOnCell(1,iCell)) = xCell(iCell) - 0.5*dc
+ yEdge(edgesOnCell(1,iCell)) = yCell(iCell)
+ zEdge(edgesOnCell(1,iCell)) = 0.0
+
+ xEdge(edgesOnCell(2,iCell)) = xCell(iCell)
+ yEdge(edgesOnCell(2,iCell)) = yCell(iCell) - 0.5*dc
+ zEdge(edgesOnCell(2,iCell)) = 0.0
+
+ xVertex(verticesOnCell(1,iCell)) = xCell(iCell) - 0.5*dc
+ yVertex(verticesOnCell(1,iCell)) = yCell(iCell) + 0.5*dc
+ zVertex(verticesOnCell(1,iCell)) = 0.0
+
+ angleEdge(edgesOnCell(1,iCell)) = pi / TWO
+ angleEdge(edgesOnCell(2,iCell)) = 0.0
+
+ end do
+ end do
+
+ write(6,*) 'maxval(xVertex), minval(xVertex)',maxval(xVertex), minval(xVertex)
+ write(6,*) 'maxval(yVertex), minval(yVertex)',maxval(yVertex), minval(yVertex)
+
+ do i=1,nEdges
+ indexToEdgeID(i) = i
+! mrp 100309 was:
+! nEdgesOnEdge(i) = 2*maxEdges
+! mrp 100309 change to:
+ nEdgesOnEdge(i) = 6
+! mrp 100309 end
+ dcEdge(i) = dc
+ dvEdge(i) = dc
+ end do
+
+ do i=1,nVertices
+ indexToVertexID(i) = i
+ areaTriangle(i) = dc*dc
+ do j=1,nVertexDegree
+ kiteAreasOnVertex(j,i) = dc*dc/FOUR
+ end do
+ end do
+
+
+ !
+ ! fill in initial conditions below
+ ! NOTE: these initial conditions will likely be removed
+ ! from the grid.nc files at some point (soon).
+ ! Initialize fields in grid
+ !
+
+ fEdge(:) = 0.0
+ fVertex(:) = 0.0
+ h_s(:) = 0.0
+ u(:,:,:) = u_max
+ v(:,:,:) = u_max
+ vh(:,:,:) = 0.0
+ circulation(:,:,:) = 0.0
+ vorticity(:,:,:) = 0.0
+ ke(:,:,:) = 0.0
+ tracers(:,:,:,:) = 0.0
+ h(:,:,:) = h_total_max/nVertLevels
+! do i=1,nCells
+! r = sqrt((xCell(i) - (nx/2)*dc)**2.0 + (yCell(i) - (ny/2)*dc)**2.0)
+! if (r < 10.0*dc) then
+! tracers(1,1,i,1) = (20.0 / 2.0) * (1.0 + cos(pi*r/(10.0*dc))) + 0.0
+! h(1,i,1) = 1.0 + 0.1*cos(pi*r/(20.0*dc))
+! else
+! tracers(1,1,i,1) = 0.0
+! h(1,i,1) = 1.0
+! end if
+! end do
+
+ ymid = dc*real(ny/2)
+
+ 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) - 0.5*dc
+ ytmp = dc*real(iRow) -0.5*dc
+ utmp(iCol, iRow, 1) = - u_src_max * cos ( TWO * pi * ytmp / ymax)
+ utmp(iCol, iRow, 2) = 0.0
+
+ !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=2*(iCell-1)+1
+ iEdge2=2*(iCell-1)+2
+ u_src(1,iEdge1) = utmp(iCol,iRow,1)
+ u_src(1,iEdge2) = utmp(iCol,iRow,2)
+ 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
+ !
+ call write_netcdf_init( nCells, nEdges, nVertices, maxEdges, nVertLevels, nTracers, nVertexDegree)
+
+ call write_netcdf_fields( 1, &
+ latCell, lonCell, xCell, yCell, zCell, indexToCellID, &
+ latEdge, lonEdge, xEdge, yEdge, zEdge, indexToEdgeID, &
+ latVertex, lonVertex, xVertex, yVertex, zVertex, indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex, &
+ h_s, &
+ uBC, &
+ u_src, &
+ u, &
+ v, &
+ h, &
+ vh, &
+ circulation, &
+ vorticity, &
+ ke, &
+ tracers &
+ )
+
+ call write_netcdf_finalize()
+
+ !
+ ! Write a graph.info file to be partitioned by kmetis
+ !
+ np = 1
+ do while (nproc(np) > 0)
+ call decompose_nproc(nproc(np), nprocx, nprocy)
+ if (nproc(np) < 10) then
+ write(decomp_fname,'(a,i1)') 'graph.info.part.',nproc(np)
+ else if (nproc(np) < 100) then
+ write(decomp_fname,'(a,i2)') 'graph.info.part.',nproc(np)
+ else if (nproc(np) < 1000) then
+ write(decomp_fname,'(a,i3)') 'graph.info.part.',nproc(np)
+ else if (nproc(np) < 10000) then
+ write(decomp_fname,'(a,i4)') 'graph.info.part.',nproc(np)
+ end if
+ indexToCellID(:) = -1
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellIdx(iCol, iRow)
+ ii = nprocx*real(iCol-1)/real(nx)
+ jj = nprocy*real(iRow-1)/real(ny)
+ indexToCellID(iCell) = jj*nprocx+ii
+ end do
+ end do
+ open(21,file=trim(decomp_fname),status='unknown')
+ do i=1,nCells
+ write(21,*) indexToCellID(i)
+ end do
+ close(21)
+ np = np + 1
+ end do
+
+end program square_periodic_grid
+
+
+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.01*dc) then
+ uBC(:,iEdge) = 1
+ u(:,iEdge,:) = 0.0
+ endif
+ 10 format(2f10.2)
+enddo
+
+end subroutine enforce_uBC
+
+
+subroutine decompose_nproc(nproc, nprocx, nprocy)
+
+ implicit none
+
+ integer, intent(in) :: nproc
+ integer, intent(out) :: nprocx, nprocy
+
+ do nprocx=int(sqrt(real(nproc))),1,-1
+ nprocy = nproc / nprocx
+ if (nprocy == ceiling(real(nproc)/real(nprocx))) return
+ end do
+
+end subroutine decompose_nproc
</font>
</pre>