<p><b>mhoffman@lanl.gov</b> 2013-02-28 12:12:01 -0700 (Thu, 28 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Created a modified version of the ocean core's tool 'basin' to use on land ice meshes for removing non-ice cells.<br>
<br>
To use, first create a land ice mesh (e.g. with periodic_hex + add_land_ice_variables_to_mpas_grid.py) and then populate it with thickness data (e.g. with copy_CISM_grid_to_MPAS_grid.py). Then run define_keepCellMask.py which will add a new field call keepCellMask to the file. There are options in that script to add a buffer around the ice-covered portion of the domain. Next compile and run crop_landice_grid. This will return a cropped domain that is an empty grid, so the steps for converting it to a land ice grid and populating it with data will need to be repeated. You will also want to run pmetis manually on the graph.info file that is generated.<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/grid_tools/crop_landice_grid/Makefile
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/Makefile (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/Makefile 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,69 @@
+# 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 #-g -traceback -check all
+#CFLAGS = #-g
+#LDFLAGS = #-g -traceback -check all
+
+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
+#NETCDF = /Users/maltrud/local
+
+
+CPP = cpp -C -P -traditional
+CPPFLAGS =
+CPPINCLUDES =
+INCLUDES = -I$(NETCDF)/include
+LIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdff
+
+RM = rm -f
+
+##########################
+
+.SUFFIXES: .F .o
+
+OBJS = utilities.o \
+ module_read_netcdf.o \
+ module_cullLoops.o \
+ module_write_netcdf.o \
+ crop_landice_grid.o
+
+all: crop_landice_grid
+
+crop_landice_grid.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_cullLoops.o
+
+crop_landice_grid: $(OBJS)
+ $(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
+
+clean:
+ $(RM) *.o *.mod pop *.f90
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(INCLUDES)
+ #$(RM) $*.f90
Added: branches/land_ice_projects/grid_tools/crop_landice_grid/crop_landice_grid.F
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/crop_landice_grid.F (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/crop_landice_grid.F 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,1133 @@
+program crop_landice_grid
+
+use read_netcdf
+use write_netcdf
+use utilities
+use cullLoops
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Program: crop_landice_grid.F
+!
+! This program culls cells from a mesh. Presumably you will use it to remove non-ice cells.
+! You define what cells to keep with a 1 value in an integer field called keepCellMask that you must add to a grid file named land_ice_grid.nc.
+! It generates a culled version called land_ice_grid_culled.nc.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+implicit none
+
+! original grid variables
+integer :: time, nCells, nEdges, nVertices
+integer :: maxEdges, maxEdges2, TWO, vertexDegree, nVertLevels
+integer, allocatable, dimension(:) :: indexToCellID, indexToEdgeID, indexToVertexID
+real, allocatable, dimension(:) :: xCell, yCell, zCell, latCell, lonCell, meshDensity
+real, allocatable, dimension(:) :: xEdge, yEdge, zEdge, latEdge, lonEdge
+real, allocatable, dimension(:) :: xVertex, yVertex, zVertex, latVertex, lonVertex
+integer, allocatable, dimension(:) :: nEdgesOnCell, nEdgesOnEdge
+integer, allocatable, dimension(:,:) :: cellsOnCell, edgesOnCell, verticesOnCell
+integer, allocatable, dimension(:,:) :: cellsOnEdge, verticesOnEdge, edgesOnEdge
+integer, allocatable, dimension(:,:) :: cellsOnVertex, edgesOnVertex
+real, allocatable, dimension(:) :: areaCell, areaTriangle, dcEdge, dvEdge, angleEdge
+real, allocatable, dimension(:,:) :: kiteAreasOnVertex, weightsOnEdge
+
+integer nlon, nlat, ndepth
+real(kind=4), allocatable, dimension(:) :: t_lon, t_lat
+
+real, dimension(40) :: dz
+
+ real (kind=8) :: ymid, ytmp, ymax, xmid, xloc, yloc, pert, ymin, distance, r, c1(3), c2(3)
+ real (kind=8) :: latmid, lattmp, latmax, latmin
+ integer :: cell1, cell2
+
+! new grid variables
+integer :: nCellsNew, nEdgesNew, nVerticesNew
+integer :: maxEdgesNew, maxEdges2New, TWONew, vertexDegreeNew, nVertLevelsNew
+integer, allocatable, dimension(:) :: indexToCellIDNew, indexToEdgeIDNew, indexToVertexIDNew
+real, allocatable, dimension(:) :: xCellNew, yCellNew, zCellNew, latCellNew, lonCellNew, meshDensityNew, meshSpacingNew
+real, allocatable, dimension(:) :: xEdgeNew, yEdgeNew, zEdgeNew, latEdgeNew, lonEdgeNew
+real, allocatable, dimension(:) :: xVertexNew, yVertexNew, zVertexNew, latVertexNew, lonVertexNew
+integer, allocatable, dimension(:) :: nEdgesOnCellNew, nEdgesOnEdgeNew, flipVerticesOnEdgeOrdering
+integer, allocatable, dimension(:,:) :: cellsOnCellNew, edgesOnCellNew, verticesOnCellNew
+integer, allocatable, dimension(:,:) :: cellsOnEdgeNew, verticesOnEdgeNew, edgesOnEdgeNew
+integer, allocatable, dimension(:,:) :: cellsOnVertexNew, edgesOnVertexNew
+integer, allocatable, dimension(:,:) :: boundaryEdgeNew, boundaryVertexNew
+real, allocatable, dimension(:) :: areaCellNew, areaTriangleNew, dcEdgeNew, dvEdgeNew, angleEdgeNew
+real, allocatable, dimension(:,:) :: kiteAreasOnVertexNew, weightsOnEdgeNew, normalsNew
+
+integer, allocatable, dimension(:) :: keepCellMask
+
+! mapping variables
+integer, allocatable, dimension(:) :: kmt
+integer, allocatable, dimension(:) :: cellMap, edgeMap, vertexMap
+
+! work variables
+integer :: i,j,jNew,k,jEdge,jEdgeNew,iVertex1New,iVertex2New,iCell1New,iCell2New
+integer :: iCell, iCell1, iCell2, iCell3, iEdge, iVertex, iVertex1, iVertex2
+integer :: iCellNew, iEdgeNew, iVertexNew, ndata, jCell1, jCell2, jCell, iter
+real :: xin, yin, zin, ulon, ulat, ux, uy, uz, rlon, rlat, temp_t, temp_s
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Namelist variables
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Variables in namelist file
+character (len=32) :: on_a_sphere, zLevel_thickness,bottom_topography, initial_conditions
+logical :: expand_from_unit_sphere, eliminate_inland_seas, &
+ check_mesh, &
+ cut_domain_from_sphere, solid_boundary_in_y, solid_boundary_in_x
+
+integer :: nVertLevelsMOD
+real (kind=8) :: sphere_radius
+
+!! specify namelist
+!namelist /landicecrop/ on_a_sphere, sphere_radius, &
+! expand_from_unit_sphere, &
+! eliminate_inland_seas, check_mesh, &
+! cut_domain_from_sphere, solid_boundary_in_y, solid_boundary_in_x
+!! Lx
+
+! Default namelist values. Default set for realistic global IC.
+on_a_sphere = 'NO'
+sphere_radius = 6.37122e6
+expand_from_unit_sphere = .true.
+
+eliminate_inland_seas=.false.
+solid_boundary_in_y = .false.
+solid_boundary_in_x = .false.
+
+! This needs to be changed for correct periodic boundaries
+! Lx is the TOTAL domain width, and needs to be exact for correct periodic
+! boundaries in x.
+!Lx = 3200.0e3 ! 40x80km=3200km
+
+
+!! Read in namelist
+! open(20,file='namelist.landicecrop',status='old')
+! read(20,landicecrop)
+! close(20)
+
+
+! get to work
+write(6,*) ' starting'
+write(6,*)
+
+! get grid
+write(6,*) ' calling read_grid'
+write(6,*)
+call read_grid
+write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
+
+nVertLevelsMOD = nVertLevels
+
+! copy dimensions
+write(6,*) ' copy dimensions'
+write(6,*)
+call copy_dimensions
+write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
+
+! define the kmt array
+write(6,*) ' calling define_kmt'
+write(6,*)
+call define_kmt
+
+! define the mapping between original and new cells, edges and vertices
+write(6,*) ' calling define_mapping'
+write(6,*)
+call define_mapping
+
+! copy the vector arrays form the original to new arrays
+write(6,*) ' calling map_vectors'
+write(6,*)
+call map_vectors
+
+! define the new connectivity variables
+write(6,*) ' calling map_connectivity'
+write(6,*)
+call map_connectivity
+
+! check the mesh
+if (check_mesh) then
+ call error_checking
+endif
+
+! dump new grid to netCDF
+write(6,*) ' calling write_grid'
+write(6,*)
+call write_grid
+
+! dump graph for partioning
+write(6,*) ' call write_graph'
+write(6,*)
+call write_graph
+
+!! write OpenDx file
+!if (write_OpenDX_flag) then
+! write(6,*) ' calling write_OpenDX'
+! write(6,*)
+! call write_OpenDX( on_a_sphere, &
+! nCellsNew, &
+! nVerticesNew, &
+! nEdgesNew, &
+! vertexDegreeNew, &
+! maxEdgesNew, &
+! xCellNew, &
+! yCellNew, &
+! zCellNew, &
+! xVertexNew, &
+! yVertexNew, &
+! zVertexNew, &
+! xEdgeNew, &
+! yEdgeNew, &
+! zEdgeNew, &
+! nEdgesOnCellNew, &
+! verticesOnCellNew, &
+! verticesOnEdgeNew, &
+! cellsOnVertexNew, &
+! edgesOnCellNew, &
+! areaCellNew, &
+! maxLevelCellNew, &
+! meshDensityNew, &
+! bottomDepthNew, &
+! temperatureNew(1,1,:), &
+! kiteAreasOnVertexNew )
+!endif
+
+!!do iCell=1,nCellsNew
+! !ulon = 1.0; ulat = 0.0
+! !xin = xCellNew(iCell); yin = yCellNew(iCell); zin = zCellNew(iCell)
+! !call transform_from_lonlat_to_xyz(xin, yin, zin, ulon, ulat, ux, uy, uz)
+! !if(abs(ux).lt.1.0e-10) ux=0.0
+! !if(abs(uy).lt.1.0e-10) uy=0.0
+! !if(abs(uz).lt.1.0e-10) uz=0.0
+! !write(20,10) ux, uy, uz
+! !10 format(3e25.10)
+!!enddo
+
+write(6,*) ' finished'
+
+contains
+
+subroutine write_graph
+implicit none
+integer :: m,itmp(maxEdgesNew),k
+
+ m=nEdgesNew
+ do i=1,nCellsNew
+ do j=1,nEdgesOnCellNew(i)
+ if(cellsOnCellNew(j,i).eq.0) m=m-1
+ enddo
+ enddo
+
+ open(42,file='graph.info',form='formatted')
+ write(42,*) nCellsNew, m
+ do i=1,nCellsNew
+ itmp = 0; k = 0;
+ do j=1,nEdgesOnCellNew(i)
+ if(cellsOnCellNew(j,i).gt.0) then
+ k=k+1; itmp(k)=cellsOnCellNew(j,i)
+ endif
+ enddo
+ write(42,'(1x,12i8)',advance='no') (itmp(m),m=1,k)
+ write(42,'(1x)')
+ end do
+ close(42)
+end subroutine write_graph
+
+
+subroutine error_checking
+real :: p(3), q(3), r(3), angle, s(3), t(3), dot, mindot, maxdot, b(vertexDegree)
+real :: work(nCellsNew)
+
+
+! write
+write(6,*)
+write(6,*) ' error checking '
+write(6,*)
+
+! check to see if every edge is normal to associated cells
+mindot = 2
+maxdot = -2
+do iEdge=1,nEdgesNew
+ if(boundaryEdgeNew(1,iEdge).eq.1) cycle
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ p(1)=xCellNew(iCell1); p(2)=yCellNew(iCell1); p(3)=zCellNew(iCell1)
+ q(1)=xCellNew(iCell2); q(2)=yCellNew(iCell2); q(3)=zCellNew(iCell2)
+ r(1)=xEdgeNew(iEdge); r(2)=yEdgeNew(iEdge); r(3)=zEdgeNew(iEdge)
+ call unit_vector_in_3space(p)
+ call unit_vector_in_3space(q)
+ call unit_vector_in_3space(r)
+ t = q - p
+ s = r - p
+ call unit_vector_in_3space(t)
+ call unit_vector_in_3space(s)
+ dot = s(1)*t(1)+s(2)*t(2)+s(3)*t(3)
+ if(dot.lt.mindot) mindot=dot
+ if(dot.gt.maxdot) maxdot=dot
+enddo
+write(6,10) 'alignment of edges and cells (should be ones)', mindot, maxdot
+10 format(a60,5x,2e15.5)
+
+! check to see if every segments connecting cells and vertices are orothogonal'
+mindot = 2
+maxdot = -2
+do iEdge=1,nEdgesNew
+ if(boundaryEdgeNew(1,iEdge).eq.1) cycle
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ iVertex1 = verticesOnEdgeNew(1,iEdge)
+ iVertex2 = verticesOnEdgeNew(2,iEdge)
+ p(1)=xCellNew(iCell1); p(2)=yCellNew(iCell1); p(3)=zCellNew(iCell1)
+ q(1)=xCellNew(iCell2); q(2)=yCellNew(iCell2); q(3)=zCellNew(iCell2)
+ r(1)=xVertexNew(iVertex1); r(2)=yVertexNew(iVertex1); r(3)=zVertexNew(iVertex1)
+ s(1)=xVertexNew(iVertex2); s(2)=yVertexNew(iVertex2); s(3)=zVertexNew(iVertex2)
+ call unit_vector_in_3space(p)
+ call unit_vector_in_3space(q)
+ call unit_vector_in_3space(r)
+ call unit_vector_in_3space(s)
+ t = q - p
+ s = s - r
+ call unit_vector_in_3space(t)
+ call unit_vector_in_3space(s)
+ dot = s(1)*t(1)+s(2)*t(2)+s(3)*t(3)
+ if(dot.lt.mindot) mindot=dot
+ if(dot.gt.maxdot) maxdot=dot
+enddo
+write(6,10) 'orthogonality of cell and vertex edges (should be zeros)', mindot, maxdot
+
+! check that the kiteareas sum to the areatriangle
+mindot = 2
+maxdot = -2
+do iVertex=1,nVerticesNew
+ b = 0
+ do i=1,vertexDegree
+ b(i) = kiteAreasOnVertexNew(i,iVertex)
+ enddo
+ angle = sum(b)
+ if(angle - areaTriangleNew(iVertex).lt.mindot) mindot = angle - areaTriangleNew(iVertex)
+ if(angle - areaTriangleNew(iVertex).gt.maxdot) maxdot = angle - areaTriangleNew(iVertex)
+enddo
+write(6,10) ' error in sum of kites and triangles (should be zeroes)', mindot, maxdot
+
+! check that the kiteareas sum to the areaCell
+mindot = 2
+maxdot = -2
+work = 0
+do iVertex=1,nVerticesNew
+ iCell1 = cellsOnVertexNew(1,iVertex)
+ iCell2 = cellsOnVertexNew(2,iVertex)
+ iCell3 = cellsOnVertexNew(3,iVertex)
+ if(iCell1.ne.0) work(iCell1) = work(iCell1) + kiteAreasOnVertexNew(1,iVertex)
+ if(iCell2.ne.0) work(iCell2) = work(iCell2) + kiteAreasOnVertexNew(2,iVertex)
+ if(iCell3.ne.0) work(iCell3) = work(iCell3) + kiteAreasOnVertexNew(3,iVertex)
+enddo
+mindot = minval(areaCellNew - work)
+maxdot = maxval(areaCellNew - work)
+write(6,10) ' error in sum of kites and cells (should be zeroes)', mindot, maxdot
+
+!check for connectivity inverses for cells/edges
+do iCell=1,nCellsNew
+ do i=1,nEdgesOnCellNew(iCell)
+ iEdge=edgesOnCellNew(i,iCell)
+ if(iEdge.le.0) stop ' iEdge le 0'
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ if(iCell1.ne.iCell.and.iCell2.ne.iCell) stop ' cells/edges inverse failed'
+ enddo
+enddo
+write(6,*) ' cellsOnEdge and edgesOnCell are duals for every cell/edge combination'
+
+!check for connectivity inverses for cells/vertices
+do iCell=1,nCellsNew
+ do i=1,nEdgesOnCellNew(iCell)
+ iVertex = verticesOnCellNew(i,iCell)
+ if(iVertex.le.0) stop ' iVertex le 0'
+ iCell1 = cellsOnVertexNew(1,iVertex)
+ iCell2 = cellsOnVertexNew(2,iVertex)
+ iCell3 = cellsOnVertexNew(3,iVertex)
+ if(iCell1.ne.iCell.and.iCell2.ne.iCell.and.iCell3.ne.iCell) stop ' cells/vertices inverse failed'
+ enddo
+enddo
+write(6,*) ' cellsOnVertex and verticesOnCell are duals for every cell/vertex combination'
+
+!check edgesOnEdge
+do iEdge=1,nEdgesNew
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ if(nEdgesOnEdgeNew(iEdge).eq.0) then
+ if(boundaryEdgeNew(1,iEdge).ne.1) stop ' stopping boundaryEdgeNew'
+ endif
+ do i=1,nEdgesOnEdgeNew(iEdge)
+ jEdge = edgesOnEdgeNew(i,iEdge)
+ jCell1 = cellsOnEdgeNew(1,jEdge)
+ jCell2 = cellsOnEdgeNew(2,jEdge)
+ if(jCell1.ne.iCell1.and.jCell1.ne.iCell2) then
+ if(jCell2.ne.iCell1.and.jCell2.ne.iCell2) then
+ write(6,*) 'error in edgesOnEdge'
+ write(6,*) iCell1, iCell2, jCell1, jCell2
+ stop
+ endif
+ endif
+ enddo
+enddo
+write(6,*) ' edgesOnEdge is consistent with cellsOnEdge'
+
+end subroutine error_checking
+
+
+subroutine copy_dimensions
+
+maxEdgesNew = maxEdges
+maxEdges2New = maxEdges2
+TWONew = TWO
+vertexDegreeNew = vertexDegree
+nVertLevelsNew = nVertLevelsMod
+
+write(6,*)
+write(6,*) ' new dimensions '
+write(6,*) ' maxEdgesNew : ', maxEdgesNew
+write(6,*) ' maxEdges2New : ', maxEdges2New
+write(6,*) ' TWONew : ', TWONew
+write(6,*) ' vertexDegreeNew : ', vertexDegreeNew
+write(6,*) ' nVertLevelsNew : ', nVertLevelsNew
+
+end subroutine copy_dimensions
+
+
+
+subroutine read_grid
+implicit none
+
+call read_netcdf_init(nCells, nEdges, nVertices, maxEdges,maxEdges2,&
+ nVertLevels,TWO,vertexDegree)
+
+write(6,*) ' init from grid '
+write(6,*) 'nCells :', nCells
+write(6,*) 'nEdges :', nEdges
+write(6,*) 'nVertices :', nVertices
+write(6,*) 'maxEdges :', maxEdges
+write(6,*) 'maxEdges2 :', maxEdges2
+write(6,*) 'nVertLevels :', nVertLevels
+write(6,*) 'vertexDegree :', vertexDegree
+write(6,*) 'TWO :', TWO
+
+allocate(xCell(nCells))
+allocate(yCell(nCells))
+allocate(zCell(nCells))
+allocate(latCell(nCells))
+allocate(lonCell(nCells))
+allocate(meshDensity(nCells))
+allocate(xEdge(nEdges))
+allocate(yEdge(nEdges))
+allocate(zEdge(nEdges))
+allocate(latEdge(nEdges))
+allocate(lonEdge(nEdges))
+allocate(xVertex(nVertices))
+allocate(yVertex(nVertices))
+allocate(zVertex(nVertices))
+allocate(latVertex(nVertices))
+allocate(lonVertex(nVertices))
+allocate(dcEdge(nEdges))
+allocate(dvEdge(nEdges))
+
+allocate(indexToCellID(nCells))
+allocate(indexToEdgeID(nEdges))
+allocate(indexToVertexID(nVertices))
+
+allocate(cellsOnEdge(TWO,nEdges))
+allocate(nEdgesOnCell(nCells))
+allocate(nEdgesOnEdge(nEdges))
+allocate(edgesOnCell(maxEdges,nCells))
+allocate(edgesOnEdge(maxEdges2,nEdges))
+allocate(weightsOnEdge(maxEdges2,nEdges))
+
+allocate(angleEdge(nEdges))
+allocate(areaCell(nCells))
+allocate(areaTriangle(nVertices))
+allocate(cellsOnCell(maxEdges,nCells))
+allocate(verticesOnCell(maxEdges,nCells))
+allocate(verticesOnEdge(TWO,nEdges))
+allocate(edgesOnVertex(vertexDegree,nVertices))
+allocate(cellsOnVertex(vertexDegree,nVertices))
+allocate(kiteAreasOnVertex(vertexDegree,nVertices))
+
+allocate(keepCellMask(nCells))
+
+
+xCell=0; yCell=0; zCell=0; latCell=0; lonCell=0; meshDensity=1.0
+xEdge=0; yEdge=0; zEdge=0; latEdge=0; lonEdge=0
+xVertex=0; yVertex=0; zVertex=0; latVertex=0; lonVertex=0
+
+indexToCellID=0; indexToEdgeID=0; indexToVertexID=0
+cellsOnEdge=0; nEdgesOnCell=0; edgesOnCell=0
+edgesOnEdge=0; weightsOnEdge=0
+angleEdge=0; areaCell=0; areaTriangle=0
+cellsOnCell=0; verticesOnCell=0; verticesOnEdge=0
+edgesOnVertex=0; cellsOnVertex=0; kiteAreasOnVertex=0
+
+
+call read_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ meshDensity, &
+ 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, &
+ keepCellMask &
+ )
+
+write(6,*) ' values from read grid, min/max'
+write(6,*) ' latCell : ', minval(latCell), maxval(latCell)
+write(6,*) ' lonCell : ', minval(lonCell), maxval(lonCell)
+write(6,*) ' meshDensity : ', minval(meshDensity),maxval(meshDensity)
+write(6,*) ' xCell : ', minval(xCell), maxval(xCell)
+write(6,*) ' yCell : ', minval(yCell), maxval(yCell)
+write(6,*) ' zCell : ', minval(zCell), maxval(zCell)
+write(6,*) ' indexToCellID : ', minval(indexToCellID), maxval(indexToCellID)
+write(6,*) ' latEdge : ', minval(latEdge), maxval(latEdge)
+write(6,*) ' lonEdge : ', minval(lonEdge), maxval(lonEdge)
+write(6,*) ' xEdge : ', minval(xEdge), maxval(xEdge)
+write(6,*) ' yEdge : ', minval(yEdge), maxval(yEdge)
+write(6,*) ' zEdge : ', minval(zEdge), maxval(zEdge)
+write(6,*) ' indexToEdgeID : ', minval(indexToEdgeID), maxval(indexToEdgeID)
+write(6,*) ' latVertex : ', minval(latVertex), maxval(latVertex)
+write(6,*) ' lonVertex : ', minval(lonVertex), maxval(lonVertex)
+write(6,*) ' xVertex : ', minval(xVertex), maxval(xVertex)
+write(6,*) ' yVertex : ', minval(yVertex), maxval(yVertex)
+write(6,*) ' zVertex : ', minval(zVertex), maxval(zVertex)
+write(6,*) ' indexToVertexID : ', minval(indexToVertexID), maxval(indexToVertexID)
+write(6,*) ' cellsOnEdge : ', minval(cellsOnEdge), maxval(cellsOnEdge)
+write(6,*) ' nEdgesOnCell : ', minval(nEdgesOnCell), maxval(nEdgesOnCell)
+write(6,*) ' nEdgesOnEdge : ', minval(nEdgesOnEdge), maxval(nEdgesOnEdge)
+write(6,*) ' edgesOnCell : ', minval(edgesOnCell), maxval(edgesOnCell)
+write(6,*) ' edgesOnEdge : ', minval(edgesOnEdge), maxval(edgesOnEdge)
+write(6,*) ' weightsOnEdge : ', minval(weightsOnEdge), maxval(weightsOnEdge)
+write(6,*) ' dvEdge : ', minval(dvEdge), maxval(dvEdge)
+write(6,*) ' dcEdge : ', minval(dcEdge), maxval(dcEdge)
+write(6,*) ' angleEdge : ', minval(angleEdge), maxval(angleEdge)
+write(6,*) ' areaCell : ', minval(areaCell), maxval(areaCell)
+write(6,*) ' areaTriangle : ', minval(areaTriangle), maxval(areaTriangle)
+write(6,*) ' cellsOnCell : ', minval(cellsOnCell), maxval(cellsOnCell)
+write(6,*) ' verticesOnCell : ', minval(verticesOnCell), maxval(verticesOnCell)
+write(6,*) ' verticesOnEdge : ', minval(verticesOnEdge), maxval(verticesOnEdge)
+write(6,*) ' edgesOnVertex : ', minval(edgesOnVertex), maxval(edgesOnVertex)
+write(6,*) ' cellsOnVertex : ', minval(cellsOnVertex), maxval(cellsOnVertex)
+write(6,*) ' kiteAreasOnVertex : ', minval(kiteAreasOnVertex), maxval(kiteAreasOnVertex)
+write(6,*) ' keepCellMask : sum (should be number of 1 values) ', sum(keepCellMask)
+
+
+end subroutine read_grid
+
+
+subroutine write_grid
+implicit none
+
+if (expand_from_unit_sphere) then
+ xCellNew = xCellNew * sphere_radius
+ yCellNew = yCellNew * sphere_radius
+ zCellNew = zCellNew * sphere_radius
+ xEdgeNew = xEdgeNew * sphere_radius
+ yEdgeNew = yEdgeNew * sphere_radius
+ zEdgeNew = zEdgeNew * sphere_radius
+ xVertexNew = xVertexNew * sphere_radius
+ yVertexNew = yVertexNew * sphere_radius
+ zVertexNew = zVertexNew * sphere_radius
+ dcEdgeNew = dcEdgeNew * sphere_radius
+ dvEdgeNew = dvEdgeNew * sphere_radius
+ areaCellNew = areaCellNew * (sphere_radius)**2
+ areaTriangleNew = areaTriangleNew * (sphere_radius)**2
+ kiteAreasOnVertexNew = kiteAreasOnVertexNew * (sphere_radius)**2
+endif
+
+call write_netcdf_init( &
+ nCellsNew, &
+ nEdgesNew, &
+ nVerticesNew, &
+ maxEdgesNew, &
+ nVertLevelsNew, &
+ vertexDegreeNew, &
+ sphere_radius, &
+ on_a_sphere &
+ )
+
+call write_netcdf_fields( &
+ 1, &
+ latCellNew, &
+ lonCellNew, &
+ meshDensityNew, &
+ xCellNew, &
+ yCellNew, &
+ zCellNew, &
+ indexToCellIDNew, &
+ latEdgeNew, &
+ lonEdgeNew, &
+ xEdgeNew, &
+ yEdgeNew, &
+ zEdgeNew, &
+ indexToEdgeIDNew, &
+ latVertexNew, &
+ lonVertexNew, &
+ xVertexNew, &
+ yVertexNew, &
+ zVertexNew, &
+ indexToVertexIDNew, &
+ cellsOnEdgeNew, &
+ nEdgesOnCellNew, &
+ nEdgesOnEdgeNew, &
+ edgesOnCellNew, &
+ edgesOnEdgeNew, &
+ weightsOnEdgeNew, &
+ dvEdgeNew, &
+ dcEdgeNew, &
+ angleEdgeNew, &
+ areaCellNew, &
+ areaTriangleNew, &
+ cellsOnCellNew, &
+ verticesOnCellNew, &
+ verticesOnEdgeNew, &
+ edgesOnVertexNew, &
+ cellsOnVertexNew, &
+ kiteAreasOnVertexNew &
+ )
+
+call write_netcdf_finalize
+
+if (expand_from_unit_sphere) then
+ xCellNew = xCellNew / sphere_radius
+ yCellNew = yCellNew / sphere_radius
+ zCellNew = zCellNew / sphere_radius
+ xEdgeNew = xEdgeNew / sphere_radius
+ yEdgeNew = yEdgeNew / sphere_radius
+ zEdgeNew = zEdgeNew / sphere_radius
+ xVertexNew = xVertexNew / sphere_radius
+ yVertexNew = yVertexNew / sphere_radius
+ zVertexNew = zVertexNew / sphere_radius
+ dcEdgeNew = dcEdgeNew / sphere_radius
+ dvEdgeNew = dvEdgeNew / sphere_radius
+ areaCellNew = areaCellNew / (sphere_radius)**2
+ areaTriangleNew = areaTriangleNew / (sphere_radius)**2
+ kiteAreasOnVertexNew = kiteAreasOnVertexNew / (sphere_radius)**2
+endif
+
+end subroutine write_grid
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Step 4: Check define_kmt routine for bottomDepth and kmt (maxLevelCell) variables
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine define_kmt
+implicit none
+real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
+real (kind=4), allocatable, dimension(:,:) :: ztopo
+integer :: nx, ny, inx, iny, ix, iy, kmt_neighbor_max
+real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax, xmin, xmax
+real :: latmin, latmax, lonmin, lonmax, ridgeDepth, maxdc
+logical :: flag, kmt_flag
+
+!pi = 4.0*atan(1.0)
+!dtr = pi / 180.0
+
+allocate(kmt(nCells))
+kmt = 0
+
+
+do iCell=1,nCells
+ if (keepCellMask(iCell) > 0) then
+ kmt(iCell) = 1 ! any positive value will keep the cell, 0 cells will be culled
+ endif
+enddo
+
+
+!if (cut_domain_from_sphere) then
+! latmin = -30*dtr
+! latmax = +30*dtr
+! lonmin = +10*dtr
+! lonmax = +70*dtr
+! write(6,*) ' lat min ', latmin
+! write(6,*) ' lat max ', latmax
+! where(latCell.lt.latmin) kmt = 0
+! where(latCell.gt.latmax) kmt = 0
+! where(lonCell.lt.lonmin) kmt = 0
+! where(lonCell.gt.lonmax) kmt = 0
+!endif
+
+!if (solid_boundary_in_y) then
+! ymin = minval(yCell)
+! write(6,*) ' minimum yCell ', ymin
+! ymax = maxval(yCell)
+! write(6,*) ' maximum yCell ', ymax
+! where(yCell.lt.1.001*ymin) kmt = 0
+! where(yCell.gt.0.999*ymax) kmt = 0
+!endif
+
+!if (solid_boundary_in_x) then
+! maxdc = maxval(dcEdge)
+! xmin = minval(xCell)
+! write(6,*) ' minimum xCell ', xmin
+! xmax = maxval(xCell)
+! write(6,*) ' maximum xCell ', xmax
+! where(xCell.lt.xmin+maxdc/1.5) kmt = 0
+! where(xCell.gt.xmax-maxdc/1.5) kmt = 0
+!endif
+
+
+allocate(work_kmt(nCells))
+work_kmt = 0.0
+where(kmt.eq.0) work_kmt=1.0
+write(6,*) 'number of cells culled ',sum(work_kmt)
+deallocate(work_kmt)
+
+
+!! Eliminate isolated ocean cells, and make these isolated deep cells
+!! flush with the deepest neighbor.
+!do iCell=1,nCells
+! kmt_neighbor_max = 0
+! do j=1,nEdgesOnCell(iCell)
+! iCell1 = cellsOnCell(j,iCell)
+! kmt_neighbor_max = max(kmt_neighbor_max,kmt(iCell1))
+! enddo
+! if (kmt(iCell).gt.kmt_neighbor_max) then
+! kmt(iCell) = kmt_neighbor_max
+! bottomDepth(iCell) = refBottomDepth(kmt(iCell))
+! endif
+!enddo
+
+!if(eliminate_inland_seas) then
+!call eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+! nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
+! xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
+! KMT)
+!endif
+
+!! do not allow land or PBCs in top layers
+!k = min(top_layers_without_land,nVertLevelsMOD)
+!where(kmt.gt.0.and.kmt.le.k)
+! bottomDepth = refBottomDepth(k)
+! kmt=k
+!endwhere
+
+end subroutine define_kmt
+
+
+
+subroutine define_mapping
+implicit none
+
+allocate(cellMap(nCells))
+allocate(edgeMap(nEdges))
+allocate(vertexMap(nVertices))
+cellMap = 0; edgeMap = 0; vertexMap = 0
+
+j=1
+do i=1,nCells
+if(kmt(i).ne.0) then
+ cellMap(i) = j
+ j=j+1
+endif
+write(10,*) i, cellMap(i)
+enddo
+
+j=1
+do i=1,nEdges
+iCell1 = cellsOnEdge(1,i)
+iCell2 = cellsOnEdge(2,i)
+if(kmt(iCell1).ne.0.or.kmt(iCell2).ne.0) then
+ edgeMap(i)=j
+ j=j+1
+endif
+write(11,*) i,edgeMap(i)
+enddo
+
+j=1
+do i=1,nVertices
+iCell1 = cellsOnVertex(1,i)
+iCell2 = cellsOnVertex(2,i)
+iCell3 = cellsOnVertex(3,i)
+if(kmt(iCell1).ne.0.or.kmt(iCell2).ne.0.or.kmt(iCell3).ne.0) then
+ vertexMap(i)=j
+ j=j+1
+endif
+write(12,*) i,vertexMap(i)
+enddo
+
+nCellsNew = 0
+do i=1,nCells
+if(cellMap(i).ne.0) nCellsNew = nCellsNew + 1
+enddo
+
+nEdgesNew = 0
+do i=1,nEdges
+if(edgeMap(i).ne.0) nEdgesNew = nEdgesNew + 1
+enddo
+
+nVerticesNew = 0
+do i=1,nVertices
+if(vertexMap(i).ne.0) nVerticesNew = nVerticesNew + 1
+enddo
+
+write(6,*) ' mesh mapping found '
+write(6,*) nCells, nCellsNew
+write(6,*) nEdges, nEdgesNew
+write(6,*) nVertices, nVerticesNew
+
+allocate(indexToCellIDNew(nCellsNew))
+allocate(indexToEdgeIDNew(nEdgesNew))
+allocate(indexToVertexIDNew(nVerticesNew))
+indextoCellIDNew = 0; indexToEdgeIDNew = 0; indexToVertexIDNew = 0
+
+do i=1,nCellsNew
+indexToCellIDNew(i)=i
+enddo
+
+do i=1,nEdgesNew
+indexToEdgeIDNew(i)=i
+enddo
+
+do i=1,nVerticesNew
+indexToVertexIDNew(i)=i
+enddo
+
+end subroutine define_mapping
+
+
+subroutine map_vectors
+implicit none
+
+allocate(xCellNew(nCellsNew))
+allocate(yCellNew(nCellsNew))
+allocate(zCellNew(nCellsNew))
+allocate(normalsNew(3,nEdgesNew))
+allocate(latCellNew(nCellsNew))
+allocate(lonCellNew(nCellsNew))
+allocate(meshDensityNew(nCellsNew))
+allocate(meshSpacingNew(nCellsNew))
+allocate(xEdgeNew(nEdgesNew))
+allocate(yEdgeNew(nEdgesNew))
+allocate(zEdgeNew(nEdgesNew))
+allocate(latEdgeNew(nEdgesNew))
+allocate(lonEdgeNew(nEdgesNew))
+allocate(xVertexNew(nVerticesNew))
+allocate(yVertexNew(nVerticesNew))
+allocate(zVertexNew(nVerticesNew))
+allocate(latVertexNew(nVerticesNew))
+allocate(lonVertexNew(nVerticesNew))
+allocate(dcEdgeNew(nEdgesNew))
+allocate(dvEdgeNew(nEdgesNew))
+allocate(angleEdgeNew(nEdgesNew))
+allocate(areaCellNew(nCellsNew))
+allocate(areaTriangleNew(nVerticesNew))
+
+
+
+
+xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0; meshDensityNew=1.0; meshSpacingNew=0.0
+xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
+xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
+
+
+do i=1,nCells
+jNew = cellMap(i)
+if(jNew.ne.0) then
+ xCellNew(jNew)=xCell(i)
+ yCellNew(jNew)=yCell(i)
+ zCellNew(jNew)=zCell(i)
+ latCellNew(jNew)=latCell(i)
+ lonCellNew(jNew)=lonCell(i)
+ meshDensityNew(jNew)=meshDensity(i)
+ areaCellNew(jNew)=areaCell(i)
+endif
+enddo
+
+do i=1,nEdges
+jNew = edgeMap(i)
+if(jNew.ne.0) then
+ xEdgeNew(jNew)=xEdge(i)
+ yEdgeNew(jNew)=yEdge(i)
+ zEdgeNew(jNew)=zEdge(i)
+ latEdgeNew(jNew)=latEdge(i)
+ lonEdgeNew(jNew)=lonEdge(i)
+ dcEdgeNew(jNew) = dcEdge(i)
+ dvEdgeNew(jNew) = dvEdge(i)
+ angleEdgeNew(jNew) = angleEdge(i)
+endif
+enddo
+
+do i=1,nVertices
+jNew = vertexMap(i)
+if(jNew.ne.0) then
+ xVertexNew(jNew)=xVertex(i)
+ yVertexNew(jNew)=yVertex(i)
+ zVertexNew(jNew)=zVertex(i)
+ latVertexNew(jNew)=latVertex(i)
+ lonVertexNew(jNew)=lonVertex(i)
+ areaTriangleNew(jNew)=areaTriangle(i)
+endif
+enddo
+
+deallocate(xCell)
+deallocate(yCell)
+deallocate(zCell)
+deallocate(latCell)
+deallocate(lonCell)
+deallocate(meshDensity)
+deallocate(xEdge)
+deallocate(yEdge)
+deallocate(zEdge)
+deallocate(latEdge)
+deallocate(lonEdge)
+deallocate(xVertex)
+deallocate(yVertex)
+deallocate(zVertex)
+deallocate(latVertex)
+deallocate(lonVertex)
+deallocate(dcEdge)
+deallocate(dvEdge)
+
+end subroutine map_vectors
+
+
+
+subroutine map_connectivity
+implicit none
+
+allocate(cellsOnEdgeNew(TWONew,nEdgesNew))
+allocate(boundaryEdgeNew(nVertLevelsNew,nEdgesNew))
+allocate(flipVerticesOnEdgeOrdering(nEdgesNew))
+cellsOnEdgeNew(:,:) = 0
+boundaryEdgeNew(:,:) = 0
+flipVerticesOnEdgeOrdering(:) = 0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+iCell1 = cellsOnEdge(1,iEdge)
+iCell2 = cellsOnEdge(2,iEdge)
+iCell1New = cellMap(iCell1)
+iCell2New = cellMap(iCell2)
+cellsOnEdgeNew(1,iEdgeNew) = iCell1New
+cellsOnEdgeNew(2,iEdgeNew) = iCell2New
+if(iCell1New.eq.0.or.iCell2New.eq.0) boundaryEdgeNew(:,iEdgeNew) = 1
+if(iCell1New.eq.0.and.iCell2New.eq.0) stop "cellsOnEdge"
+if(iCell1New.eq.0) then
+ cellsOnEdgeNew(1,iEdgeNew) = iCell2New
+ cellsOnEdgeNew(2,iEdgeNew) = iCell1New
+ flipVerticesOnEdgeOrdering(iEdgeNew) = 1
+endif
+enddo
+deallocate(cellsOnEdge)
+
+allocate(verticesOnEdgeNew(TWONew,nEdgesNew))
+allocate(boundaryVertexNew(nVertLevelsNew,nVerticesNew))
+verticesOnEdgeNew(:,:) = 0
+boundaryVertexNew(:,:) = 0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+iVertex1 = VerticesOnEdge(1,iEdge)
+iVertex2 = VerticesOnEdge(2,iEdge)
+iVertex1New = vertexMap(iVertex1)
+iVertex2New = vertexMap(iVertex2)
+if(iVertex1New.eq.0.or.iVertex2New.eq.0) stop "verticesOnEdge"
+if(flipVerticesOnEdgeOrdering(iEdgeNew).eq.0) then
+ verticesOnEdgeNew(1,iEdgeNew) = iVertex1New
+ verticesOnEdgeNew(2,iEdgeNew) = iVertex2New
+else
+ verticesOnEdgeNew(1,iEdgeNew) = iVertex2New
+ verticesOnEdgeNew(2,iEdgeNew) = iVertex1New
+endif
+if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
+ boundaryVertexNew(:,iVertex1New)=1
+ boundaryVertexNew(:,iVertex2New)=1
+endif
+enddo
+deallocate(verticesOnEdge)
+
+allocate(nEdgesOnEdgeNew(nEdgesNew))
+allocate(edgesOnEdgeNew(maxEdges2,nEdgesNew))
+allocate(weightsOnEdgeNew(maxEdges2,nEdgesNew))
+nEdgesOnEdgeNew(:) = 0
+edgesOnEdgeNew(:,:) = 0
+weightsOnEdgeNew(:,:) = 0.0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
+ nEdgesOnEdgeNew(iEdgeNew) = 0
+ edgesOnEdgeNew(:,iEdgeNew) = 0
+ weightsOnEdgeNew(:,iEdgeNew) = 0.0
+else
+ nEdgesOnEdgeNew(iEdgeNew) = nEdgesOnEdge(iEdge)
+ do i=1,nEdgesOnEdgeNew(iEdgeNew)
+ jEdge = edgesOnEdge(i,iEdge)
+ jEdgeNew = edgeMap(jEdge)
+ if(jEdgeNew.eq.0) stop "jEdgeNew"
+ edgesOnEdgeNew(i,iEdgeNew)=jEdgeNew
+ weightsOnEdgeNew(i,iEdgeNew) = weightsOnEdge(i,iEdge)
+ enddo
+endif
+enddo
+deallocate(nEdgesOnEdge)
+deallocate(edgesOnEdge)
+deallocate(weightsOnEdge)
+
+allocate(cellsOnCellNew(maxEdges,nCellsNew))
+allocate(nEdgesOnCellNew(nCellsNew))
+cellsOnCellNew = 0
+nEdgesOnCellNew = 0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+nEdgesOnCellNew(iCellNew)=nEdgesOnCell(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = cellsOnCell(i,iCell)
+jNew = cellMap(j)
+cellsOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
+deallocate(cellsOnCell)
+deallocate(nEdgesOnCell)
+
+allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
+edgesOnCellNew(:,:) = 0
+meshSpacingNew(:) = 0.0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = edgesOnCell(i,iCell)
+jNew = edgeMap(j)
+if(jNew.eq.0) stop "edgesOnCell"
+edgesOnCellNew(i,iCellNew) = jNew
+meshSpacingNew(iCellNew) = meshSpacingNew(iCellNew) + dcEdgeNew(jNew)/nEdgesOnCellNew(iCellNew)
+enddo
+enddo
+deallocate(edgesOnCell)
+
+allocate(verticesOnCellNew(maxEdgesNew,nCellsNew))
+verticesOnCellNew(:,:)=0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j=verticesOnCell(i,iCell)
+jNew = vertexMap(j)
+if(jNew.eq.0) stop "verticesOnCell"
+verticesOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
+deallocate(verticesOnCell)
+
+allocate(cellsOnVertexNew(vertexDegreeNew,nVerticesNew))
+allocate(kiteAreasOnVertexNew(vertexDegreeNew,nVerticesNew))
+cellsOnVertexNew = 0
+kiteAreasOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=cellsOnVertex(i,iVertex)
+jNew=cellMap(j)
+if(jNew.eq.0) then
+ kiteAreasOnVertexNew(i,iVertexNew)=0
+else
+ kiteAreasOnVertexNew(i,iVertexNew)=kiteAreasOnVertex(i,iVertex)
+endif
+cellsOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
+deallocate(cellsOnVertex)
+deallocate(kiteAreasOnVertex)
+
+areaTriangleNew = 0
+do iVertex=1,nVerticesNew
+do i=1,vertexDegree
+areaTriangleNew(iVertex) = areaTriangleNew(iVertex) + kiteAreasOnVertexNew(i,iVertex)
+enddo
+enddo
+
+allocate(edgesOnVertexNew(vertexDegreeNew, nVerticesNew))
+edgesOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=edgesOnVertex(i,iVertex)
+jNew=edgeMap(j)
+edgesOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
+deallocate(edgesOnVertex)
+
+! find normals
+normalsNew = 0.0
+do iEdge=1,nEdgesNew
+cell1 = cellsOnEdgeNew(1,iEdge)
+cell2 = cellsOnEdgeNew(2,iEdge)
+if(cell1.eq.0.or.cell2.eq.0) cycle
+c1(1) = xCellNew(cell1); c1(2) = yCellNew(cell1); c1(3) = zCellNew(cell1)
+c2(1) = xCellNew(cell2); c2(2) = yCellNew(cell2); c2(3) = zCellNew(cell2)
+distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+
+if(on_a_sphere.eq.'YES') then
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+else
+! if(distance.gt.0.5*Lx) then
+! write(6,*) ' periodic edge ', iEdge, distance
+! write(6,10) ' c1 ', c1(:)
+! write(6,10) ' c2 ', c2(:)
+! r = c2(1) - c1(1)
+! if(r.gt.0) c2(1) = c2(1) - Lx
+! if(r.lt.0) c2(1) = c2(1) + Lx
+! distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+! write(6,*) ' periodic edge fix ', iEdge, r, distance
+! endif
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+endif
+enddo
+10 format(a20,3e15.5)
+
+end subroutine map_connectivity
+
+
+end program crop_landice_grid
Added: branches/land_ice_projects/grid_tools/crop_landice_grid/define_keepCellMask.py
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/define_keepCellMask.py (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/define_keepCellMask.py 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,91 @@
+#!/usr/bin/python
+# Script for adding a field named cullMask to an MPAS land ice grid for use with the crop_landice_grid.F tool that actually culls the unwanted cells.
+# Matt Hoffman, February 28, 2013
+
+import sys
+import numpy as np
+from optparse import OptionParser
+import matplotlib.pyplot as plt
+try:
+ from Scientific.IO.NetCDF import NetCDFFile
+ netCDF_module = 'Scientific.IO.NetCDF'
+except ImportError:
+ try:
+ from netCDF4 import Dataset as NetCDFFile
+ netCDF_module = 'netCDF4'
+ except ImportError:
+ print 'Unable to import any of the following python modules:'
+ print ' Scientific.IO.NetCDF </font>
<font color="gray"> netcdf4 '
+ print 'One of them must be installed.'
+ raise ImportError('No netCDF module found')
+
+
+print "** Gathering information."
+parser = OptionParser()
+parser.add_option("-f", "--file", dest="file", help="grid file to modify; default: land_ice_grid.nc", metavar="FILE")
+
+
+options, args = parser.parse_args()
+
+if not options.file:
+ print "No grid filename provided. Using land_ice_grid.nc."
+ options.file = "land_ice_grid.nc"
+
+try:
+ f = NetCDFFile(options.file,'r+')
+except:
+ sys.exit('Error loading grid file.')
+
+
+# get the right representation of the datatype
+if netCDF_module == 'Scientific.IO.NetCDF':
+ datatype = 'i'
+else:
+ datatype = f.variables['indexToCellID'].dtype # Get the datatype for double precision float
+
+# Try to add the new variable
+if 'keepCellMask' not in f.variables:
+ f.createVariable('keepCellMask', datatype, ('nCells',))
+
+keepCellMask = f.variables['keepCellMask'][:]
+thickness = f.variables['thickness'][:]
+
+keepCellMask[:] = 0
+
+# ===== Various methods for defining the mask ====
+maskmeth = 2
+
+# 1. only keep cells with ice
+if maskmeth == 1:
+ keepCellMask[thickness[0,:] > 0.0] = 1
+
+# 2. add a buffer of X cells around the ice
+elif maskmeth == 2:
+
+ buffersize=4 # number of cells to expand
+
+ keepCellMask[thickness[0,:] > 0.0] = 1 # start with the mask being where ice is
+ print 'Num of cells to keep:', sum(keepCellMask)
+ cellsOnCell = f.variables['cellsOnCell']
+ for i in range(buffersize):
+ print 'Starting buffer loop ', i+1
+ keepCellMaskNew = keepCellMask[:]+0 # make a copy to edit --- the +0 makes an actual copy rather than pointing to the memory
+ for iCell in range(len(keepCellMask)):
+ for neighbor in cellsOnCell[iCell,:]-1: # the -1 converts from the fortran indexing in the variable to python indexing
+ if keepCellMask[neighbor] == 1: # if any neighbors are already being kept on the old mask then keep this cell too. This will get all (or almost all) of the ice cells, but they have already been assigned so they don't matter. What we care about is getting cells that are not ice that have on or more ice neighbors.
+ keepCellMaskNew[iCell] = 1
+ keepCellMask = keepCellMaskNew[:]+0 # after we've looped over all cells assign the new mask to the variable we need (either for another loop around the domain or to write out)
+ print 'Num of cells to keep:', sum(keepCellMask)
+fig = plt.figure(1, facecolor='w')
+ax = fig.add_subplot(111, aspect='equal')
+plt.scatter(f.variables['xCell'][:], f.variables['yCell'][:], 50, keepCellMask, edgecolors='none') #, vmin=minval, vmax=maxval)
+plt.colorbar()
+plt.draw()
+plt.show()
+
+f.variables['keepCellMask'][:] = keepCellMask
+
+#f.close()
+
+
+
Added: branches/land_ice_projects/grid_tools/crop_landice_grid/module_cullLoops.F
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/module_cullLoops.F (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/module_cullLoops.F 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,84 @@
+module cullLoops
+
+ public :: eliminateLoops
+
+ contains
+
+ subroutine eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+ nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
+ xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
+ KMT)
+
+ implicit none
+
+ ! intent (in)
+ integer :: nCells, nEdges, nVertices, maxEdges, vertexDegree
+ integer :: nEdgesOnCell(nCells), cellsOnCell(maxEdges,nCells), verticesOnEdge(2,nEdges)
+ integer :: cellsOnVertex(vertexDegree,nVertices), edgesOnCell(maxEdges,nCells)
+ real :: lonCell(nCells), latCell(nCells)
+ real :: xCell(nCells), yCell(nCells), zCell(nCells)
+ real :: xEdge(nEdges), yEdge(nEdges), zEdge(nEdges)
+ real :: xVertex(nVertices), yVertex(nVertices), zVertex(nVertices)
+ integer :: edgeList(nEdges), iCellMask(nCells)
+
+ ! intent(inout)
+ integer, intent(inout) :: KMT(ncells)
+
+ ! local workspace
+ integer :: iCell, jCell, oCell, lCell, iEdge, i, kCell, iSharedEdge, iStartEdge, iSave, iSweep
+ integer :: iEdgeCounter, nEdgesInLoop(nCells), iCellAhead, LeftTurns, RightTurns
+ logical :: connected, atBoundary, moveSouth, moveEast, atGrenwich
+ real :: lat, rlat, rlon, rCenter(3), s(3), t(3), q(3), rCross, mylon, mylat, pi
+
+ integer, dimension(:), pointer :: cellStack
+ integer, dimension(:), pointer :: oceanMask
+ integer :: iCellStart, nStack, addedCells
+ real :: latStart, lonStart
+
+ write(6,*) 'Culling inland seas.....'
+
+ allocate(cellStack(nCells/2))
+ allocate(oceanMask(nCells))
+
+ oceanMask = 0
+ addedCells = 0
+
+ iCellStart = maxloc(kmt, dim=1)
+
+ write(6,*) 'Starting index. ', iCellStart
+ write(6,*) 'lat, lon: ', latCell(iCellStart), lonCell(iCellStart)
+ write(6,*) 'Starting kmt: ', kmt(iCellStart)
+
+ nStack = 1
+ cellStack(nStack) = iCellStart
+ oceanMask(iCellStart) = 1
+ addedCells = 1
+
+ do while(nStack > 0)
+ oCell = cellStack(nStack)
+ nStack = nStack - 1
+ !write(6,*) ' Working on cell ', oCell, addedCells, nStack
+
+ do i = 1, nEdgesOnCell(oCell)
+ iCell = cellsOnCell(i, oCell)
+
+ if(kmt(iCell) > 0 .and. oceanMask(iCell) == 0) then
+ nStack = nStack + 1
+ cellStack(nStack) = iCell
+ oceanMask(iCell) = 1
+ addedCells = addedCells + 1
+ end if
+ end do
+ end do
+
+ where(oceanMask == 0) kmt(:) = 0
+
+ write(6,*) addedCells, ' total cells have been in the stack.'
+ write(6,*) 'Done culling inland seas.....'
+
+ deallocate(cellStack)
+ deallocate(oceanMask)
+
+ end subroutine eliminateLoops
+
+end module cullLoops
Added: branches/land_ice_projects/grid_tools/crop_landice_grid/module_read_netcdf.F
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/module_read_netcdf.F (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/module_read_netcdf.F 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,476 @@
+module read_netcdf
+
+ integer :: rd_ncid
+ integer :: rdDimIDTime
+ integer :: rdDimIDnCells
+ integer :: rdDimIDnEdges
+ integer :: rdDimIDnVertices
+ integer :: rdDimIDmaxEdges
+ integer :: rdDimIDmaxEdges2
+ integer :: rdDimIDnVertLevels
+ integer :: rdDimIDTWO
+ integer :: rdDimIDvertexDegree
+ integer :: rdVarIDlatCell
+ integer :: rdVarIDlonCell
+ integer :: rdVarIDmeshDensity
+ integer :: rdVarIDxCell
+ integer :: rdVarIDyCell
+ integer :: rdVarIDzCell
+ integer :: rdVarIDindexToCellID
+ integer :: rdVarIDlatEdge
+ integer :: rdVarIDlonEdge
+ integer :: rdVarIDxEdge
+ integer :: rdVarIDyEdge
+ integer :: rdVarIDzEdge
+ integer :: rdVarIDindexToEdgeID
+ integer :: rdVarIDlatVertex
+ integer :: rdVarIDlonVertex
+ integer :: rdVarIDxVertex
+ integer :: rdVarIDyVertex
+ integer :: rdVarIDzVertex
+ integer :: rdVarIDindexToVertexID
+ integer :: rdVarIDcellsOnEdge
+ integer :: rdVarIDnEdgesOnCell
+ integer :: rdVarIDnEdgesOnEdge
+ integer :: rdVarIDedgesOnCell
+ integer :: rdVarIDedgesOnEdge
+ integer :: rdVarIDweightsOnEdge
+ integer :: rdVarIDdvEdge
+ integer :: rdVarIDdcEdge
+ integer :: rdVarIDangleEdge
+ integer :: rdVarIDareaCell
+ integer :: rdVarIDareaTriangle
+ integer :: rdVarIDcellsOnCell
+ integer :: rdVarIDverticesOnCell
+ integer :: rdVarIDverticesOnEdge
+ integer :: rdVarIDedgesOnVertex
+ integer :: rdVarIDcellsOnVertex
+ integer :: rdVarIDkiteAreasOnVertex
+ integer :: rdVarIDkeepCellMask
+
+
+ integer :: rdLocalnCells
+ integer :: rdLocalnEdges
+ integer :: rdLocalnVertices
+ integer :: rdLocalmaxEdges
+ integer :: rdLocalmaxEdges2
+ integer :: rdLocalnVertLevels
+ integer :: rdLocalTWO
+ integer :: rdLocalvertexDegree
+
+ contains
+
+ subroutine read_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ maxEdges2, &
+ nVertLevels, &
+ TWO, &
+ vertexDegree &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(out) :: nCells
+ integer, intent(out) :: nEdges
+ integer, intent(out) :: nVertices
+ integer, intent(out) :: maxEdges
+ integer, intent(out) :: maxEdges2
+ integer, intent(out) :: nVertLevels
+ integer, intent(out) :: TWO
+ integer, intent(out) :: vertexDegree
+
+ integer :: nferr
+
+
+ nferr = nf_open('land_ice_grid.nc', NF_SHARE, rd_ncid)
+
+ !
+ ! Get IDs for variable dimensions
+ !
+ nferr = nf_inq_unlimdim(rd_ncid, rdDimIDTime)
+ nferr = nf_inq_dimid(rd_ncid, 'nCells', rdDimIDnCells)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnCells, rdLocalnCells)
+ nferr = nf_inq_dimid(rd_ncid, 'nEdges', rdDimIDnEdges)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnEdges, rdLocalnEdges)
+ nferr = nf_inq_dimid(rd_ncid, 'nVertices', rdDimIDnVertices)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertices, rdLocalnVertices)
+ nferr = nf_inq_dimid(rd_ncid, 'maxEdges', rdDimIDmaxEdges)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges, rdLocalmaxEdges)
+ nferr = nf_inq_dimid(rd_ncid, 'maxEdges2', rdDimIDmaxEdges2)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges2, rdLocalmaxEdges2)
+ nferr = nf_inq_dimid(rd_ncid, 'nVertLevels', rdDimIDnVertLevels)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertLevels, rdLocalnVertLevels)
+ nferr = nf_inq_dimid(rd_ncid, 'vertexDegree', rdDimIDvertexDegree)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDvertexDegree, rdLocalvertexDegree)
+ nferr = nf_inq_dimid(rd_ncid, 'TWO', rdDimIDTWO)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDTWO, rdLocalTWO)
+
+
+ nCells = rdLocalnCells
+ nEdges = rdLocalnEdges
+ nVertices = rdLocalnVertices
+ maxEdges = rdLocalmaxEdges
+ maxEdges2 = rdLocalmaxEdges2
+ nVertLevels = rdLocalnVertLevels
+ vertexDegree = rdLocalvertexDegree
+ TWO = rdLocalTWO
+
+ !
+ ! Get IDs for variables
+ !
+ nferr = nf_inq_varid(rd_ncid, 'latCell', rdVarIDlatCell)
+ nferr = nf_inq_varid(rd_ncid, 'lonCell', rdVarIDlonCell)
+ nferr = nf_inq_varid(rd_ncid, 'meshDensity', rdVarIDmeshDensity)
+ nferr = nf_inq_varid(rd_ncid, 'xCell', rdVarIDxCell)
+ nferr = nf_inq_varid(rd_ncid, 'yCell', rdVarIDyCell)
+ nferr = nf_inq_varid(rd_ncid, 'zCell', rdVarIDzCell)
+ nferr = nf_inq_varid(rd_ncid, 'indexToCellID', rdVarIDindexToCellID)
+ nferr = nf_inq_varid(rd_ncid, 'latEdge', rdVarIDlatEdge)
+ nferr = nf_inq_varid(rd_ncid, 'lonEdge', rdVarIDlonEdge)
+ nferr = nf_inq_varid(rd_ncid, 'xEdge', rdVarIDxEdge)
+ nferr = nf_inq_varid(rd_ncid, 'yEdge', rdVarIDyEdge)
+ nferr = nf_inq_varid(rd_ncid, 'zEdge', rdVarIDzEdge)
+ nferr = nf_inq_varid(rd_ncid, 'indexToEdgeID', rdVarIDindexToEdgeID)
+ nferr = nf_inq_varid(rd_ncid, 'latVertex', rdVarIDlatVertex)
+ nferr = nf_inq_varid(rd_ncid, 'lonVertex', rdVarIDlonVertex)
+ nferr = nf_inq_varid(rd_ncid, 'xVertex', rdVarIDxVertex)
+ nferr = nf_inq_varid(rd_ncid, 'yVertex', rdVarIDyVertex)
+ nferr = nf_inq_varid(rd_ncid, 'zVertex', rdVarIDzVertex)
+ nferr = nf_inq_varid(rd_ncid, 'indexToVertexID', rdVarIDindexToVertexID)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnEdge', rdVarIDcellsOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'nEdgesOnCell', rdVarIDnEdgesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'nEdgesOnEdge', rdVarIDnEdgesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnCell', rdVarIDedgesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnEdge', rdVarIDedgesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'weightsOnEdge', rdVarIDweightsOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'dvEdge', rdVarIDdvEdge)
+ nferr = nf_inq_varid(rd_ncid, 'dcEdge', rdVarIDdcEdge)
+ nferr = nf_inq_varid(rd_ncid, 'angleEdge', rdVarIDangleEdge)
+ nferr = nf_inq_varid(rd_ncid, 'areaCell', rdVarIDareaCell)
+ nferr = nf_inq_varid(rd_ncid, 'areaTriangle', rdVarIDareaTriangle)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnCell', rdVarIDcellsOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'verticesOnCell', rdVarIDverticesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'verticesOnEdge', rdVarIDverticesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnVertex', rdVarIDedgesOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnVertex', rdVarIDcellsOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'kiteAreasOnVertex', rdVarIDkiteAreasOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'keepCellMask', rdVarIDkeepCellMask)
+ end subroutine read_netcdf_init
+
+
+ subroutine read_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ meshDensity, &
+ 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, &
+ keepCellMask &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: time
+ real (kind=8), dimension(:), intent(out) :: latCell
+ real (kind=8), dimension(:), intent(out) :: lonCell
+ real (kind=8), dimension(:), intent(out) :: meshDensity
+ real (kind=8), dimension(:), intent(out) :: xCell
+ real (kind=8), dimension(:), intent(out) :: yCell
+ real (kind=8), dimension(:), intent(out) :: zCell
+ integer, dimension(:), intent(out) :: indexToCellID
+ real (kind=8), dimension(:), intent(out) :: latEdge
+ real (kind=8), dimension(:), intent(out) :: lonEdge
+ real (kind=8), dimension(:), intent(out) :: xEdge
+ real (kind=8), dimension(:), intent(out) :: yEdge
+ real (kind=8), dimension(:), intent(out) :: zEdge
+ integer, dimension(:), intent(out) :: indexToEdgeID
+ real (kind=8), dimension(:), intent(out) :: latVertex
+ real (kind=8), dimension(:), intent(out) :: lonVertex
+ real (kind=8), dimension(:), intent(out) :: xVertex
+ real (kind=8), dimension(:), intent(out) :: yVertex
+ real (kind=8), dimension(:), intent(out) :: zVertex
+ integer, dimension(:), intent(out) :: indexToVertexID
+ integer, dimension(:,:), intent(out) :: cellsOnEdge
+ integer, dimension(:), intent(out) :: nEdgesOnCell
+ integer, dimension(:), intent(out) :: nEdgesOnEdge
+ integer, dimension(:,:), intent(out) :: edgesOnCell
+ integer, dimension(:,:), intent(out) :: edgesOnEdge
+ real (kind=8), dimension(:,:), intent(out) :: weightsOnEdge
+ real (kind=8), dimension(:), intent(out) :: dvEdge
+ real (kind=8), dimension(:), intent(out) :: dcEdge
+ real (kind=8), dimension(:), intent(out) :: angleEdge
+ real (kind=8), dimension(:), intent(out) :: areaCell
+ real (kind=8), dimension(:), intent(out) :: areaTriangle
+ integer, dimension(:,:), intent(out) :: cellsOnCell
+ integer, dimension(:,:), intent(out) :: verticesOnCell
+ integer, dimension(:,:), intent(out) :: verticesOnEdge
+ integer, dimension(:,:), intent(out) :: edgesOnVertex
+ integer, dimension(:,:), intent(out) :: cellsOnVertex
+ real (kind=8), dimension(:,:), intent(out) :: kiteAreasOnVertex
+ integer, dimension(:), intent(out) :: keepCellMask
+
+ logical :: meshDensityPresent
+
+ integer :: nferr
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ meshDensityPresent = .false.
+
+ start1(1) = 1
+
+ start2(1) = 1
+ start2(2) = 1
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatCell, start1, count1, latCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonCell, start1, count1, lonCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_inq_varid(rd_ncid, 'meshDensity', rdVarIDmeshDensity)
+ if(nferr.eq.0) then
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDmeshDensity, start1, count1, meshDensity)
+ else
+ meshDensity=1.0
+ write(6,*) ' mesh density not present ', nferr, rdVarIDmeshDensity
+ endif
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxCell, start1, count1, xCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyCell, start1, count1, yCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzCell, start1, count1, zCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToCellID, start1, count1, indexToCellID)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatEdge, start1, count1, latEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonEdge, start1, count1, lonEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxEdge, start1, count1, xEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyEdge, start1, count1, yEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzEdge, start1, count1, zEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToEdgeID, start1, count1, indexToEdgeID)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatVertex, start1, count1, latVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonVertex, start1, count1, lonVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxVertex, start1, count1, xVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyVertex, start1, count1, yVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzVertex, start1, count1, zVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToVertexID, start1, count1, indexToVertexID)
+
+ start2(2) = 1
+ count2( 1) = rdLocalTWO
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnEdge, start2, count2, cellsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnCell, start2, count2, edgesOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges2
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnEdge, start2, count2, edgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges2
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDweightsOnEdge, start2, count2, weightsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDdvEdge, start1, count1, dvEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDdcEdge, start1, count1, dcEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDangleEdge, start1, count1, angleEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDareaCell, start1, count1, areaCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDareaTriangle, start1, count1, areaTriangle)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnCell, start2, count2, cellsOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnCell, start2, count2, verticesOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalTWO
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnEdge, start2, count2, verticesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDkeepCellMask, start1, count1, keepCellMask)
+
+ end subroutine read_netcdf_fields
+
+
+ subroutine read_netcdf_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(rd_ncid)
+
+ end subroutine read_netcdf_finalize
+
+end module read_netcdf
Added: branches/land_ice_projects/grid_tools/crop_landice_grid/module_write_netcdf.F
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/module_write_netcdf.F (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/module_write_netcdf.F 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,507 @@
+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 :: wrDimIDnMonths
+ integer :: wrVarIDlatCell
+ integer :: wrVarIDlonCell
+ integer :: wrVarIDmeshDensity
+ 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 :: wrVarIDmaxLevelCell
+ 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 :: wrVarIDbottomDepth
+ integer :: wrVarIDu
+ integer :: wrVarIDboundaryEdge
+ integer :: wrVarIDboundaryVertex
+ integer :: wrVarIDu_src
+ integer :: wrVarIDwindStressMonthly
+ integer :: wrVarIDh
+ integer :: wrVarIDrho
+ integer :: wrVarIDtemperature
+ integer :: wrVarIDsalinity
+ integer :: wrVarIDtracer1
+ integer :: wrVarIDtemperatureRestore
+ integer :: wrVarIDsalinityRestore
+ integer :: wrVarIDtemperatureRestoreMonthly
+ integer :: wrVarIDsalinityRestoreMonthly
+ integer :: wrVarIDhZLevel
+ integer :: wrVarIDrefBottomDepth
+
+ integer :: wrLocalnCells
+ integer :: wrLocalnEdges
+ integer :: wrLocalnVertices
+ integer :: wrLocalmaxEdges
+ integer :: wrLocalnVertLevels
+ integer :: wrLocalvertexDegree
+
+
+ contains
+
+ subroutine write_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ nVertLevels, &
+ vertexDegree, &
+ sphere_radius, &
+ on_a_sphere &
+ )
+
+ 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) :: vertexDegree
+ character (len=16) :: on_a_sphere
+ real*8 :: sphere_radius
+
+
+ integer :: nferr
+ integer, dimension(10) :: dimlist
+
+
+ wrLocalnCells = nCells
+ wrLocalnEdges = nEdges
+ wrLocalnVertices = nVertices
+ wrLocalmaxEdges = maxEdges
+ wrLocalnVertLevels = nVertLevels
+ wrLocalvertexDegree = vertexDegree
+
+ nferr = nf_create('land_ice_grid_culled.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, '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, 'meshDensity', NF_DOUBLE, 1, dimlist, wrVarIDmeshDensity)
+ 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) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'maxLevelCell', NF_INT, 1, dimlist, wrVarIDmaxLevelCell)
+ 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)
+
+
+
+
+ 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)
+
+ nferr = nf_enddef(wr_ncid)
+
+ end subroutine write_netcdf_init
+
+
+ subroutine write_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ meshDensity, &
+ 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 &
+ )
+
+ 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) :: meshDensity
+ 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
+
+
+
+ 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, wrVarIDmeshDensity, start1, count1, meshDensity)
+
+ 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) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+
+
+
+ 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_projects/grid_tools/crop_landice_grid/utilities.F
===================================================================
--- branches/land_ice_projects/grid_tools/crop_landice_grid/utilities.F (rev 0)
+++ branches/land_ice_projects/grid_tools/crop_landice_grid/utilities.F 2013-02-28 19:12:01 UTC (rev 2521)
@@ -0,0 +1,781 @@
+module utilities
+
+contains
+
+subroutine write_OpenDX( on_a_sphere, &
+ nCells, &
+ nVertices, &
+ nEdges, &
+ vertexDegree, &
+ maxEdges, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ cellsOnVertex, &
+ edgesOnCell, &
+ areaCell, &
+ maxLevelCell, &
+ meshSpacing, &
+ depthCell, &
+ SST, &
+ kiteAreasOnVertex )
+
+ implicit none
+
+ character (len=16), intent(in) :: on_a_sphere
+ integer, intent(in) :: nCells, nVertices, vertexDegree, nEdges, maxEdges
+ real (kind=8), dimension(nCells), intent(inout) :: xCell
+ real (kind=8), dimension(nCells), intent(inout) :: yCell
+ real (kind=8), dimension(nCells), intent(inout) :: zCell
+ real (kind=8), dimension(nVertices), intent(inout) :: xVertex
+ real (kind=8), dimension(nVertices), intent(inout) :: yVertex
+ real (kind=8), dimension(nVertices), intent(inout) :: zVertex
+ real (kind=8), dimension(nEdges), intent(inout) :: xEdge
+ real (kind=8), dimension(nEdges), intent(inout) :: yEdge
+ real (kind=8), dimension(nEdges), intent(inout) :: zEdge
+ integer, dimension(nCells), intent(in) :: nEdgesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: edgesOnCell
+ integer, dimension(2,nEdges), intent(in) :: verticesOnEdge
+ integer, dimension(vertexDegree, nVertices), intent(in) :: cellsOnVertex
+ integer, dimension(nCells), intent(in) :: maxLevelCell
+ real (kind=8), dimension(nCells), intent(in) :: areaCell
+ real (kind=8), dimension(nCells), intent(in) :: depthCell, SST, meshSpacing
+ real (kind=8), dimension(vertexDegree,nVertices), intent(in) :: kiteAreasOnVertex
+
+ character(len=80) :: a, b, c, d, e, f
+ integer :: i, j, k, nVerticesTotal, iEdge, iLoop, iFace, Vert(4), Edge(4), iVertex, i1, i2, jp1
+ integer :: nKitesTotal, iCell, iEdge1, iEdge2, iVertex11, iVertex12, iVertex21, iVertex22, ksave
+ real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xscale, work(nCells), work1(nCells), work2(nCells)
+ real (kind=8) :: xv, yv, zv, xc, yc, zc, dist
+ logical (kind=8) :: eflag
+
+ if(on_a_sphere.eq.'NO ') then
+ write(6,*) ' write_dx, not on a sphere '
+ endif
+
+ xscale = 1.00
+ xCell = xCell*xscale
+ yCell = yCell*xscale
+ zCell = zCell*xscale
+ xVertex = xVertex*xscale
+ yVertex = yVertex*xscale
+ zVertex = zVertex*xscale
+ xEdge = xEdge*xscale
+ yEdge = yEdge*xscale
+ zEdge = zEdge*xscale
+
+ write(6,*) 'xCell', minval(xCell), maxval(xCell)
+ write(6,*) ' nCells', nCells
+ write(6,*) ' nEdges', nEdges
+ write(6,*) ' nVertices', nVertices
+ write(6,*) ' nEdgesOnCell',minval(nEdgesOnCell), maxval(nEdgesOnCell)
+
+ open(unit=1,file='dx/vector.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file vector.position.data')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,*)
+
+ a = trim('object 0 class array type float rank 1 shape 3 items')
+ b = trim('ascii data file vector.data')
+ c = trim('attribute "dep" string "positions"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "vector" class field')
+ b = trim('component "positions" "positions list"')
+ c = trim('component "data" 0')
+ write(1,10) a
+ write(1,10) b
+ write(1,10) c
+
+ close(1)
+
+ open(unit=14,file='dx/vector.position.data',form='formatted',status='unknown')
+ do i=1,nCells
+ write(14,22) xCell(i), yCell(i), zCell(i)
+ enddo
+ close(14)
+
+
+
+ nVerticesTotal = 0
+ do i=1,nCells
+ nVerticesTotal = nVerticesTotal + nEdgesOnCell(i)
+ enddo
+ write(6,*) 'total number of vertices', nVerticesTotal
+
+ open(unit=1,file='dx/ocean.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file ocean.position.data')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,*)
+ 10 format(a70,i10)
+
+ a = trim('object "edge list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.edge.data')
+ c = trim('attribute "ref" string "positions"')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "loops list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.loop.data')
+ c = trim('attribute "ref" string "edges"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "face list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.face.data')
+ c = trim('attribute "ref" string "loops"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object 0 class array type float rank 0 items')
+ b = trim('data file ocean.meshSpacing.data')
+ c = trim('attribute "dep" string "faces"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "area" class field')
+ b = trim('component "positions" "positions list"')
+ c = trim('component "edges" "edge list"')
+ d = trim('component "loops" "loops list"')
+ e = trim('component "faces" "face list"')
+ f = trim('component "data" 0')
+ write(1,10) a
+ write(1,10) b
+ write(1,10) c
+ write(1,10) d
+ write(1,10) e
+ write(1,10) f
+
+ close(1)
+
+
+ work2 = meshSpacing
+ work1 = depthCell
+ work = SST
+
+ open(unit= 8,file='dx/ocean.meshSpacing.data',form='formatted',status='unknown')
+ open(unit= 9,file='dx/ocean.depth.data',form='formatted',status='unknown')
+ open(unit=10,file='dx/ocean.area.data',form='formatted',status='unknown')
+ open(unit=11,file='dx/ocean.face.data',form='formatted',status='unknown')
+ open(unit=12,file='dx/ocean.loop.data',form='formatted',status='unknown')
+ open(unit=13,file='dx/ocean.edge.data',form='formatted',status='unknown')
+ open(unit=14,file='dx/ocean.position.data',form='formatted',status='unknown')
+
+ iLoop = 0
+ iEdge = 0
+ do i=1,nCells
+ write(8,20) work2(i)
+ write(9,20) work1(i)
+ write(10,20) work(i)
+ write(11,21) i-1
+ write(12,21) iLoop
+ iLoop = iLoop + nEdgesOnCell(i)
+
+ eflag = .false.
+ do j=1,nEdgesOnCell(i)
+ k = verticesOnCell(j,i)
+ xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
+ xc = xCell(i); yc = yCell(i); zc = zCell(i)
+ dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
+ if(dist.gt.5.0e5.and.on_a_sphere.eq.'NO ') then
+ eflag = .true.
+ endif
+ enddo
+
+ if(eflag) then
+
+ do j=1,nEdgesOnCell(i)
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ k = verticesOnCell(j,i)
+ xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
+ xc = xCell(i); yc = yCell(i); zc = zCell(i)
+ dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
+ if(dist.gt.5.0e5) then
+ write(14,22) xc, yc, zc
+ else
+ write(14,22) xv, yv, zv
+ endif
+ enddo
+
+ else
+
+ do j=1,nEdgesOnCell(i)
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ k = verticesOnCell(j,i)
+ if(k.le.0) write(6,*) ' vert1 ',k, verticesOnCell(:,i)
+ write(14,22) xVertex(k), yVertex(k), zVertex(k)
+ write(15,23) j,i,k,xVertex(k), yVertex(k), zVertex(k)
+ enddo
+ endif
+ enddo
+
+ 20 format(e20.10)
+ 21 format(i20)
+ 22 format(3e20.10)
+ 23 format(3i8, 3e20.10)
+
+ close(9)
+ close(10)
+ close(11)
+ close(12)
+ close(13)
+ close(14)
+
+ ! nVerticesTotal = 0
+ ! nKitesTotal = 0
+ ! do i=1,nCells
+ ! nKitesTotal = nKitesTotal + nEdgesOnCell(i)
+ ! enddo
+ ! nVerticesTotal = nKitesTotal*4
+ ! write(6,*) nKitesTotal, nVerticesTotal
+
+ ! open(unit=1,file='dx/kite.dx',form='formatted',status='unknown')
+
+ ! a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ ! b = trim('ascii data file kite.position.data')
+ ! write(1,10) a, nVerticesTotal
+ ! write(1,10) b
+ ! write(1,*)
+
+ ! a = trim('object "edge list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.edge.data')
+ ! c = trim('attribute "ref" string "positions"')
+ ! write(1,10) a, nVerticesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "loops list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.loop.data')
+ ! c = trim('attribute "ref" string "edges"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "face list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.face.data')
+ ! c = trim('attribute "ref" string "loops"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object 0 class array type float rank 0 items')
+ ! b = trim('data file kite.area.data')
+ ! c = trim('attribute "dep" string "faces"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "area" class field')
+ ! b = trim('component "positions" "positions list"')
+ ! c = trim('component "edges" "edge list"')
+ ! d = trim('component "loops" "loops list"')
+ ! e = trim('component "faces" "face list"')
+ ! f = trim('component "data" 0')
+ ! write(1,10) a
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,10) d
+ ! write(1,10) e
+ ! write(1,10) f
+
+ ! close(1)
+
+ ! open(unit=10,file='dx/kite.area.data',form='formatted',status='unknown')
+ ! open(unit=11,file='dx/kite.face.data',form='formatted',status='unknown')
+ ! open(unit=12,file='dx/kite.loop.data',form='formatted',status='unknown')
+ ! open(unit=13,file='dx/kite.edge.data',form='formatted',status='unknown')
+ ! open(unit=14,file='dx/kite.position.data',form='formatted',status='unknown')
+
+ ! iLoop = 0
+ ! iEdge = 0
+ ! iFace = 0
+
+ ! do iCell=1,nCells
+ ! do j=1,nEdgesOnCell(iCell)
+ ! iEdge1 = edgesOnCell(j,iCell)
+ ! jp1 = j+1
+ ! if(j.eq.nEdgesOnCell(iCell)) jp1=1
+ ! iEdge2 = edgesOnCell(jp1,iCell)
+
+ ! iVertex11 = verticesOnEdge(1,iEdge1)
+ ! iVertex21 = verticesOnEdge(2,iEdge1)
+ ! iVertex12 = verticesOnEdge(1,iEdge2)
+ ! ivertex22 = verticesOnEdge(2,iEdge2)
+
+ ! if(iVertex11.eq.iVertex12.or.iVertex11.eq.iVertex22) then
+ ! iVertex = iVertex11
+ ! elseif(iVertex21.eq.iVertex12.or.iVertex21.eq.iVertex22) then
+ ! iVertex = iVertex21
+ ! else
+ ! write(6,*) iVertex11, iVertex21, iVertex12, iVertex22
+ ! stop
+ ! endif
+
+ ! ksave = 0
+ ! do k=1,vertexDegree
+ ! if(cellsOnVertex(k,iVertex).eq.iCell) ksave=k
+ ! enddo
+ ! if(ksave.eq.0) then
+ ! write(6,*) ' can not find iCell'
+ ! write(6,*) cellsOnVertex(:,iVertex)
+ ! write(6,*) iCell
+ ! write(6,*) iEdge1, iEdge2
+ ! write(6,*) iVertex11, iVertex21, iVertex21, iVertex22
+ ! write(6,*) iVertex
+ ! stop
+ ! endif
+
+ ! write(11,21) iFace
+ ! write(12,21) iLoop
+ ! iFace = iFace + 1
+ ! iLoop = iLoop + 4
+ ! do k=1,4
+ ! write(13,21) iEdge
+ ! iEdge = iEdge + 1
+ ! enddo
+ !
+ ! x1 = xCell(iCell) ; y1 = yCell(iCell) ; z1 = zCell(iCell)
+ ! x2 = xEdge(iEdge1) ; y2 = yEdge(iEdge1) ; z2 = zEdge(iEdge1)
+ ! x3 = xVertex(iVertex); y3 = yVertex(iVertex); z3 = zVertex(iVertex)
+ ! x4 = xEdge(iEdge2) ; y4 = yEdge(iEdge2) ; z4 = zEdge(iEdge2)
+ !
+ ! write(14,22) x1, y1, z1
+ ! write(14,22) x2, y2, z2
+ ! write(14,22) x3, y3, z3
+ ! write(14,22) x4, y4, z4
+ ! write(10,22) kiteAreasOnVertex(ksave,iVertex)
+
+ ! enddo
+ ! enddo
+
+end subroutine write_OpenDX
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE CONVERT_LX
+!
+! Convert (lat,lon) to an (x, y, z) location on a sphere with specified radius.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine convert_lx(x, y, z, radius, lat, lon)
+
+ implicit none
+
+ real, intent(in) :: radius
+ real, intent(in) :: lat, lon
+ real, intent(out) :: x, y, z
+
+ z = radius * sin(lat)
+ x = radius * cos(lon) * cos(lat)
+ y = radius * sin(lon) * cos(lat)
+
+end subroutine convert_lx
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE CONVERT_XL
+!
+! Convert (x, y, z) to a (lat, lon) location on a sphere with
+! radius sqrt(x^2 + y^2 + z^2).
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine convert_xl(x, y, z, lat,lon)
+
+ implicit none
+
+ real, intent(in) :: x, y, z
+ real, intent(out) :: lat, lon
+
+ real :: dl, clat, pii, rtod
+ real :: eps
+ parameter (eps=1.e-10)
+
+ pii = 2.*asin(1.0)
+ rtod=180./pii
+ dl = sqrt(x*x + y*y + z*z)
+
+ lat = asin(z/dl)
+
+! check for being close to either pole
+
+ if (abs(x) > eps) then
+
+ if (abs(y) > eps) then
+
+ lon = atan(abs(y/x))
+
+ if ((x <= 0.) .and. (y >= 0.)) then
+ lon = pii-lon
+ else if ((x <= 0.) .and. (y < 0.)) then
+ lon = lon+pii
+ else if ((x >= 0.) .and. (y <= 0.)) then
+ lon = 2*pii-lon
+ end if
+
+ else ! we're either on longitude 0 or 180
+
+ if (x > 0) then
+ lon = 0.
+ else
+ lon = pii
+ end if
+
+ end if
+
+ else if (abs(y) > eps) then
+
+ if (y > 0) then
+ lon = pii/2.
+ else
+ lon = 3.*pii/2.
+ end if
+
+ else ! we are at a pole
+
+ lon = 0.
+
+ end if
+
+end subroutine convert_xl
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine transform_from_lonlat_to_xyz(xin, yin, zin, ulon, ulat, ux, uy, uz)
+!
+! transform vector measured in latitude/longitude space to a vector measured in x,y,z
+!
+! INTENT(IN)
+! xin = x position
+! yin = y position
+! zin = z position
+! ulon = east component of vector
+! ulat = north component of vector
+!
+! INTENT(OUT)
+! ux = x component of vector
+! uy = y component of vector
+! uz = z component of vector
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+real, intent(in) :: xin, yin, zin, ulon, ulat
+real, intent(out) :: ux, uy, uz
+real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X2_to_X1(3,3), r
+integer :: i,j,k
+logical :: l_Pole
+real, parameter :: epsvt = 1.0e-10
+
+!-----------------------------------------------------------------------
+! define the e1, e2, and e3 directions
+!-----------------------------------------------------------------------
+ X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
+ X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
+ X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
+
+!-----------------------------------------------------------------------
+! find the vectors (measured in X1) that point in the local
+! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
+!-----------------------------------------------------------------------
+ h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
+ call unit_vector_in_3space(h(3,:))
+
+!-----------------------------------------------------------------------
+! g(:) is a work array and holds the vector pointing to the North Pole.
+! measured in X1
+!-----------------------------------------------------------------------
+ g(:) = X1(3,:)
+
+!-----------------------------------------------------------------------
+! determine if the local vertical hits a pole
+!-----------------------------------------------------------------------
+ l_Pole = .false.
+ r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
+ r = abs(r) + epsvt
+ if(r.gt.1.0) then
+ l_Pole = .true.
+ h(3,:) = h(3,:) + epsvt
+ call unit_vector_in_3space(h(3,:))
+ endif
+
+!-----------------------------------------------------------------------
+! find the vector that is perpendicular to the local vertical vector
+! and points in the direction of of the North pole, this defines the local
+! north direction. measured in X1
+!-----------------------------------------------------------------------
+ call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
+
+!-----------------------------------------------------------------------
+! take the cross product of the local North direction and the local vertical
+! to find the local east vector. still in X1
+!-----------------------------------------------------------------------
+ call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
+
+!-----------------------------------------------------------------------
+! put these 3 vectors into a matrix X2
+!-----------------------------------------------------------------------
+ X2(1,:) = h(1,:) ! local east (measured in X1)
+ X2(2,:) = h(2,:) ! local north (measured in X1)
+ X2(3,:) = h(3,:) ! local vertical (measured in X1)
+
+!-----------------------------------------------------------------------
+! compute the transformation matrix
+!-----------------------------------------------------------------------
+ trans_X2_to_X1(:,:) = matmul(X1,transpose(X2))
+
+!-----------------------------------------------------------------------
+! transform (ulon, ulat) into (x,y,z)
+!-----------------------------------------------------------------------
+ p(1) = ulon; p(2) = ulat; p(3) = 0
+ g(:) = matmul(trans_X2_to_X1(:, :), p(:))
+ ux = g(1); uy = g(2); uz = g(3)
+
+end subroutine transform_from_lonlat_to_xyz
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine transform_from_xyz_to_lonlat(xin, yin, zin, ux, uy, uz, ulon, ulat)
+!
+! transform vector measured in x,y,z space to a vector measured in latitude/longitude space
+!
+! INTENT(IN)
+! xin = x position
+! yin = y position
+! zin = z position
+! ux = x component of vector
+! uy = y component of vector
+! uz = z component of vector
+!
+! INTENT(OUT)
+! ulon = east component of vector
+! ulat = north component of vector
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+real, intent(in) :: xin, yin, zin, ux, uy, uz
+real, intent(out) :: ulon, ulat
+real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X1_to_X2(3,3), r
+integer :: i,j,k
+logical :: l_Pole
+real, parameter :: epsvt = 1.0e-10
+
+!-----------------------------------------------------------------------
+! define the e1, e2, and e3 directions
+!-----------------------------------------------------------------------
+ X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
+ X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
+ X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
+
+!-----------------------------------------------------------------------
+! find the vectors (measured in X1) that point in the local
+! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
+!-----------------------------------------------------------------------
+ h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
+ call unit_vector_in_3space(h(3,:))
+
+!-----------------------------------------------------------------------
+! g(:) is a work array and holds the vector pointing to the North Pole.
+! measured in X1
+!-----------------------------------------------------------------------
+ g(:) = X1(3,:)
+
+!-----------------------------------------------------------------------
+! determine if the local vertical hits a pole
+!-----------------------------------------------------------------------
+ l_Pole = .false.
+ r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
+ r = abs(r) + epsvt
+ if(r.gt.1.0) then
+ l_Pole = .true.
+ h(3,:) = h(3,:) + epsvt
+ call unit_vector_in_3space(h(3,:))
+ endif
+
+!-----------------------------------------------------------------------
+! find the vector that is perpendicular to the local vertical vector
+! and points in the direction of of the North pole, this defines the local
+! north direction. measured in X1
+!-----------------------------------------------------------------------
+ call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
+
+!-----------------------------------------------------------------------
+! take the cross product of the local North direction and the local vertical
+! to find the local east vector. still in X1
+!-----------------------------------------------------------------------
+ call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
+
+!-----------------------------------------------------------------------
+! put these 3 vectors into a matrix X2
+!-----------------------------------------------------------------------
+ X2(1,:) = h(1,:) ! local east (measured in X1)
+ X2(2,:) = h(2,:) ! local north (measured in X1)
+ X2(3,:) = h(3,:) ! local vertical (measured in X1)
+
+!-----------------------------------------------------------------------
+! compute the transformation matrix
+!-----------------------------------------------------------------------
+ trans_X1_to_X2(:,:) = matmul(X2,transpose(X1))
+
+!-----------------------------------------------------------------------
+! transform (ulon, ulat) into (x,y,z)
+!-----------------------------------------------------------------------
+ p(1) = ux; p(2) = uy; p(3) = uz
+ g(:) = matmul(trans_X1_to_X2(:, :), p(:))
+ ulon = g(1); ulat= g(2);
+
+end subroutine transform_from_xyz_to_lonlat
+
+!======================================================================
+! BEGINNING OF UNIT_VECTOR_IN_3SPACE
+!======================================================================
+ subroutine unit_vector_in_3space (p_1)
+
+!-----------------------------------------------------------------------
+! PURPOSE : normalize p_1 to unit length and overwrite p_1
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(inout)
+!-----------------------------------------------------------------------
+ real , intent(inout) :: &
+ p_1 (:)
+
+!-----------------------------------------------------------------------
+! local
+!-----------------------------------------------------------------------
+ real :: length
+
+ length = SQRT (p_1(1)**2 + p_1(2)**2 + p_1(3)**2 )
+ length = 1.0/length
+ p_1(1) = p_1(1)*length
+ p_1(2) = p_1(2)*length
+ p_1(3) = p_1(3)*length
+
+ end subroutine unit_vector_in_3space
+!======================================================================
+! END OF UNIT_VECTOR_IN_3SPACE
+!======================================================================
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_3SPACE
+!======================================================================
+ subroutine cross_product_in_3space(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real , intent(in) :: &
+ p_1 (:), &
+ p_2 (:)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real , intent(out) :: &
+ p_out (:)
+
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+ end subroutine cross_product_in_3space
+!======================================================================
+! END OF CROSS_PRODUCT_IN_3SPACE
+!======================================================================
+
+!======================================================================
+! BEGINNING OF VECTOR_ON_TANGENT_PLANE
+!======================================================================
+ subroutine vector_on_tangent_plane(p_1, p_2, p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE : given two points measured in (x,y,z) and lying on
+! the unit sphere, find the vector (p_out) that lies on the plane
+! perpendicular to the p_1 vector and points in the direction of
+! the projection of p_2 onto the tangent plane.
+!
+! NOTE : p_1 and p_2 are assumed to be of unit length
+! NOTE : p_out is normalized to unit length
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real , intent(in) :: &
+ p_1 (:), &
+ p_2 (:)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real , intent(out) :: &
+ p_out (:)
+
+!-----------------------------------------------------------------------
+! local
+!-----------------------------------------------------------------------
+ real :: &
+ work (3), t1(3), t2(3)
+
+! work (1) = - p_1(2) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
+! + p_1(3) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) )
+
+! work (2) = + p_1(1) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
+! - p_1(3) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
+
+! work (3) = - p_1(1) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) ) &
+! + p_1(2) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
+
+
+ t1(:) = p_2(:) - p_1(:)
+ t2(:) = p_1
+
+ call unit_vector_in_3space (t1)
+ call unit_vector_in_3space (t2)
+
+ call cross_product_in_3space(t1(:), t2(:), work(:))
+ call unit_vector_in_3space (work)
+ call cross_product_in_3space(t2(:),work(:),p_out(:))
+ call unit_vector_in_3space (p_out)
+
+ end subroutine vector_on_tangent_plane
+!======================================================================
+! END OF VECTOR_ON_TANGENT_PLANE
+!======================================================================
+
+end module utilities
</font>
</pre>