<p><b>xylar@lanl.gov</b> 2012-01-20 13:43:50 -0700 (Fri, 20 Jan 2012)</p><p>BRANC COMMIT<br>
<br>
Adding a version of the periodic_hex grid generator that uses land ice variables<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice/icesheet/periodic_hex/Makefile
===================================================================
--- branches/land_ice/icesheet/periodic_hex/Makefile (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/Makefile 2012-01-20 20:43:50 UTC (rev 1402)
@@ -0,0 +1,68 @@
+# 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
+
+# gfortran
+FC = gfortran
+CC = gcc
+FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian
+CFLAGS = -O3 -m64
+LDFLAGS = -O3 -m64
+
+# 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 -lnetcdff
+
+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/land_ice/icesheet/periodic_hex/module_cell_indexing.F
===================================================================
--- branches/land_ice/icesheet/periodic_hex/module_cell_indexing.F (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/module_cell_indexing.F 2012-01-20 20:43:50 UTC (rev 1402)
@@ -0,0 +1,164 @@
+module cell_indexing
+
+! this subroutine provide index mapping for hexagon meshes dimensioned (nx, ny)
+
+ integer, parameter :: maxEdges = 6
+
+ integer :: nx, ny, nVertLevels, nTracers, vertexDegree
+ real (kind=8) :: dc
+ integer, dimension(20) :: nproc
+
+
+ contains
+
+
+ subroutine cell_indexing_read_nl()
+
+ implicit none
+
+ namelist /periodic_grid/ nx, ny, dc, nVertLevels, nTracers, nproc, vertexDegree
+
+ nx = 200
+ ny = 200
+ dc = 10000.
+ nVertLevels = 1
+ nTracers = 2
+ nproc(:) = -1
+ vertexDegree = 3
+
+ 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 (mod(iRow,2) == 1) then
+ if (neighborNumber == 1) then
+ cellOnCell = cellIdx(mx, iRow)
+ else if (neighborNumber == 2) then
+ cellOnCell = cellIdx(mx, my)
+ else if (neighborNumber == 3) then
+ cellOnCell = cellIdx(iCol, my)
+ else if (neighborNumber == 4) then
+ cellOnCell = cellIdx(px, iRow)
+ else if (neighborNumber == 5) then
+ cellOnCell = cellIdx(iCol, py)
+ else if (neighborNumber == 6) then
+ cellOnCell = cellIdx(mx, py)
+ end if
+ else
+ 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, my)
+ else if (neighborNumber == 4) then
+ cellOnCell = cellIdx(px, iRow)
+ else if (neighborNumber == 5) then
+ cellOnCell = cellIdx(px, py)
+ else if (neighborNumber == 6) then
+ cellOnCell = cellIdx(iCol, py)
+ end if
+ end if
+
+ 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 = 3*(iCell - 1) + 1
+ else if (neighborNumber == 2) then
+ edgeOnCell = 3*(iCell - 1) + 2
+ else if (neighborNumber == 3) then
+ edgeOnCell = 3*(iCell - 1) + 3
+ else if (neighborNumber == 4) then
+ edgeOnCell = 3*(cellOnCell(myCol, myRow, 4) - 1) + 1
+ else if (neighborNumber == 5) then
+ edgeOnCell = 3*(cellOnCell(myCol, myRow, 5) - 1) + 2
+ else if (neighborNumber == 6) then
+ edgeOnCell = 3*(cellOnCell(myCol, myRow, 6) - 1) + 3
+ end if
+
+ end function edgeOnCell
+
+
+ integer function vertexOnCell(iCell, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCell, neighborNumber
+
+ integer :: myRow, myCol
+
+ call cellColRow(iCell, myCol, myRow)
+
+ if (neighborNumber == 1) then
+ vertexOnCell = 2*(iCell - 1) + 1
+ else if (neighborNumber == 2) then
+ vertexOnCell = 2*(iCell - 1) + 2
+ else if (neighborNumber == 3) then
+ vertexOnCell = 2*(cellOnCell(myCol, myRow, 3) - 1) + 1
+ else if (neighborNumber == 4) then
+ vertexOnCell = 2*(cellOnCell(myCol, myRow, 4) - 1) + 2
+ else if (neighborNumber == 5) then
+ vertexOnCell = 2*(cellOnCell(myCol, myRow, 4) - 1) + 1
+ else if (neighborNumber == 6) then
+ vertexOnCell = 2*(cellOnCell(myCol, myRow, 5) - 1) + 2
+ end if
+
+ end function vertexOnCell
+
+
+end module cell_indexing
Added: branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F
===================================================================
--- branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F 2012-01-20 20:43:50 UTC (rev 1402)
@@ -0,0 +1,542 @@
+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 :: wrDimIDnVertLevelsPlus2
+ 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 :: wrVarIDthickness
+ integer :: wrVarIDbedTopography
+ integer :: wrVarIDbeta
+ integer :: wrVarIDnormalVelocity
+ integer :: wrVarIDtracers
+
+ integer :: wrLocalnCells
+ integer :: wrLocalnEdges
+ integer :: wrLocalnVertices
+ integer :: wrLocalmaxEdges
+ integer :: wrLocalnVertLevels
+ integer :: wrLocalnVertLevelsPlus2
+ integer :: wrLocalnTracers
+
+ contains
+
+ subroutine write_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ nVertLevels, &
+ nVertLevelsPlus2, &
+ nTracers, &
+ vertexDegree &
+ )
+
+ 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) :: nVertLevelsPlus2
+ integer, intent(in) :: nTracers
+ integer, intent(in) :: vertexDegree
+
+ integer :: nferr
+ integer, dimension(10) :: dimlist
+ character (len=16) :: on_a_sphere
+ real (kind=8) :: sphere_radius
+
+
+ wrLocalnCells = nCells
+ wrLocalnEdges = nEdges
+ wrLocalnVertices = nVertices
+ wrLocalmaxEdges = maxEdges
+ wrLocalnVertLevels = nVertLevels
+ wrLocalnTracers = nTracers
+
+ on_a_sphere = 'NO '
+ sphere_radius = 0.0
+
+ 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, 'vertexDegree', vertexDegree, wrDimIDvertexDegree)
+ nferr = nf_def_dim(wr_ncid, 'nVertLevels', nVertLevels, wrDimIDnVertLevels)
+ nferr = nf_def_dim(wr_ncid, 'nVertLevelsPlus2', nVertLevelsPlus2, wrDimIDnVertLevelsPlus2)
+ nferr = nf_def_dim(wr_ncid, 'nTracers', nTracers, wrDimIDnTracers)
+ nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
+
+
+ !
+ ! Define attributes
+ !
+ nferr = nf_put_att_text(wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, on_a_sphere)
+ nferr = nf_put_att_double(wr_ncid, NF_GLOBAL, 'sphere_radius', NF_DOUBLE, 1, sphere_radius)
+
+
+ !
+ ! 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) = wrDimIDnCells
+ dimlist( 2) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'thickness', NF_DOUBLE, 2, dimlist, wrVarIDthickness)
+ dimlist( 1) = wrDimIDnCells
+ dimlist( 2) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'bedTopography', NF_DOUBLE, 2, dimlist, wrVarIDbedTopography)
+ dimlist( 1) = wrDimIDnCells
+ dimlist( 2) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'beta', NF_DOUBLE, 2, dimlist, wrVarIDbeta)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'normalVelocity', NF_DOUBLE, 3, dimlist, wrVarIDnormalVelocity)
+ dimlist( 1) = wrDimIDnTracers
+ dimlist( 2) = wrDimIDnVertLevelsPlus2
+ 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, &
+ thickness, &
+ bedTopography, &
+ beta, &
+ normalVelocity, &
+ 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) :: thickness
+ real (kind=8), dimension(:,:), intent(in) :: bedTopography
+ real (kind=8), dimension(:,:), intent(in) :: beta
+ real (kind=8), dimension(:,:,:), intent(in) :: normalVelocity
+ 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) = 3
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = 3
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = 3
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+ start2(2) = time
+ count2( 2) = wrLocalnCells
+ count2( 2) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDthickness, start2, count2, thickness)
+
+ start2(2) = time
+ count2( 2) = wrLocalnCells
+ count2( 2) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDbedTopography, start2, count2, bedTopography)
+
+ start2(2) = time
+ count2( 2) = wrLocalnCells
+ count2( 2) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDbeta, start2, count2, beta)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDnormalVelocity, start3, count3, normalVelocity)
+
+ start4(4) = time
+ count4( 1) = wrLocalnTracers
+ count4( 2) = wrLocalnVertLevelsPlus2
+ 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/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km
===================================================================
--- branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km 2012-01-20 20:43:50 UTC (rev 1402)
@@ -0,0 +1,8 @@
+&periodic_grid
+ nx = 30,
+ ny = 30,
+ dc = 2000.,
+ nVertLevels = 9,
+ nTracers = 1,
+ nproc = 2,
+/
Property changes on: branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km
___________________________________________________________________
Added: svn:executable
+ *
Added: branches/land_ice/icesheet/periodic_hex/periodic_grid.F
===================================================================
--- branches/land_ice/icesheet/periodic_hex/periodic_grid.F (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/periodic_grid.F 2012-01-20 20:43:50 UTC (rev 1402)
@@ -0,0 +1,363 @@
+program hexagonal_periodic_grid
+
+ use cell_indexing
+ use write_netcdf
+
+ implicit none
+
+ 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 :: THREE = 3.0_8
+ real (kind=8), parameter :: FOUR = 4.0_8
+ real (kind=8), parameter :: SIX = 6.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
+ 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
+ real (kind=8), allocatable, dimension(:,:) :: thickness, bedTopography, beta
+ real (kind=8), allocatable, dimension(:,:,:) :: normalVelocity
+ real (kind=8), allocatable, dimension(:,:,:,:) :: tracers
+
+ integer :: i, j, np, iCell
+ integer :: nCells, nEdges, nVertices, nVertLevelsPlus2
+ integer :: iRow, iCol, ii, jj
+ integer :: nprocx, nprocy
+ real (kind=8) :: r
+ character (len=32) :: decomp_fname
+
+ call cell_indexing_read_nl()
+
+ nCells = nx*ny
+ nEdges = 3*nCells
+ nVertices = 2*nCells
+ nVertLevelsPlus2 = nVertLevels+2
+
+ 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(3,nVertices))
+ allocate(cellsOnVertex(3,nVertices))
+ allocate(kiteAreasOnVertex(3,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(thickness(nCells,1))
+ allocate(bedTopography(nCells,1))
+ allocate(beta(nCells,1))
+ allocate(normalVelocity(nVertLevels,nEdges,1))
+ allocate(tracers(nTracers,nVertLevelsPlus2,nCells,1))
+
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellIdx(iCol,iRow)
+ nEdgesOnCell(iCell) = 6
+ 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,3
+ cellsOnEdge(2,edgesOnCell(j,iCell)) = iCell ! Edges owned by this cell
+ end do
+ do j=4,6
+ 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)
+ verticesOnEdge(1,edgesOnCell(3,iCell)) = verticesOnCell(4,iCell)
+ verticesOnEdge(2,edgesOnCell(3,iCell)) = verticesOnCell(3,iCell)
+
+ edgesOnEdge(1,edgesOnCell(4,iCell)) = edgesOnCell(5,iCell)
+ edgesOnEdge(2,edgesOnCell(4,iCell)) = edgesOnCell(6,iCell)
+ edgesOnEdge(3,edgesOnCell(4,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(4,edgesOnCell(4,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(5,edgesOnCell(4,iCell)) = edgesOnCell(3,iCell)
+
+ edgesOnEdge(1,edgesOnCell(5,iCell)) = edgesOnCell(6,iCell)
+ edgesOnEdge(2,edgesOnCell(5,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(3,edgesOnCell(5,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(4,edgesOnCell(5,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(5,edgesOnCell(5,iCell)) = edgesOnCell(4,iCell)
+
+ edgesOnEdge(1,edgesOnCell(6,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(2,edgesOnCell(6,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(3,edgesOnCell(6,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(4,edgesOnCell(6,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(5,edgesOnCell(6,iCell)) = edgesOnCell(5,iCell)
+
+ edgesOnEdge(6,edgesOnCell(1,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(7,edgesOnCell(1,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(8,edgesOnCell(1,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(9,edgesOnCell(1,iCell)) = edgesOnCell(5,iCell)
+ edgesOnEdge(10,edgesOnCell(1,iCell)) = edgesOnCell(6,iCell)
+
+ edgesOnEdge(6,edgesOnCell(2,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(7,edgesOnCell(2,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(8,edgesOnCell(2,iCell)) = edgesOnCell(5,iCell)
+ edgesOnEdge(9,edgesOnCell(2,iCell)) = edgesOnCell(6,iCell)
+ edgesOnEdge(10,edgesOnCell(2,iCell)) = edgesOnCell(1,iCell)
+
+ edgesOnEdge(6,edgesOnCell(3,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(7,edgesOnCell(3,iCell)) = edgesOnCell(5,iCell)
+ edgesOnEdge(8,edgesOnCell(3,iCell)) = edgesOnCell(6,iCell)
+ edgesOnEdge(9,edgesOnCell(3,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(10,edgesOnCell(3,iCell)) = edgesOnCell(2,iCell)
+
+ weightsOnEdge(1,edgesOnCell(4,iCell)) = ONE / THREE
+ weightsOnEdge(2,edgesOnCell(4,iCell)) = ONE / SIX
+ weightsOnEdge(3,edgesOnCell(4,iCell)) = 0.0
+ weightsOnEdge(4,edgesOnCell(4,iCell)) = ONE / SIX
+ weightsOnEdge(5,edgesOnCell(4,iCell)) = ONE / THREE
+
+ weightsOnEdge(1,edgesOnCell(5,iCell)) = ONE / THREE
+ weightsOnEdge(2,edgesOnCell(5,iCell)) = -ONE / SIX
+ weightsOnEdge(3,edgesOnCell(5,iCell)) = 0.0
+ weightsOnEdge(4,edgesOnCell(5,iCell)) = ONE / SIX
+ weightsOnEdge(5,edgesOnCell(5,iCell)) = -ONE / THREE
+
+ weightsOnEdge(1,edgesOnCell(6,iCell)) = -ONE / THREE
+ weightsOnEdge(2,edgesOnCell(6,iCell)) = -ONE / SIX
+ weightsOnEdge(3,edgesOnCell(6,iCell)) = 0.0
+ weightsOnEdge(4,edgesOnCell(6,iCell)) = -ONE / SIX
+ weightsOnEdge(5,edgesOnCell(6,iCell)) = -ONE / THREE
+
+ weightsOnEdge(6,edgesOnCell(1,iCell)) = ONE / THREE
+ weightsOnEdge(7,edgesOnCell(1,iCell)) = ONE / SIX
+ weightsOnEdge(8,edgesOnCell(1,iCell)) = 0.0
+ weightsOnEdge(9,edgesOnCell(1,iCell)) = ONE / SIX
+ weightsOnEdge(10,edgesOnCell(1,iCell)) = ONE / THREE
+
+ weightsOnEdge(6,edgesOnCell(2,iCell)) = ONE / THREE
+ weightsOnEdge(7,edgesOnCell(2,iCell)) = -ONE / SIX
+ weightsOnEdge(8,edgesOnCell(2,iCell)) = 0.0
+ weightsOnEdge(9,edgesOnCell(2,iCell)) = ONE / SIX
+ weightsOnEdge(10,edgesOnCell(2,iCell)) = -ONE / THREE
+
+ weightsOnEdge(6,edgesOnCell(3,iCell)) = -ONE / THREE
+ weightsOnEdge(7,edgesOnCell(3,iCell)) = -ONE / SIX
+ weightsOnEdge(8,edgesOnCell(3,iCell)) = 0.0
+ weightsOnEdge(9,edgesOnCell(3,iCell)) = -ONE / SIX
+ weightsOnEdge(10,edgesOnCell(3,iCell)) = -ONE / THREE
+
+ cellsOnVertex(3,verticesOnCell(2,iCell)) = iCell
+ cellsOnVertex(1,verticesOnCell(4,iCell)) = iCell
+ cellsOnVertex(2,verticesOnCell(6,iCell)) = iCell
+ cellsOnVertex(1,verticesOnCell(1,iCell)) = iCell
+ cellsOnVertex(2,verticesOnCell(3,iCell)) = iCell
+ cellsOnVertex(3,verticesOnCell(5,iCell)) = iCell
+
+ edgesOnVertex(1,verticesOnCell(1,iCell)) = edgesOnCell(1,iCell)
+ edgesOnVertex(1,verticesOnCell(2,iCell)) = edgesOnCell(1,iCell)
+ edgesOnVertex(3,verticesOnCell(3,iCell)) = edgesOnCell(2,iCell)
+ edgesOnVertex(3,verticesOnCell(2,iCell)) = edgesOnCell(2,iCell)
+ edgesOnVertex(2,verticesOnCell(3,iCell)) = edgesOnCell(3,iCell)
+ edgesOnVertex(2,verticesOnCell(4,iCell)) = edgesOnCell(3,iCell)
+ end do
+ end do
+
+ weightsOnEdge(:,:) = weightsOnEdge(:,:) * ONE/sqrt(THREE)
+
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellIdx(iCol, iRow)
+ indexToCellID(iCell) = iCell
+ areaCell = dc*dc*sqrt(THREE) / TWO
+ latCell(iCell) = 0.0
+ lonCell(iCell) = 0.0
+
+ if (mod(iRow,2) == 1) then
+ xCell(iCell) = dc*real(iCol) - 0.5*dc
+ yCell(iCell) = dc*real(iRow)*sqrt(THREE) / TWO
+ zCell(iCell) = 0.0
+ else
+ xCell(iCell) = dc*real(iCol)
+ yCell(iCell) = dc*real(iRow)*sqrt(THREE) / TWO
+ zCell(iCell) = 0.0
+ end if
+
+ latEdge(iCell) = 0.0
+ lonEdge(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) - 0.5 * dc * cos(pi/THREE)
+ yEdge(edgesOnCell(2,iCell)) = yCell(iCell) - 0.5 * dc * sin(pi/THREE)
+ zEdge(edgesOnCell(2,iCell)) = 0.0
+
+ xEdge(edgesOnCell(3,iCell)) = xCell(iCell) + 0.5 * dc * cos(pi/THREE)
+ yEdge(edgesOnCell(3,iCell)) = yCell(iCell) - 0.5 * dc * sin(pi/THREE)
+ zEdge(edgesOnCell(3,iCell)) = 0.0
+
+ latVertex(iCell) = 0.0
+ lonVertex(iCell) = 0.0
+ xVertex(verticesOnCell(1,iCell)) = xCell(iCell) - 0.5*dc
+ yVertex(verticesOnCell(1,iCell)) = yCell(iCell) + dc * sqrt(THREE) / SIX
+ zVertex(verticesOnCell(1,iCell)) = 0.0
+
+ xVertex(verticesOnCell(2,iCell)) = xCell(iCell) - 0.5*dc
+ yVertex(verticesOnCell(2,iCell)) = yCell(iCell) - dc * sqrt(THREE) / SIX
+ zVertex(verticesOnCell(2,iCell)) = 0.0
+
+ angleEdge(edgesOnCell(1,iCell)) = 0.0
+ angleEdge(edgesOnCell(2,iCell)) = pi / THREE
+ angleEdge(edgesOnCell(3,iCell)) = TWO*pi/THREE
+ end do
+ end do
+
+ do i=1,nEdges
+ indexToEdgeID(i) = i
+ nEdgesOnEdge(i) = 10
+ dcEdge(i) = dc
+ dvEdge(i) = dcEdge(i) * sqrt(THREE)/THREE
+ end do
+
+ do i=1,nVertices
+ indexToVertexID(i) = i
+ areaTriangle(i) = dc*dc*sqrt(THREE)/FOUR
+ do j=1,3
+ kiteAreasOnVertex(j,i) = dc*dc*sqrt(THREE)/(TWO*SIX)
+ 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
+ !
+
+ normalVelocity(:,:,:) = 0.0
+ tracers(:,:,:,:) = 1.0
+ thickness(:,:) = 1.0
+ bedTopography(:,:) = 0.0
+
+
+ !
+ ! Write grid to grid.nc file
+ !
+ call write_netcdf_init( nCells, nEdges, nVertices, maxEdges, nVertLevels, nVertLevelsPlus2, nTracers, vertexDegree )
+
+ 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, &
+ thickness, &
+ bedTopography, &
+ beta, &
+ normalVelocity, &
+ 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 hexagonal_periodic_grid
+
+
+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>