<p><b>sprice@lanl.gov</b> 2011-08-19 08:55:41 -0600 (Fri, 19 Aug 2011)</p><p>BRANCH COMMIT (land ice): More updates to code for extracting ice sheet domain from periodic hex grid. Partially working for Greenland.<br>
</p><hr noshade><pre><font color="gray">Deleted: branches/land_ice/icesheet/src/icesheet.F
===================================================================
--- branches/land_ice/icesheet/src/icesheet.F        2011-08-18 23:07:37 UTC (rev 949)
+++ branches/land_ice/icesheet/src/icesheet.F        2011-08-19 14:55:41 UTC (rev 950)
@@ -1,1314 +0,0 @@
-program map_to_icesheet 
-
-use read_netcdf
-use read_topo
-use read_TS
-use write_netcdf
-use utilities
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Program: basin.F
-!
-! This program is meant to add land to grids, as well as initial conditions.
-! 
-! This program is used to take a specific mesh, and remove Cells from it
-! It can be used to change a planar grid into a Channel or a basin grid, or to 
-! Change a spherical grid into a Limited area spherical grid.
-!
-! How to use:
-! Step 1: Set the number of Vertical levels
-! Step 2: Set if the grid is on a sphere or not, and it's radius
-! Step 3: Specify some Parameters
-! Step 4: Check the Initial conditions routine get_init_conditions
-! Step 5: Check the depth routine define_kmt
-!
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-implicit none
-
-integer, parameter :: nx = 50
-real, parameter :: Lx = 2.0e6
-
-! 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
-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
-
-real, allocatable, dimension(:) :: fEdge, fVertex, h_s, work1
-real, allocatable, dimension(:,:) :: u_src
-real, allocatable, dimension(:,:,:) :: u, v, h
-real, allocatable, dimension(:,:,:) :: rho
-
-integer nlon, nlat, ndepth
-real(kind=4), allocatable, dimension(:) :: t_lon, t_lat, depth_t
-real(kind=4), allocatable, dimension(:,:) :: mTEMP, mSALT
-real(kind=4), allocatable, dimension(:,:,:) :: TEMP, SALT
-real(kind=4), allocatable, dimension(:,:) :: TAUX, TAUY
-
-real, dimension(1) :: dz
-
-! Step 1: Set the number of Vertical levels
-!integer, parameter :: nVertLevelsMOD = 40
-integer, parameter :: nVertLevelsMOD = 1
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! basin-mod
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-! Need to uncomment the options appropriate for the input grid file. If it's on
-! a sphere, specify the flag &quot;on_a_sphere&quot; and &quot;sphere_radius&quot;. Otherwise set
-! them equal to NO and 0.0 respectively
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-! Step 2: Set if the grid is on a sphere or not, and it's radius
-!character (len=16) :: on_a_sphere = 'YES              '
-!real*8, parameter :: sphere_radius = 6.37122e6
-!!real*8, parameter :: sphere_radius = 1.0
-
-character (len=16) :: on_a_sphere = 'NO               '
-real*8, parameter :: sphere_radius = 0.0
-
-logical, parameter :: real_bathymetry=.false.
-logical, parameter :: l_woce = .false.
-
-
-! Step 3: Specify some Parameters
-   real (kind=8), parameter :: &amp;
-    h_total_max = 2000.0, &amp;
-    u_max = 0.0, &amp;
-    u_src_max = 0.1, &amp; ! max wind stress, N/m2
-    beta = 1.4e-11, &amp;
-    f0 = -1.1e-4, &amp;
-    omega = 7.29212e-5
-
-   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
-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
-
-real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew
-real, allocatable, dimension(:,:) :: u_srcNew
-real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew
-real, allocatable, dimension(:,:,:) :: rhoNew, temperatureNew, salinityNew, tracer1New
-real, allocatable, dimension(:) :: temperatureRestoreNew, salinityRestoreNew
-
-! mapping variables
-integer, allocatable, dimension(:) :: kmt, maxLevelCellNew
-integer, allocatable, dimension(:) :: cellMap, edgeMap, vertexMap
-real, allocatable, dimension(:) :: depthCell
-
-! 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
-
-! get to work
-write(6,*) ' starting'
-write(6,*)
-
-! get depth profile for later
-call get_dz
-
-! get grid
-write(6,*) ' calling read_grid'
-write(6,*)
-call read_grid
-write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
-
-! 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
-call error_checking
-
-write(6,*) ' getting woce t and s '
-call read_TS_init(nlon, nlat, ndepth)
-write(6,*) ' TS INIT ', nlon, nlat, ndepth
-allocate(t_lon(nlon), t_lat(nlat), depth_t(ndepth), TEMP(nlon,nlat,ndepth), SALT(nlon,nlat,ndepth))
-allocate(TAUX(nlon,nlat), TAUY(nlon,nlat))
-allocate(mTEMP(nlat,ndepth), mSALT(nlat,ndepth))
-call read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
-call read_TS_finalize()
-do k=1,ndepth
-     ndata = 0; temp_t=0; temp_s=0
-     do j=1,nlat
-     do i=1,nlon
-       if(TEMP(i,j,k).gt.-10.0) then
-         ndata = ndata + 1
-         temp_t = temp_t + TEMP(i,j,k)
-         temp_s = temp_s + SALT(i,j,k)
-       endif
-     enddo
-     enddo
-     mTEMP(:,k) = temp_t / float(ndata)
-     mSALT(:,k) = temp_s / float(ndata)
-     write(6,*) ndata,mTemp(1,k),mSalt(1,k)
-enddo
-
-! generate initial conditions
-call get_init_conditions
-
-! 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
-write(6,*) ' calling write_OpenDX'
-write(6,*)
-call write_OpenDX(           on_a_sphere, &amp;
-                             nCellsNew, &amp;
-                             nVerticesNew, &amp;
-                             nEdgesNew, &amp;
-                             vertexDegreeNew, &amp;
-                             maxEdgesNew, &amp;
-                             xCellNew, &amp;
-                             yCellNew, &amp;
-                             zCellNew, &amp;
-                             xVertexNew, &amp;
-                             yVertexNew, &amp;
-                             zVertexNew, &amp;
-                             xEdgeNew, &amp;
-                             yEdgeNew, &amp;
-                             zEdgeNew, &amp;
-                             nEdgesOnCellNew, &amp;
-                             verticesOnCellNew, &amp;
-                             verticesOnEdgeNew, &amp;
-                             cellsOnVertexNew, &amp;
-                             edgesOnCellNew, &amp;
-                             areaCellNew, &amp;
-                             maxLevelCellNew, &amp;
-                             depthCell, &amp;
-                             kiteAreasOnVertexNew )
-
-
-!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
-
-
-!Step 4: Check the Initial conditions routine get_init_conditions
-subroutine get_init_conditions
-implicit none
-real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
-real :: dotProd
-integer :: iTracer, ix, iy, ndata, i, j, k, ixt, iyt, ncull, jcount, iNoData, kdata(nVertLevelsMod)
-logical :: flag_lat
-
-pi = 4.0*atan(1.0)
-dtr = pi/180.0
-
-hNew = 100.0
-temperatureNew = 1.0
-salinityNew = 1.0
-tracer1New = 1.0
-uNew = 0
-vNew = 0
-
-write(6,*) ' done get_init_conditions'
-
-end subroutine get_init_conditions
-
-
-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,&amp;
-                       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(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(fEdge(nEdges))
-allocate(fVertex(nVertices))
-allocate(h_s(nCells))
-allocate(work1(nCells))
-allocate(u_src(nVertLevels,nEdges))
-allocate(u(1,nVertLevels,nEdges))
-allocate(v(1,nVertLevels,nEdges))
-allocate(h(1,nVertLevels,nCells))
-allocate(rho(1,nVertLevels,nCells))
-
-xCell=0; yCell=0; zCell=0; latCell=0; lonCell=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
-
-fEdge=0; fVertex=0; h_s=0; u_src=0; work1=0
-u=0; v=0; h=0; rho=0
-
-
-call  read_netcdf_fields( &amp;
-                    time, &amp;
-                    latCell, &amp;
-                    lonCell, &amp;
-                    xCell, &amp;
-                    yCell, &amp;
-                    zCell, &amp;
-                    indexToCellID, &amp;
-                    latEdge, &amp;
-                    lonEdge, &amp;
-                    xEdge, &amp;
-                    yEdge, &amp;
-                    zEdge, &amp;
-                    indexToEdgeID, &amp;
-                    latVertex, &amp;
-                    lonVertex, &amp;
-                    xVertex, &amp;
-                    yVertex, &amp;
-                    zVertex, &amp;
-                    indexToVertexID, &amp;
-                    cellsOnEdge, &amp;
-                    nEdgesOnCell, &amp;
-                    nEdgesOnEdge, &amp;
-                    edgesOnCell, &amp;
-                    edgesOnEdge, &amp;
-                    weightsOnEdge, &amp;
-                    dvEdge, &amp;
-                    dcEdge, &amp;
-                    angleEdge, &amp;
-                    areaCell, &amp;
-                    areaTriangle, &amp;
-                    cellsOnCell, &amp;
-                    verticesOnCell, &amp;
-                    verticesOnEdge, &amp;
-                    edgesOnVertex, &amp;
-                    cellsOnVertex, &amp;
-                    kiteAreasOnVertex, &amp;
-                    fEdge, &amp;
-                    fVertex, &amp;
-                    h_s, &amp;
-                    u, &amp;
-                    v, &amp;
-                    h &amp;
-                   )
-
-write(6,*) ' values from read grid, min/max'
-write(6,*) ' latCell : ', minval(latCell), maxval(latCell)
-write(6,*) ' lonCell : ', minval(lonCell), maxval(lonCell)
-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,*) ' fEdge : ', minval(fEdge), maxval(fEdge)
-write(6,*) ' fVertex : ', minval(fVertex), maxval(fVertex)
-write(6,*) ' h_s : ', minval(h_s), maxval(h_s)
-write(6,*) ' u : ', minval(u), maxval(u)
-write(6,*) ' v : ', minval(v), maxval(v)
-write(6,*) ' h : ', minval(h), maxval(h)
-
-end subroutine read_grid
-
-
-subroutine write_grid
-implicit none
-
-if(on_a_sphere.eq.'YES              ') 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( &amp;
-                nCellsNew, &amp;
-                nEdgesNew, &amp;
-                nVerticesNew, &amp;
-                maxEdgesNew, &amp;
-                nVertLevelsNew, &amp;
-                vertexDegreeNew, &amp;
-                sphere_radius, &amp;
-                on_a_sphere &amp;
-                )
-
-call write_netcdf_fields( &amp;
-                    1, &amp;
-                    latCellNew, &amp;
-                    lonCellNew, &amp;
-                    xCellNew, &amp;
-                    yCellNew, &amp;
-                    zCellNew, &amp;
-                    indexToCellIDNew, &amp;
-                    latEdgeNew, &amp;
-                    lonEdgeNew, &amp;
-                    xEdgeNew, &amp;
-                    yEdgeNew, &amp;
-                    zEdgeNew, &amp;
-                    indexToEdgeIDNew, &amp;
-                    latVertexNew, &amp;
-                    lonVertexNew, &amp;
-                    xVertexNew, &amp;
-                    yVertexNew, &amp;
-                    zVertexNew, &amp;
-                    indexToVertexIDNew, &amp;
-                    maxLevelCellNew, &amp;
-                    cellsOnEdgeNew, &amp;
-                    nEdgesOnCellNew, &amp;
-                    nEdgesOnEdgeNew, &amp;
-                    edgesOnCellNew, &amp;
-                    edgesOnEdgeNew, &amp;
-                    weightsOnEdgeNew, &amp;
-                    dvEdgeNew, &amp;
-                    dcEdgeNew, &amp;
-                    angleEdgeNew, &amp;
-                    areaCellNew, &amp;
-                    areaTriangleNew, &amp;
-                    cellsOnCellNew, &amp;
-                    verticesOnCellNew, &amp;
-                    verticesOnEdgeNew, &amp;
-                    edgesOnVertexNew, &amp;
-                    cellsOnVertexNew, &amp;
-                    kiteAreasOnVertexNew, &amp;
-                    fEdgeNew, &amp;
-                    fVertexNew, &amp;
-                    h_sNew, &amp;
-                    boundaryEdgeNew, &amp;
-                    boundaryVertexNew, &amp;
-                    u_srcNew, &amp;
-                    uNew, &amp;
-                    vNew, &amp;
-                    hNew, &amp;
-                    rhoNew, &amp;
-                    temperatureNew, &amp;
-                    salinityNew, &amp;
-                    tracer1New, &amp;
-                    temperatureRestoreNew, &amp;
-                    salinityRestoreNew &amp;
-                   )
-
-call write_netcdf_finalize
-
-if(on_a_sphere.eq.'YES              ') 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 5: Check the depth routine define_kmt
-subroutine define_kmt
-implicit none
-real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
-real (kind=4), allocatable, dimension(:,:) :: thck  
-integer :: nx, ny, inx, iny, ix, iy
-real :: pi, dtr, thkdata, rlon, rlat, r, ymin, ymax, xmin, xmax
-real :: latmin, latmax, lonmin, lonmax
-logical :: flag, kmt_flag
-pi = 4.0*atan(1.0)
-dtr = pi / 180.0
-
-allocate(kmt(nCells))
-kmt = 0
-
-!    nx = 301        ! 5 km  
-!    ny = 561
-    nx = 76         ! 20 km  
-    ny = 141
-    allocate(x(nx))
-    allocate(y(ny))
-    x = 0.0
-    y = 0.0
-
-    allocate(thck(nx,ny))
-    thck = 0.0
-
-    call read_topo_init( inx, iny)
-    write(6,*) ' inx, iny ', inx, iny                     
-    if(inx.ne.nx) stop ' nx topo'
-    if(iny.ne.ny) stop ' ny topo'
-
-    call read_topo_fields(x,y,thck)
-
-    call read_topo_finalize()
-
-    !!! debugging output !!!
-    write(6,*) ' '
-    print *, 'min(x), max(x)'
-    write(6,*) minval(x), maxval(x)
-    print *, 'min(y), max(y)'
-    write(6,*) minval(y), maxval(y)
-    print *, 'x(1), y(1)'
-    write(6,*) x(1), y(1)
-    write(6,*) ' '
-
-    print *, 'min(xCell), max(xCell)' 
-    write(6,*) minval(xCell), maxval(xCell)
-    print *, 'min(yCell), max(yCell)' 
-    write(6,*) minval(yCell), maxval(yCell)
-    print *, 'xCell(1), yCell(1)'
-    write(6,*) xCell(1), yCell(1)
-    write(6,*) ' '
-
-!    pause
-
-    do iCell=1,nCells
-
-        ix = nint( xCell(iCell) / 20.0d3 )
-        iy = nint( yCell(iCell) / 20.0d3 )
-!        ix = nint( xCell(iCell) / 5.0d3 )
-!        iy = nint( yCell(iCell) / 5.0d3 )
-
-        ix = min(nx,ix); iy=min(ny,iy)
-        ix = max(1,ix); iy=max(1,iy)
-
-        print *, '- new cell -'
-        print *, 'xCell(iCell)', xCell(iCell)
-        print *, 'yCell(iCell)', yCell(iCell)
-        print *, 'ix,iy = ', ix,iy 
-
-        thkdata = thck(ix,iy)
-        print *, 'x1(ix,iy) = ', x(ix) 
-        print *, 'y1(ix,iy) = ', y(iy) 
-        print *, 'thck(ix,iy) = ', thkdata
-
-         if( thkdata .gt. 0 )then
-           kmt(iCell) = nVertLevelsMod
-           print *, 'cell kept'
-           print *, ' '
-         else
-           kmt(iCell) = 0
-           print *, 'cell removed', ix, iy
-           print *, ' '
-         endif
-
-    enddo
-
-    ! output no. of cells removed
-    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)
-
-    deallocate(x)
-    deallocate(y)
-    deallocate(thck)
-
-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(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))
-allocate(maxLevelCellNew(nCellsNew))
-allocate(depthCell(nCellsNew))
-
-allocate(fEdgeNew(nEdgesNew))
-allocate(fVertexNew(nVerticesNew))
-allocate(h_sNew(nCellsNew))
-allocate(u_srcNew(nVertLevelsNew,nEdgesNew))
-allocate(uNew(1,nVertLevelsNew,nEdgesNew))
-allocate(vNew(1,nVertLevelsNew,nEdgesNew))
-allocate(hNew(1,nVertLevelsNew,nCellsNew))
-allocate(rhoNew(1,nVertLevelsNew,nCellsNew))
-allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
-allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
-allocate(tracer1New(1,nVertLevelsNew,nCellsNew))
-
-allocate(temperatureRestoreNew(nCellsNew))
-allocate(salinityRestoreNew(nCellsNew))
-
-xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0
-xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
-xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
-
-fEdgeNew=0; fVertexNew=0; h_sNew=0; u_srcNew=0
-uNew=0; vNew=0; hNew=0; rhoNew=0
-temperatureNew=0; salinityNew=0; tracer1New=0;
-
-temperatureRestoreNew = 0.0
-salinityRestoreNew = 0.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)
-    areaCellNew(jNew)=areaCell(i)
-    maxLevelCellNew(jNew) = kmt(i)
-    depthCell(jNew) = kmt(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)
-    fEdgeNew(jNew) = fEdge(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)
-    fVertexNew(jNew)=fVertex(i)
-    areaTriangleNew(jNew)=areaTriangle(i)
-endif
-enddo
-
-deallocate(xCell)
-deallocate(yCell)
-deallocate(zCell)
-deallocate(latCell)
-deallocate(lonCell)
-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 &quot;cellsOnEdge&quot;
-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) then
-    write(6,*) iEdge
-    write(6,*) iEdgeNew
-    write(6,*) iVertex1New, iVertex1
-    write(6,*) iVertex2New, iVertex2
-    stop &quot;verticesOnEdge&quot;
-endif
-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 &quot;jEdgeNew&quot;
-    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
-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 &quot;edgesOnCell&quot;
-edgesOnCellNew(i,iCellNew) = jNew
-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 &quot;verticesOnCell&quot;
-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
-
-
-subroutine get_dz
-integer k
-
-  dz( 1) = 10.0   !   5.006218       10.01244
-!  dz( 1) = 1001.244   !   5.006218       10.01244
-!  dz( 2) = 1011.258   !   15.06873       20.12502
-!  dz( 3) = 1031.682   !   25.28342       30.44183
-!  dz( 4) = 1063.330   !   35.75848       41.07513
-!  dz( 5) = 1107.512   !   46.61269       52.15025
-!  dz( 6) = 1166.145   !   57.98098       63.81171
-!  dz( 7) = 1241.928   !   70.02135       76.23099
-!  dz( 8) = 1338.612   !   82.92405       89.61711
-!  dz( 9) = 1461.401   !   96.92412       104.2311
-!  dz(10) = 1617.561   !   112.3189       120.4067
-!  dz(11) = 1817.368   !   129.4936       138.5804
-!  dz(12) = 2075.558   !   148.9582       159.3360
-!  dz(13) = 2413.680   !   171.4044       183.4728
-!  dz(14) = 2863.821   !   197.7919       212.1110
-!  dz(15) = 3474.644   !   229.4842       246.8575
-!  dz(16) = 4320.857   !   268.4617       290.0660
-!  dz(17) = 5516.812   !   317.6501       345.2342
-!  dz(18) = 7230.458   !   381.3865       417.5388
-!  dz(19) = 9674.901   !   465.9133       514.2878
-!  dz(20) = 13003.92   !   579.3074       644.3270
-!  dz(21) = 17004.89   !   729.3514       814.3759
-!  dz(22) = 20799.33   !   918.3725       1022.369
-!  dz(23) = 23356.94   !   1139.154       1255.939
-!  dz(24) = 24527.19   !   1378.574       1501.210
-!  dz(25) = 24898.04   !   1625.701       1750.191
-!  dz(26) = 24983.22   !   1875.107       2000.023
-!  dz(27) = 24997.87   !   2125.012       2250.002
-!  dz(28) = 24999.79   !   2375.000       2500.000
-!  dz(29) = 24999.98   !   2625.000       2749.999
-!  dz(30) = 25000.00   !   2874.999       2999.999
-!  dz(31) = 25000.00   !   3124.999       3249.999
-!  dz(32) = 25000.00   !   3374.999       3499.999
-!  dz(33) = 25000.00   !   3624.999       3749.999
-!  dz(34) = 25000.00   !   3874.999       3999.999
-!  dz(35) = 25000.00   !   4124.999       4249.999
-!  dz(36) = 25000.00   !   4374.999       4499.999
-!  dz(37) = 25000.00   !   4624.999       4749.999
-!  dz(38) = 25000.00   !   4874.999       4999.999
-!  dz(39) = 25000.00   !   5124.999       5249.999
-!  dz(40) = 25000.00   !   5374.999       5499.999
-
-!  dz = dz / 100.0
-
-  write(6,*)
-!  do k=1,40
-  do k=1,1
-    print *, 'k, dz = ' 
-    write(6,*) k,dz(k)
-  enddo
-  write(6,*)
-
-end subroutine get_dz
-end program map_to_icesheet

Added: branches/land_ice/icesheet/src/icesheet.F
===================================================================
--- branches/land_ice/icesheet/src/icesheet.F                                (rev 0)
+++ branches/land_ice/icesheet/src/icesheet.F        2011-08-19 14:55:41 UTC (rev 950)
@@ -0,0 +1,1379 @@
+program map_to_icesheet 
+
+use read_netcdf
+use read_topo
+use read_TS
+use write_netcdf
+use utilities
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Program: basin.F
+!
+! This program is meant to add land to grids, as well as initial conditions.
+! 
+! This program is used to take a specific mesh, and remove Cells from it
+! It can be used to change a planar grid into a Channel or a basin grid, or to 
+! Change a spherical grid into a Limited area spherical grid.
+!
+! How to use:
+! Step 1: Set the number of Vertical levels
+! Step 2: Set if the grid is on a sphere or not, and it's radius
+! Step 3: Specify some Parameters
+! Step 4: Check the Initial conditions routine get_init_conditions
+! Step 5: Check the depth routine define_kmt
+!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+implicit none
+
+integer, parameter :: nx = 50
+real, parameter :: Lx = 2.0e6
+
+! 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
+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
+
+real, allocatable, dimension(:) :: fEdge, fVertex, h_s, work1
+real, allocatable, dimension(:,:) :: u_src
+real, allocatable, dimension(:,:,:) :: u, v, h
+real, allocatable, dimension(:,:,:) :: rho
+
+! ice sheet variables
+real(kind=4), allocatable, dimension(:,:) :: thck, topg
+
+integer nlon, nlat, ndepth
+real(kind=4), allocatable, dimension(:) :: t_lon, t_lat, depth_t
+real(kind=4), allocatable, dimension(:,:) :: mTEMP, mSALT
+real(kind=4), allocatable, dimension(:,:,:) :: TEMP, SALT
+real(kind=4), allocatable, dimension(:,:) :: TAUX, TAUY
+
+
+real, dimension(1) :: dz
+
+! Step 1: Set the number of Vertical levels
+!integer, parameter :: nVertLevelsMOD = 40
+integer, parameter :: nVertLevelsMOD = 1
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! basin-mod
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! Need to uncomment the options appropriate for the input grid file. If it's on
+! a sphere, specify the flag &quot;on_a_sphere&quot; and &quot;sphere_radius&quot;. Otherwise set
+! them equal to NO and 0.0 respectively
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! Step 2: Set if the grid is on a sphere or not, and it's radius
+!character (len=16) :: on_a_sphere = 'YES              '
+!real*8, parameter :: sphere_radius = 6.37122e6
+!!real*8, parameter :: sphere_radius = 1.0
+
+character (len=16) :: on_a_sphere = 'NO               '
+real*8, parameter :: sphere_radius = 0.0
+
+logical, parameter :: real_bathymetry=.false.
+logical, parameter :: l_woce = .false.
+
+
+! Step 3: Specify some Parameters
+   real (kind=8), parameter :: &amp;
+    h_total_max = 2000.0, &amp;
+    u_max = 0.0, &amp;
+    u_src_max = 0.1, &amp; ! max wind stress, N/m2
+    beta = 1.4e-11, &amp;
+    f0 = -1.1e-4, &amp;
+    omega = 7.29212e-5
+
+   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
+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
+
+real, allocatable, dimension(:) :: fEdgeNew, fVertexNew, h_sNew
+real, allocatable, dimension(:,:) :: u_srcNew
+real, allocatable, dimension(:,:,:) :: uNew, vNew, hNew
+real, allocatable, dimension(:,:,:) :: rhoNew, temperatureNew, salinityNew, tracer1New
+real, allocatable, dimension(:) :: temperatureRestoreNew, salinityRestoreNew
+
+!ice sheet variables
+real(kind=8), allocatable, dimension(:,:,:) :: thckNew, topgNew  
+!real(kind=4), allocatable, dimension(:,:,:) :: thckNew0, topgNew0   ! temp work variables - figure out how to do w/o these  
+real(kind=4), allocatable, dimension(:) :: thckdata, topgdata
+
+! mapping variables
+integer, allocatable, dimension(:) :: cellMap, edgeMap, vertexMap
+real, allocatable, dimension(:) :: depthCell
+integer, allocatable, dimension(:) :: kmt, maxLevelCellNew 
+
+! 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
+
+! get to work
+write(6,*) ' starting'
+write(6,*)
+
+! get depth profile for later
+call get_dz
+
+! get grid
+write(6,*) ' calling read_grid'
+write(6,*)
+call read_grid
+write(6,*) ' xCell 1: ',minval(xCell), maxval(xCell)
+
+! 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 from 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
+call error_checking
+
+write(6,*) ' getting woce t and s '
+call read_TS_init(nlon, nlat, ndepth)
+write(6,*) ' TS INIT ', nlon, nlat, ndepth
+allocate(t_lon(nlon), t_lat(nlat), depth_t(ndepth), TEMP(nlon,nlat,ndepth), SALT(nlon,nlat,ndepth))
+allocate(TAUX(nlon,nlat), TAUY(nlon,nlat))
+allocate(mTEMP(nlat,ndepth), mSALT(nlat,ndepth))
+call read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
+call read_TS_finalize()
+do k=1,ndepth
+     ndata = 0; temp_t=0; temp_s=0
+     do j=1,nlat
+     do i=1,nlon
+       if(TEMP(i,j,k).gt.-10.0) then
+         ndata = ndata + 1
+         temp_t = temp_t + TEMP(i,j,k)
+         temp_s = temp_s + SALT(i,j,k)
+       endif
+     enddo
+     enddo
+     mTEMP(:,k) = temp_t / float(ndata)
+     mSALT(:,k) = temp_s / float(ndata)
+     write(6,*) ndata,mTemp(1,k),mSalt(1,k)
+enddo
+
+! generate initial conditions
+call get_init_conditions
+
+! 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
+write(6,*) ' calling write_OpenDX'
+write(6,*)
+call write_OpenDX(           on_a_sphere, &amp;
+                             nCellsNew, &amp;
+                             nVerticesNew, &amp;
+                             nEdgesNew, &amp;
+                             vertexDegreeNew, &amp;
+                             maxEdgesNew, &amp;
+                             xCellNew, &amp;
+                             yCellNew, &amp;
+                             zCellNew, &amp;
+                             xVertexNew, &amp;
+                             yVertexNew, &amp;
+                             zVertexNew, &amp;
+                             xEdgeNew, &amp;
+                             yEdgeNew, &amp;
+                             zEdgeNew, &amp;
+                             nEdgesOnCellNew, &amp;
+                             verticesOnCellNew, &amp;
+                             verticesOnEdgeNew, &amp;
+                             cellsOnVertexNew, &amp;
+                             edgesOnCellNew, &amp;
+                             areaCellNew, &amp;
+                             maxLevelCellNew, &amp;
+                             depthCell, &amp;
+                             kiteAreasOnVertexNew )
+
+
+!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
+
+
+!Step 4: Check the Initial conditions routine get_init_conditions
+subroutine get_init_conditions
+implicit none
+real :: halfwidth, dtr, pi, p(3), q(3), xin, yin, zin, ulon, ulat, stress, n1, n2, distance, r, temp_t, temp_s
+real :: dotProd
+integer :: iTracer, ix, iy, ndata, i, j, k, ixt, iyt, ncull, jcount, iNoData, kdata(nVertLevelsMod)
+logical :: flag_lat
+
+pi = 4.0*atan(1.0)
+dtr = pi/180.0
+
+!hNew = 10.0
+!salinityNew = 1.0
+
+hNew(1,:,:) = dble(thckNew(1,:,:))
+salinityNew(1,:,:) = dble(topgNew(1,:,:))
+
+temperatureNew = 1.0
+tracer1New = 1.0
+uNew = 0
+vNew = 0
+
+! ice sheet vars
+!thckNew(1,:,:) = dble(thckNew0(1,:,:))
+!topgNew(1,:,:) = dble(topgNew0(1,:,:))
+!deallocate( thckNew0 )
+!deallocate( topgNew0 )
+
+write(6,*) ' done get_init_conditions'
+
+end subroutine get_init_conditions
+
+
+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,&amp;
+                       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(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(fEdge(nEdges))
+allocate(fVertex(nVertices))
+allocate(h_s(nCells))
+allocate(work1(nCells))
+allocate(u_src(nVertLevels,nEdges))
+allocate(u(1,nVertLevels,nEdges))
+allocate(v(1,nVertLevels,nEdges))
+allocate(h(1,nVertLevels,nCells))
+allocate(rho(1,nVertLevels,nCells))
+
+xCell=0; yCell=0; zCell=0; latCell=0; lonCell=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
+
+fEdge=0; fVertex=0; h_s=0; u_src=0; work1=0
+u=0; v=0; h=0; rho=0
+
+
+call  read_netcdf_fields( &amp;
+                    time, &amp;
+                    latCell, &amp;
+                    lonCell, &amp;
+                    xCell, &amp;
+                    yCell, &amp;
+                    zCell, &amp;
+                    indexToCellID, &amp;
+                    latEdge, &amp;
+                    lonEdge, &amp;
+                    xEdge, &amp;
+                    yEdge, &amp;
+                    zEdge, &amp;
+                    indexToEdgeID, &amp;
+                    latVertex, &amp;
+                    lonVertex, &amp;
+                    xVertex, &amp;
+                    yVertex, &amp;
+                    zVertex, &amp;
+                    indexToVertexID, &amp;
+                    cellsOnEdge, &amp;
+                    nEdgesOnCell, &amp;
+                    nEdgesOnEdge, &amp;
+                    edgesOnCell, &amp;
+                    edgesOnEdge, &amp;
+                    weightsOnEdge, &amp;
+                    dvEdge, &amp;
+                    dcEdge, &amp;
+                    angleEdge, &amp;
+                    areaCell, &amp;
+                    areaTriangle, &amp;
+                    cellsOnCell, &amp;
+                    verticesOnCell, &amp;
+                    verticesOnEdge, &amp;
+                    edgesOnVertex, &amp;
+                    cellsOnVertex, &amp;
+                    kiteAreasOnVertex, &amp;
+                    fEdge, &amp;
+                    fVertex, &amp;
+                    h_s, &amp;
+                    u, &amp;
+                    v, &amp;
+                    h &amp;
+                   )
+
+write(6,*) ' values from read grid, min/max'
+write(6,*) ' latCell : ', minval(latCell), maxval(latCell)
+write(6,*) ' lonCell : ', minval(lonCell), maxval(lonCell)
+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,*) ' fEdge : ', minval(fEdge), maxval(fEdge)
+write(6,*) ' fVertex : ', minval(fVertex), maxval(fVertex)
+write(6,*) ' h_s : ', minval(h_s), maxval(h_s)
+write(6,*) ' u : ', minval(u), maxval(u)
+write(6,*) ' v : ', minval(v), maxval(v)
+write(6,*) ' h : ', minval(h), maxval(h)
+
+end subroutine read_grid
+
+
+subroutine write_grid
+implicit none
+
+if(on_a_sphere.eq.'YES              ') 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( &amp;
+                nCellsNew, &amp;
+                nEdgesNew, &amp;
+                nVerticesNew, &amp;
+                maxEdgesNew, &amp;
+                nVertLevelsNew, &amp;
+                vertexDegreeNew, &amp;
+                sphere_radius, &amp;
+                on_a_sphere &amp;
+                )
+
+call write_netcdf_fields( &amp;
+                    1, &amp;
+                    latCellNew, &amp;
+                    lonCellNew, &amp;
+                    xCellNew, &amp;
+                    yCellNew, &amp;
+                    zCellNew, &amp;
+                    indexToCellIDNew, &amp;
+                    latEdgeNew, &amp;
+                    lonEdgeNew, &amp;
+                    xEdgeNew, &amp;
+                    yEdgeNew, &amp;
+                    zEdgeNew, &amp;
+                    indexToEdgeIDNew, &amp;
+                    latVertexNew, &amp;
+                    lonVertexNew, &amp;
+                    xVertexNew, &amp;
+                    yVertexNew, &amp;
+                    zVertexNew, &amp;
+                    indexToVertexIDNew, &amp;
+                    maxLevelCellNew, &amp;
+                    cellsOnEdgeNew, &amp;
+                    nEdgesOnCellNew, &amp;
+                    nEdgesOnEdgeNew, &amp;
+                    edgesOnCellNew, &amp;
+                    edgesOnEdgeNew, &amp;
+                    weightsOnEdgeNew, &amp;
+                    dvEdgeNew, &amp;
+                    dcEdgeNew, &amp;
+                    angleEdgeNew, &amp;
+                    areaCellNew, &amp;
+                    areaTriangleNew, &amp;
+                    cellsOnCellNew, &amp;
+                    verticesOnCellNew, &amp;
+                    verticesOnEdgeNew, &amp;
+                    edgesOnVertexNew, &amp;
+                    cellsOnVertexNew, &amp;
+                    kiteAreasOnVertexNew, &amp;
+                    fEdgeNew, &amp;
+                    fVertexNew, &amp;
+                    h_sNew, &amp;
+                    boundaryEdgeNew, &amp;
+                    boundaryVertexNew, &amp;
+                    u_srcNew, &amp;
+                    uNew, &amp;
+                    vNew, &amp;
+                    hNew, &amp;
+                    rhoNew, &amp;
+                    temperatureNew, &amp;
+                    salinityNew, &amp;
+                    tracer1New, &amp;
+                    temperatureRestoreNew, &amp;
+                    salinityRestoreNew &amp;
+!                    thckNew, &amp;
+!                    topgNew &amp;
+                   )
+
+call write_netcdf_finalize
+
+if(on_a_sphere.eq.'YES              ') 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 5: Check the depth routine define_kmt
+subroutine define_kmt
+implicit none
+
+real (kind=4), allocatable, dimension(:) :: x,y, work_kmt
+
+integer :: nx, ny, inx, iny, ix, iy, nNewCells, counter
+
+real :: pi, dtr, rlon, rlat, r, ymin, ymax, xmin, xmax
+real :: latmin, latmax, lonmin, lonmax
+logical :: flag, kmt_flag
+
+pi = 4.0*atan(1.0)
+dtr = pi / 180.0
+
+allocate(kmt(nCells))
+kmt = 0
+allocate(thckdata(nCells))
+thckdata = 0
+allocate(topgdata(nCells))
+topgdata = 0
+
+!    nx = 301        ! 5 km  
+!    ny = 561
+    nx = 76         ! 20 km  
+    ny = 141
+    allocate(x(nx))
+    allocate(y(ny))
+    x = 0.0
+    y = 0.0
+
+    allocate(thck(nx,ny))
+    thck = 0.0
+    allocate(topg(nx,ny))
+    topg = 0.0
+
+    call read_topo_init( inx, iny)
+    write(6,*) ' inx, iny ', inx, iny                     
+    if(inx.ne.nx) stop ' nx topo'
+    if(iny.ne.ny) stop ' ny topo'
+
+    call read_topo_fields(x,y,thck,topg)
+!    write(6,*) 'thck = ', thck                         
+!    write(6,*) 'topg = ', topg                  
+
+    call read_topo_finalize( )
+
+    !!! debugging output !!!
+    write(6,*) ' '
+    print *, 'min(x), max(x)'
+    write(6,*) minval(x), maxval(x)
+    print *, 'min(y), max(y)'
+    write(6,*) minval(y), maxval(y)
+    print *, 'x(1), y(1)'
+    write(6,*) x(1), y(1)
+    write(6,*) ' '
+
+    print *, 'min(xCell), max(xCell)' 
+    write(6,*) minval(xCell), maxval(xCell)
+    print *, 'min(yCell), max(yCell)' 
+    write(6,*) minval(yCell), maxval(yCell)
+    print *, 'xCell(1), yCell(1)'
+    write(6,*) xCell(1), yCell(1)
+    write(6,*) ' '
+
+!    pause
+
+    do iCell=1,nCells
+
+        ix = nint( xCell(iCell) / 20.0d3 )
+        iy = nint( yCell(iCell) / 20.0d3 )
+!        ix = nint( xCell(iCell) / 5.0d3 )
+!        iy = nint( yCell(iCell) / 5.0d3 )
+
+        ix = min(nx,ix); iy=min(ny,iy)
+        ix = max(1,ix); iy=max(1,iy)
+
+        print *, '- new cell -'
+        print *, 'xCell(iCell)', xCell(iCell)
+        print *, 'yCell(iCell)', yCell(iCell)
+        print *, 'ix,iy = ', ix,iy 
+
+        thckdata(iCell) = thck(ix,iy)
+        topgdata(iCell) = topg(ix,iy)
+        print *, 'x1(ix,iy) = ', x(ix) 
+        print *, 'y1(ix,iy) = ', y(iy) 
+        print *, 'thck(ix,iy) = ', thck(ix,iy)     
+        print *, 'topg(ix,iy) = ', topg(ix,iy)   
+        print *, 'thckdata(iCell) = ', thckdata(iCell)
+        print *, 'topgdata(iCell) = ', topgdata(iCell)
+
+        ! create new mask for where geometry data exists
+         if( thckdata(iCell) .gt. 0 )then
+           kmt(iCell) = nVertLevelsMod
+           print *, 'cell kept'
+           print *, ' '
+         else
+           kmt(iCell) = 0
+           print *, 'cell removed', ix, iy
+           print *, ' '
+         endif
+
+    enddo
+
+    ! output no. of cells removed
+    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)
+
+    ! fill &quot;thckNew&quot; array and &quot;topgNew&quot; array
+    nNewCells = nCells - sum(work_kmt)
+    allocate(thckNew(1,1,nNewCells))
+    allocate(topgNew(1,1,nNewCells))
+    thckNew = 0
+    topgNew = 0
+!    allocate(thckNew0(1,1,nNewCells))
+!    allocate(topgNew0(1,1,nNewCells))
+!    thckNew0 = 0
+!    topgNew0 = 0
+    counter = 1
+    do iCell=1,nCells
+      if( kmt(iCell) /= 0 )then
+!        thckNew0(1,1,counter) = thckdata(iCell)   
+!        topgNew0(1,1,counter) = topgdata(iCell)   
+        thckNew(1,1,counter) = thckdata(iCell)   
+        topgNew(1,1,counter) = topgdata(iCell)   
+        counter = counter + 1
+      end if
+    end do
+
+    deallocate(work_kmt)
+    deallocate(x)
+    deallocate(y)
+    deallocate(thck)
+    deallocate(topg)
+
+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(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))
+allocate(maxLevelCellNew(nCellsNew))
+allocate(depthCell(nCellsNew))
+
+allocate(fEdgeNew(nEdgesNew))
+allocate(fVertexNew(nVerticesNew))
+allocate(h_sNew(nCellsNew))
+allocate(u_srcNew(nVertLevelsNew,nEdgesNew))
+allocate(uNew(1,nVertLevelsNew,nEdgesNew))
+allocate(vNew(1,nVertLevelsNew,nEdgesNew))
+allocate(hNew(1,nVertLevelsNew,nCellsNew))
+allocate(rhoNew(1,nVertLevelsNew,nCellsNew))
+allocate(temperatureNew(1,nVertLevelsNew,nCellsNew))
+allocate(salinityNew(1,nVertLevelsNew,nCellsNew))
+allocate(tracer1New(1,nVertLevelsNew,nCellsNew))
+
+! ice sheet variables
+!allocate(thckNew(1,nVertLevelsNew,nCellsNew))
+!allocate(topgNew(1,nVertLevelsNew,nCellsNew))
+
+allocate(temperatureRestoreNew(nCellsNew))
+allocate(salinityRestoreNew(nCellsNew))
+
+xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0
+xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
+xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
+
+fEdgeNew=0; fVertexNew=0; h_sNew=0; u_srcNew=0
+uNew=0; vNew=0; hNew=0; rhoNew=0
+temperatureNew=0; salinityNew=0; tracer1New=0;
+
+!thckNew = 0; topgNew = 0;
+
+temperatureRestoreNew = 0.0
+salinityRestoreNew = 0.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)
+    areaCellNew(jNew)=areaCell(i)
+    maxLevelCellNew(jNew) = kmt(i)
+    depthCell(jNew) = kmt(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)
+    fEdgeNew(jNew) = fEdge(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)
+    fVertexNew(jNew)=fVertex(i)
+    areaTriangleNew(jNew)=areaTriangle(i)
+endif
+enddo
+
+deallocate(xCell)
+deallocate(yCell)
+deallocate(zCell)
+deallocate(latCell)
+deallocate(lonCell)
+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 &quot;cellsOnEdge&quot;
+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) then
+    write(6,*) iEdge
+    write(6,*) iEdgeNew
+    write(6,*) iVertex1New, iVertex1
+    write(6,*) iVertex2New, iVertex2
+    stop &quot;verticesOnEdge&quot;
+endif
+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 &quot;jEdgeNew&quot;
+    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
+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 &quot;edgesOnCell&quot;
+edgesOnCellNew(i,iCellNew) = jNew
+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 &quot;verticesOnCell&quot;
+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
+
+
+subroutine get_dz
+integer k
+
+  dz( 1) = 10.0   !   5.006218       10.01244
+!  dz( 1) = 1001.244   !   5.006218       10.01244
+!  dz( 2) = 1011.258   !   15.06873       20.12502
+!  dz( 3) = 1031.682   !   25.28342       30.44183
+!  dz( 4) = 1063.330   !   35.75848       41.07513
+!  dz( 5) = 1107.512   !   46.61269       52.15025
+!  dz( 6) = 1166.145   !   57.98098       63.81171
+!  dz( 7) = 1241.928   !   70.02135       76.23099
+!  dz( 8) = 1338.612   !   82.92405       89.61711
+!  dz( 9) = 1461.401   !   96.92412       104.2311
+!  dz(10) = 1617.561   !   112.3189       120.4067
+!  dz(11) = 1817.368   !   129.4936       138.5804
+!  dz(12) = 2075.558   !   148.9582       159.3360
+!  dz(13) = 2413.680   !   171.4044       183.4728
+!  dz(14) = 2863.821   !   197.7919       212.1110
+!  dz(15) = 3474.644   !   229.4842       246.8575
+!  dz(16) = 4320.857   !   268.4617       290.0660
+!  dz(17) = 5516.812   !   317.6501       345.2342
+!  dz(18) = 7230.458   !   381.3865       417.5388
+!  dz(19) = 9674.901   !   465.9133       514.2878
+!  dz(20) = 13003.92   !   579.3074       644.3270
+!  dz(21) = 17004.89   !   729.3514       814.3759
+!  dz(22) = 20799.33   !   918.3725       1022.369
+!  dz(23) = 23356.94   !   1139.154       1255.939
+!  dz(24) = 24527.19   !   1378.574       1501.210
+!  dz(25) = 24898.04   !   1625.701       1750.191
+!  dz(26) = 24983.22   !   1875.107       2000.023
+!  dz(27) = 24997.87   !   2125.012       2250.002
+!  dz(28) = 24999.79   !   2375.000       2500.000
+!  dz(29) = 24999.98   !   2625.000       2749.999
+!  dz(30) = 25000.00   !   2874.999       2999.999
+!  dz(31) = 25000.00   !   3124.999       3249.999
+!  dz(32) = 25000.00   !   3374.999       3499.999
+!  dz(33) = 25000.00   !   3624.999       3749.999
+!  dz(34) = 25000.00   !   3874.999       3999.999
+!  dz(35) = 25000.00   !   4124.999       4249.999
+!  dz(36) = 25000.00   !   4374.999       4499.999
+!  dz(37) = 25000.00   !   4624.999       4749.999
+!  dz(38) = 25000.00   !   4874.999       4999.999
+!  dz(39) = 25000.00   !   5124.999       5249.999
+!  dz(40) = 25000.00   !   5374.999       5499.999
+
+!  dz = dz / 100.0
+
+  write(6,*)
+!  do k=1,40
+  do k=1,1
+    print *, 'k, dz = ' 
+    write(6,*) k,dz(k)
+  enddo
+  write(6,*)
+
+end subroutine get_dz
+end program map_to_icesheet

Modified: branches/land_ice/icesheet/src/module_read_topo.F
===================================================================
--- branches/land_ice/icesheet/src/module_read_topo.F        2011-08-18 23:07:37 UTC (rev 949)
+++ branches/land_ice/icesheet/src/module_read_topo.F        2011-08-19 14:55:41 UTC (rev 950)
@@ -3,10 +3,11 @@
    integer :: rd_ncid
    integer :: rdDimIDnx
    integer :: rdDimIDny
-!   integer :: rdVarIDz
    integer :: rdVarIDx
    integer :: rdVarIDy
    integer :: rdVarIDthk
+   integer :: rdVarIDtopg
+
  
    integer :: rdLocalnx
    integer :: rdLocalny
@@ -54,26 +55,33 @@
       !
       nferr = nf_inq_varid(rd_ncid, 'x1', rdVarIDx)
       write(6,*) ' nferr xvar-id ', nferr, rdVarIDx
+
       nferr = nf_inq_varid(rd_ncid, 'y1', rdVarIDy)
       write(6,*) ' nferr yvar-id ', nferr, rdVarIDy
+
       nferr = nf_inq_varid(rd_ncid, 'thk', rdVarIDthk)
       write(6,*) ' nferr thkvar-id ', nferr, rdVarIDthk
 
+      nferr = nf_inq_varid(rd_ncid, 'topg', rdVarIDtopg)
+      write(6,*) ' nferr topgvar-id ', nferr, rdVarIDtopg
+
       nferr = nf_inq_vardimid( rd_ncid, rdVarIDthk, dimids )
       write(6,*) ' nferr thk-dim-ids ', nferr, dimids  
 
+      nferr = nf_inq_vardimid( rd_ncid, rdVarIDtopg, dimids )
+      write(6,*) ' nferr topg-dim-ids ', nferr, dimids  
  
    end subroutine read_topo_init
  
  
-   subroutine read_topo_fields(x,y,thk)
+   subroutine read_topo_fields(x,y,thk,topg)
  
       implicit none
  
       include 'netcdf.inc'
  
       real (kind=4), dimension(:), intent(out) :: x,y
-      real (kind=4), dimension(:,:), intent(out) :: thk
+      real (kind=4), dimension(:,:), intent(out) :: thk, topg
 
       integer, dimension(1) :: start1, count1
       integer, dimension(2) :: start2, count2
@@ -107,6 +115,9 @@
       nferr = nf_get_vara_real(rd_ncid, rdVarIDthk, start3, count3, thk)
       write(6,*) ' nferr ncid thkvar-id (get thk val) ', nferr, rd_ncid, rdVarIDthk
  
+      nferr = nf_get_vara_real(rd_ncid, rdVarIDtopg, start3, count3, topg)
+      write(6,*) ' nferr ncid topgvar-id (get topg val) ', nferr, rd_ncid, rdVarIDtopg
+
    end subroutine read_topo_fields
  
  

Modified: branches/land_ice/icesheet/src/module_write_netcdf.F
===================================================================
--- branches/land_ice/icesheet/src/module_write_netcdf.F        2011-08-18 23:07:37 UTC (rev 949)
+++ branches/land_ice/icesheet/src/module_write_netcdf.F        2011-08-19 14:55:41 UTC (rev 950)
@@ -61,7 +61,11 @@
    integer :: wrVarIDtracer1
    integer :: wrVarIDtemperatureRestore
    integer :: wrVarIDsalinityRestore
- 
+
+   ! ice sheet variables 
+!   integer :: wrVarIDthck    
+!   integer :: wrVarIDtopg    
+
    integer :: wrLocalnCells
    integer :: wrLocalnEdges
    integer :: wrLocalnVertices
@@ -256,6 +260,16 @@
       ! If you do not want tracer1 in your input file, simply comment out these two lines (one of two)
       nferr = nf_def_var(wr_ncid, 'tracer1', NF_DOUBLE,  3, dimlist, wrVarIDtracer1)
 
+      ! ice sheet variables (NOTE: for now, we are assuming vert dim of 1, so this works). If nVertLevels &gt; 1, this dimensioning
+      ! won't work, since thickness and topography are only 2d fields)
+!      nferr = nf_def_var(wr_ncid, 'thck', NF_DOUBLE,  3, dimlist, wrVarIDthck)
+!      dimlist( 1) = wrDimIDnVertLevels
+!      dimlist( 2) = wrDimIDnCells
+!      dimlist( 3) = wrDimIDTime
+!      nferr = nf_def_var(wr_ncid, 'topg', NF_DOUBLE,  3, dimlist, wrVarIDtopg)
+!      dimlist( 1) = wrDimIDnVertLevels
+!      dimlist( 2) = wrDimIDnCells
+!      dimlist( 3) = wrDimIDTime
  
       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)
@@ -318,7 +332,9 @@
                                   tracer1, &amp;
                                   temperatureRestore, &amp;
                                   salinityRestore &amp;
-                                 )
+!                                  thck, &amp;
+!                                  topg &amp;
+                                  )
  
       implicit none
  
@@ -377,6 +393,9 @@
       real (kind=8), dimension(:), intent(in) :: temperatureRestore
       real (kind=8), dimension(:), intent(in) :: salinityRestore
 
+      ! ice sheet fields     
+!      real (kind=8), dimension(:,:,:), intent(in) :: thck        
+!      real (kind=8), dimension(:,:,:), intent(in) :: topg    
  
       integer :: nferr
       integer, dimension(1) :: start1, count1
@@ -629,7 +648,20 @@
       count3( 3) = 1
       ! If you do not want tracer1 in your input file, simply comment out these two lines (two of two)
       nferr = nf_put_vara_double(wr_ncid, wrVarIDtracer1, start3, count3, tracer1)
- 
+
+      ! ice sheet fields 
+!      start3(3) = time
+!      count3( 1) = wrLocalnVertLevels
+!      count3( 2) = wrLocalnCells
+!      count3( 3) = 1
+!      nferr = nf_put_vara_double(wr_ncid, wrVarIDthck, start3, count3, thck)
+!
+!      start3(3) = time
+!      count3( 1) = wrLocalnVertLevels
+!      count3( 2) = wrLocalnCells
+!      count3( 3) = 1
+!      nferr = nf_put_vara_double(wr_ncid, wrVarIDtopg, start3, count3, topg)
+
    end subroutine write_netcdf_fields
  
  

</font>
</pre>