<p><b>sprice@lanl.gov</b> 2012-01-09 15:19:16 -0700 (Mon, 09 Jan 2012)</p><p><br>
Branch Commit (land ice):<br>
<br>
Add some directories to the land ice branch.<br>
</p><hr noshade><pre><font color="gray">Copied: branches/land_ice/icesheet/README (from rev 1322, branches/land_ice/icesheet/README)
===================================================================
--- branches/land_ice/icesheet/README (rev 0)
+++ branches/land_ice/icesheet/README 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,43 @@
+
+This dir contains a hacked version of the "basin" code from the "ocean_projects" dir.
+The goal is to do something similar but for ice sheet grids. That is ...
+
+1) use periodic_hex code to generate a periodic hex grid w/ appropriate overall dims that contain
+the geometry data from the ice sheet of interest (e.g. that encompass Greenland or Ant.)
+
+2) use altered version of basin code (icesheet.F) to extract cells from that grid for which there are
+non-zero geom data (e.g. for which thickness data are non-zero)
+
+The .nc files containing the original data for ice sheets are in "topo". The "grid.nc" file for the
+generic, MPAS periodic hex mesh are in the "periodic_hex" subdir. See additional notes in the README
+file there for information on creating that (necessary) "grid.nc" file.
+
+To execute the code, copy or link to the "grid.nc" file in the "periodic_hex" subdir, then execute
+the script "runit". This will build and run the code that does all of the work (reads in the topo
+data and information from the original MPAS hex mesh, compares them, and extracts only cells from the
+MPAS mesh for which data exist).
+
+sfp 8/16/11
+
+
+
+--- original contents of this readme file ---
+
+Program: basin
+
+This program reads in a MPAS grid file (grid.nc) and outputs a modified grid file
+(ocean.nc) and it's associated graph.info files for partitions from 2 up to 1024 in powers of 2.
+
+The purpose of this code is to remove grid cells from any valid MPAS grid.
+Please see source file src/basin.F to define the specifics of the output grid file.
+
+After a grid.nc file has been placed in this directory, simply run the script runit.
+
+This script will compile basin, run basin (producing an ocean.nc file) and use metis
+to partition the graph.info file.
+
+The metis portion of the script requires the executable kmetis to be in your path.
+If it is not, none of the graph.info.part.N files will be created, but can later be
+created manually using metis and the assocaite graph.info file.
+
+
Copied: branches/land_ice/icesheet/periodic_hex/README (from rev 1322, branches/land_ice/icesheet/periodic_hex/README)
===================================================================
--- branches/land_ice/icesheet/periodic_hex/README (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/README 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,24 @@
+
+This dir contains files for use in creating MPAS period hex mesh grids based on ...
+
+1) a generic perodic hex mesh constructed using the "namelist.input" files stored here
+
+2) a given input .nc dataset, stored in "../topo/"
+
+3) the "periodic_grid" executable in "../../grid_gen/periodic_hex/"
+
+To construct the mesh, first link to the necessary executable ...
+
+ln -s ../../grid_gen/periodic_hex/periodic_grid
+
+Then copy one of the "namelist.specific-file-info.input" files to "namelist.input".
+
+Next, execute "periodic_grid", which will produce a generic MPAS periodic hex mesh
+with the dimensions specified in the "namelist.input" file. This will be used along
+with the other code in "../src/" and the .nc files in "../topo/" to build a more
+specific, boundary conforming hex mesh, based on the geometry read from the input
+.nc files in "../topo/".
+
+Additional instructions for these later steps can be found one dir level up.
+
+sfp 8/16/11
Copied: branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.20km (from rev 1322, branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.20km)
===================================================================
--- branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.20km (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.20km 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,8 @@
+&periodic_grid
+ nx = 76,
+ ny = 163,
+ dc = 20000.,
+ nVertLevels = 1,
+ nTracers = 1,
+ nproc = 2,
+/
Copied: branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.5km (from rev 1322, branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.5km)
===================================================================
--- branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.5km (rev 0)
+++ branches/land_ice/icesheet/periodic_hex/namelist.input.Greenland.5km 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,8 @@
+&periodic_grid
+ nx = 301,
+ ny = 648,
+ dc = 5000.,
+ nVertLevels = 3,
+ nTracers = 1,
+ nproc = 2,
+/
Copied: branches/land_ice/icesheet/runit (from rev 1322, branches/land_ice/icesheet/runit)
===================================================================
--- branches/land_ice/icesheet/runit (rev 0)
+++ branches/land_ice/icesheet/runit 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,20 @@
+rm -f map
+cd src
+make clean
+make
+mv map ..
+cd ..
+rm -f graph*
+./map
+#cd ../metis-4.0
+kmetis graph.info 2
+kmetis graph.info 4
+kmetis graph.info 8
+kmetis graph.info 16
+kmetis graph.info 32
+kmetis graph.info 64
+kmetis graph.info 128
+kmetis graph.info 256
+kmetis graph.info 512
+kmetis graph.info 1024
+#cd ../basin
Copied: branches/land_ice/icesheet/src/Makefile (from rev 1322, branches/land_ice/icesheet/src/Makefile)
===================================================================
--- branches/land_ice/icesheet/src/Makefile (rev 0)
+++ branches/land_ice/icesheet/src/Makefile 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,71 @@
+# IBM with Xlf compilers
+#FC = xlf90
+#CC = xlc
+#FFLAGS = -qrealsize=8 -g -C
+#CFLAGS = -g
+#LDFLAGS = -g -C
+
+# pgf90
+#FC = pgf90
+#CC = pgcc
+#FFLAGS = -r8 -O3
+#CFLAGS = -O3
+#LDFLAGS = -O3
+
+# ifort
+#FC = ifort
+#CC = icc
+#FFLAGS = -real-size 64 -g -traceback -check all
+#CFLAGS = -g
+#LDFLAGS = -g -traceback -check all
+
+# gfortran
+FC = gfortran
+CC = gcc
+FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian
+CFLAGS = -O3 -m64
+LDFLAGS = -O3 -m64
+
+# absoft
+#FC = f90
+#CC = gcc
+#FFLAGS = -dp -O3
+#CFLAGS = -O3
+#LDFLAGS = -O3
+#NETCDF = /Users/maltrud/local
+
+
+CPP = cpp -C -P -traditional
+CPPFLAGS =
+CPPINCLUDES =
+INCLUDES = -I$(NETCDF)/include
+LIBS = -L$(NETCDF)/lib -lnetcdf
+
+RM = rm -f
+
+##########################
+
+.SUFFIXES: .F .o
+
+
+OBJS = icesheet.o \
+ utilities.o \
+ module_read_netcdf.o \
+ module_read_topo.o \
+ module_write_netcdf.o
+
+all: map
+
+icesheet.o: utilities.o module_write_netcdf.o module_read_netcdf.o module_read_topo.o
+
+map: $(OBJS)
+ $(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
+
+clean:
+ $(RM) *.o *.mod pop *.f90
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(INCLUDES)
+ #$(RM) $*.f90
Copied: branches/land_ice/icesheet/src/icesheet.F (from rev 1322, branches/land_ice/icesheet/src/icesheet.F)
===================================================================
--- branches/land_ice/icesheet/src/icesheet.F (rev 0)
+++ branches/land_ice/icesheet/src/icesheet.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,1259 @@
+program map_to_icesheet
+
+use read_netcdf
+use read_topo
+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, work1
+
+! ice sheet variables
+real(kind=4), allocatable, dimension(:,:) :: thck, topg
+
+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 "on_a_sphere" and "sphere_radius". 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 :: &
+ h_total_max = 2000.0, &
+ u_max = 0.0, &
+ u_src_max = 0.1, & ! max wind stress, N/m2
+ beta = 1.4e-11, &
+ f0 = -1.1e-4, &
+ 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
+
+real, allocatable, dimension(:,:,:) :: hNew
+
+!ice sheet variables
+real(kind=8), allocatable, dimension(:,:,:) :: thckNew, topgNew
+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
+
+! 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, &
+ nCellsNew, &
+ nVerticesNew, &
+ nEdgesNew, &
+ vertexDegreeNew, &
+ maxEdgesNew, &
+ xCellNew, &
+ yCellNew, &
+ zCellNew, &
+ xVertexNew, &
+ yVertexNew, &
+ zVertexNew, &
+ xEdgeNew, &
+ yEdgeNew, &
+ zEdgeNew, &
+ nEdgesOnCellNew, &
+ verticesOnCellNew, &
+ verticesOnEdgeNew, &
+ cellsOnVertexNew, &
+ edgesOnCellNew, &
+ areaCellNew, &
+ maxLevelCellNew, &
+ depthCell, &
+ 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
+
+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)
+ 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)
+ t = q - p
+ s = s - r
+ call unit_vector_in_3space(t)
+ call unit_vector_in_3space(s)
+ dot = s(1)*t(1)+s(2)*t(2)+s(3)*t(3)
+ if(dot.lt.mindot) mindot=dot
+ if(dot.gt.maxdot) maxdot=dot
+enddo
+write(6,10) 'orthogonality of cell and vertex edges (should be zeros)', mindot, maxdot
+
+! check that the kiteareas sum to the areatriangle
+mindot = 2
+maxdot = -2
+do iVertex=1,nVerticesNew
+ b = 0
+ do i=1,vertexDegree
+ b(i) = kiteAreasOnVertexNew(i,iVertex)
+ enddo
+ angle = sum(b)
+ if(angle - areaTriangleNew(iVertex).lt.mindot) mindot = angle - areaTriangleNew(iVertex)
+ if(angle - areaTriangleNew(iVertex).gt.maxdot) maxdot = angle - areaTriangleNew(iVertex)
+enddo
+write(6,10) ' error in sum of kites and triangles (should be zeroes)', mindot, maxdot
+
+! check that the kiteareas sum to the areaCell
+mindot = 2
+maxdot = -2
+work = 0
+do iVertex=1,nVerticesNew
+ iCell1 = cellsOnVertexNew(1,iVertex)
+ iCell2 = cellsOnVertexNew(2,iVertex)
+ iCell3 = cellsOnVertexNew(3,iVertex)
+ if(iCell1.ne.0) work(iCell1) = work(iCell1) + kiteAreasOnVertexNew(1,iVertex)
+ if(iCell2.ne.0) work(iCell2) = work(iCell2) + kiteAreasOnVertexNew(2,iVertex)
+ if(iCell3.ne.0) work(iCell3) = work(iCell3) + kiteAreasOnVertexNew(3,iVertex)
+enddo
+mindot = minval(areaCellNew - work)
+maxdot = maxval(areaCellNew - work)
+write(6,10) ' error in sum of kites and cells (should be zeroes)', mindot, maxdot
+
+!check for connectivity inverses for cells/edges
+do iCell=1,nCellsNew
+ do i=1,nEdgesOnCellNew(iCell)
+ iEdge=edgesOnCellNew(i,iCell)
+ if(iEdge.le.0) stop ' iEdge le 0'
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ if(iCell1.ne.iCell.and.iCell2.ne.iCell) stop ' cells/edges inverse failed'
+ enddo
+enddo
+write(6,*) ' cellsOnEdge and edgesOnCell are duals for every cell/edge combination'
+
+!check for connectivity inverses for cells/vertices
+do iCell=1,nCellsNew
+ do i=1,nEdgesOnCellNew(iCell)
+ iVertex = verticesOnCellNew(i,iCell)
+ if(iVertex.le.0) stop ' iVertex le 0'
+ iCell1 = cellsOnVertexNew(1,iVertex)
+ iCell2 = cellsOnVertexNew(2,iVertex)
+ iCell3 = cellsOnVertexNew(3,iVertex)
+ if(iCell1.ne.iCell.and.iCell2.ne.iCell.and.iCell3.ne.iCell) stop ' cells/vertices inverse failed'
+ enddo
+enddo
+write(6,*) ' cellsOnVertex and verticesOnCell are duals for every cell/vertex combination'
+
+!check edgesOnEdge
+do iEdge=1,nEdgesNew
+ iCell1 = cellsOnEdgeNew(1,iEdge)
+ iCell2 = cellsOnEdgeNew(2,iEdge)
+ if(nEdgesOnEdgeNew(iEdge).eq.0) then
+ if(boundaryEdgeNew(1,iEdge).ne.1) stop ' stopping boundaryEdgeNew'
+ endif
+ do i=1,nEdgesOnEdgeNew(iEdge)
+ jEdge = edgesOnEdgeNew(i,iEdge)
+ jCell1 = cellsOnEdgeNew(1,jEdge)
+ jCell2 = cellsOnEdgeNew(2,jEdge)
+ if(jCell1.ne.iCell1.and.jCell1.ne.iCell2) then
+ if(jCell2.ne.iCell1.and.jCell2.ne.iCell2) then
+ write(6,*) 'error in edgesOnEdge'
+ write(6,*) iCell1, iCell2, jCell1, jCell2
+ stop
+ endif
+ endif
+ enddo
+enddo
+write(6,*) ' edgesOnEdge is consistent with cellsOnEdge'
+
+end subroutine error_checking
+
+
+subroutine copy_dimensions
+
+maxEdgesNew = maxEdges
+maxEdges2New = maxEdges2
+TWONew = TWO
+vertexDegreeNew = vertexDegree
+nVertLevelsNew = nVertLevelsMod
+
+write(6,*)
+write(6,*) ' new dimensions '
+write(6,*) ' maxEdgesNew : ', maxEdgesNew
+write(6,*) ' maxEdges2New : ', maxEdges2New
+write(6,*) ' TWONew : ', TWONew
+write(6,*) ' vertexDegreeNew : ', vertexDegreeNew
+write(6,*) ' nVertLevelsNew : ', nVertLevelsNew
+
+end subroutine copy_dimensions
+
+
+
+subroutine read_grid
+implicit none
+
+call read_netcdf_init(nCells, nEdges, nVertices, maxEdges,maxEdges2,&
+ nVertLevels,TWO,vertexDegree)
+
+write(6,*) ' init from grid '
+write(6,*) 'nCells :', nCells
+write(6,*) 'nEdges :', nEdges
+write(6,*) 'nVertices :', nVertices
+write(6,*) 'maxEdges :', maxEdges
+write(6,*) 'maxEdges2 :', maxEdges2
+write(6,*) 'nVertLevels :', nVertLevels
+write(6,*) 'vertexDegree :', vertexDegree
+write(6,*) 'TWO :', TWO
+
+allocate(xCell(nCells))
+allocate(yCell(nCells))
+allocate(zCell(nCells))
+allocate(latCell(nCells))
+allocate(lonCell(nCells))
+allocate(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(work1(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
+
+
+call read_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ xCell, &
+ yCell, &
+ zCell, &
+ indexToCellID, &
+ latEdge, &
+ lonEdge, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ indexToEdgeID, &
+ latVertex, &
+ lonVertex, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex &
+ )
+
+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)
+
+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( &
+ nCellsNew, &
+ nEdgesNew, &
+ nVerticesNew, &
+ maxEdgesNew, &
+ nVertLevelsNew, &
+ vertexDegreeNew, &
+ sphere_radius, &
+ on_a_sphere &
+ )
+
+call write_netcdf_fields( &
+ 1, &
+ latCellNew, &
+ lonCellNew, &
+ xCellNew, &
+ yCellNew, &
+ zCellNew, &
+ indexToCellIDNew, &
+ latEdgeNew, &
+ lonEdgeNew, &
+ xEdgeNew, &
+ yEdgeNew, &
+ zEdgeNew, &
+ indexToEdgeIDNew, &
+ latVertexNew, &
+ lonVertexNew, &
+ xVertexNew, &
+ yVertexNew, &
+ zVertexNew, &
+ indexToVertexIDNew, &
+ maxLevelCellNew, &
+ cellsOnEdgeNew, &
+ nEdgesOnCellNew, &
+ nEdgesOnEdgeNew, &
+ edgesOnCellNew, &
+ edgesOnEdgeNew, &
+ weightsOnEdgeNew, &
+ dvEdgeNew, &
+ dcEdgeNew, &
+ angleEdgeNew, &
+ areaCellNew, &
+ areaTriangleNew, &
+ cellsOnCellNew, &
+ verticesOnCellNew, &
+ verticesOnEdgeNew, &
+ edgesOnVertexNew, &
+ cellsOnVertexNew, &
+ kiteAreasOnVertexNew, &
+ fEdgeNew, &
+ fVertexNew, &
+ boundaryEdgeNew, &
+ boundaryVertexNew, &
+ hNew, &
+ topgNew, &
+ thckNew &
+ )
+
+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 "thckNew" array and "topgNew" array
+ nNewCells = nCells - sum(work_kmt)
+ allocate(thckNew(1,1,nNewCells))
+ allocate(topgNew(1,1,nNewCells))
+ thckNew = 0
+ topgNew = 0
+ counter = 1
+ do iCell=1,nCells
+ if( kmt(iCell) /= 0 )then
+ 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(hNew(1,nVertLevelsNew,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
+
+hNew=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 "cellsOnEdge"
+if(iCell1New.eq.0) then
+ cellsOnEdgeNew(1,iEdgeNew) = iCell2New
+ cellsOnEdgeNew(2,iEdgeNew) = iCell1New
+ flipVerticesOnEdgeOrdering(iEdgeNew) = 1
+endif
+enddo
+deallocate(cellsOnEdge)
+
+allocate(verticesOnEdgeNew(TWONew,nEdgesNew))
+allocate(boundaryVertexNew(nVertLevelsNew,nVerticesNew))
+verticesOnEdgeNew(:,:) = 0
+boundaryVertexNew(:,:) = 0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+iVertex1 = VerticesOnEdge(1,iEdge)
+iVertex2 = VerticesOnEdge(2,iEdge)
+iVertex1New = vertexMap(iVertex1)
+iVertex2New = vertexMap(iVertex2)
+if(iVertex1New.eq.0.or.iVertex2New.eq.0) then
+ write(6,*) iEdge
+ write(6,*) iEdgeNew
+ write(6,*) iVertex1New, iVertex1
+ write(6,*) iVertex2New, iVertex2
+ stop "verticesOnEdge"
+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 "jEdgeNew"
+ edgesOnEdgeNew(i,iEdgeNew)=jEdgeNew
+ weightsOnEdgeNew(i,iEdgeNew) = weightsOnEdge(i,iEdge)
+ enddo
+endif
+enddo
+deallocate(nEdgesOnEdge)
+deallocate(edgesOnEdge)
+deallocate(weightsOnEdge)
+
+allocate(cellsOnCellNew(maxEdges,nCellsNew))
+allocate(nEdgesOnCellNew(nCellsNew))
+cellsOnCellNew = 0
+nEdgesOnCellNew = 0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+nEdgesOnCellNew(iCellNew)=nEdgesOnCell(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = cellsOnCell(i,iCell)
+jNew = cellMap(j)
+cellsOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
+deallocate(cellsOnCell)
+deallocate(nEdgesOnCell)
+
+allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
+edgesOnCellNew(:,:) = 0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = edgesOnCell(i,iCell)
+jNew = edgeMap(j)
+if(jNew.eq.0) stop "edgesOnCell"
+edgesOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
+deallocate(edgesOnCell)
+
+allocate(verticesOnCellNew(maxEdgesNew,nCellsNew))
+verticesOnCellNew(:,:)=0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j=verticesOnCell(i,iCell)
+jNew = vertexMap(j)
+if(jNew.eq.0) stop "verticesOnCell"
+verticesOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
+deallocate(verticesOnCell)
+
+allocate(cellsOnVertexNew(vertexDegreeNew,nVerticesNew))
+allocate(kiteAreasOnVertexNew(vertexDegreeNew,nVerticesNew))
+cellsOnVertexNew = 0
+kiteAreasOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=cellsOnVertex(i,iVertex)
+jNew=cellMap(j)
+if(jNew.eq.0) then
+ kiteAreasOnVertexNew(i,iVertexNew)=0
+else
+ kiteAreasOnVertexNew(i,iVertexNew)=kiteAreasOnVertex(i,iVertex)
+endif
+cellsOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
+deallocate(cellsOnVertex)
+deallocate(kiteAreasOnVertex)
+
+areaTriangleNew = 0
+do iVertex=1,nVerticesNew
+do i=1,vertexDegree
+areaTriangleNew(iVertex) = areaTriangleNew(iVertex) + kiteAreasOnVertexNew(i,iVertex)
+enddo
+enddo
+
+allocate(edgesOnVertexNew(vertexDegreeNew, nVerticesNew))
+edgesOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=edgesOnVertex(i,iVertex)
+jNew=edgeMap(j)
+edgesOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
+deallocate(edgesOnVertex)
+
+! find normals
+normalsNew = 0.0
+do iEdge=1,nEdgesNew
+cell1 = cellsOnEdgeNew(1,iEdge)
+cell2 = cellsOnEdgeNew(2,iEdge)
+if(cell1.eq.0.or.cell2.eq.0) cycle
+c1(1) = xCellNew(cell1); c1(2) = yCellNew(cell1); c1(3) = zCellNew(cell1)
+c2(1) = xCellNew(cell2); c2(2) = yCellNew(cell2); c2(3) = zCellNew(cell2)
+distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+
+if(on_a_sphere.eq.'YES ') then
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+else
+ if(distance.gt.0.5*Lx) then
+ write(6,*) ' periodic edge ', iEdge, distance
+ write(6,10) ' c1 ', c1(:)
+ write(6,10) ' c2 ', c2(:)
+ r = c2(1) - c1(1)
+ if(r.gt.0) c2(1) = c2(1) - Lx
+ if(r.lt.0) c2(1) = c2(1) + Lx
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ write(6,*) ' periodic edge fix ', iEdge, r, distance
+ endif
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+endif
+enddo
+10 format(a20,3e15.5)
+
+end subroutine map_connectivity
+
+
+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
Copied: branches/land_ice/icesheet/src/module_read_TS.F (from rev 1322, branches/land_ice/icesheet/src/module_read_TS.F)
===================================================================
--- branches/land_ice/icesheet/src/module_read_TS.F (rev 0)
+++ branches/land_ice/icesheet/src/module_read_TS.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,165 @@
+module read_TS
+
+ integer :: rd_ncid, rd_ncids, rd_ncidu
+ integer :: rdDimIDt_lon
+ integer :: rdDimIDt_lat
+ integer :: rdDimIDdepth_t
+ integer :: rdVarIDt_lon
+ integer :: rdVarIDt_lat
+ integer :: rdVarIDdepth_t
+ integer :: rdVarIDTEMP
+ integer :: rdVarIDSALT
+ integer :: rdVarIDTAUX
+ integer :: rdVarIDTAUY
+
+ integer :: rdLocalt_lon
+ integer :: rdLocalt_lat
+ integer :: rdLocaldepth_t
+
+ contains
+
+ subroutine read_TS_init(nx, ny, nz)
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(out) :: nx, ny, nz
+
+ integer :: nferr, nferrs, nferru
+
+ nferr = nf_open('TS/woce_t_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncid)
+ write(6,*) ' nferr ', nferr, rd_ncid
+
+ !
+ ! Get IDs for variable dimensions
+ !
+ nferr = nf_inq_dimid(rd_ncid, 't_lon', rdDimIDt_lon)
+ write(6,*) ' nferr ', nferr, rdDimIDt_lon
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lon, rdLocalt_lon)
+ write(6,*) ' nferr ', nferr, rdLocalt_lon
+ nferr = nf_inq_dimid(rd_ncid, 't_lat', rdDimIDt_lat)
+ write(6,*) ' nferr ', nferr, rdDimIDt_lat
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDt_lat, rdLocalt_lat)
+ write(6,*) ' nferr ', nferr, rdLocalt_lat
+ nferr = nf_inq_dimid(rd_ncid, 'depth_t', rdDimIDdepth_t)
+ write(6,*) ' nferr ', nferr, rdDimIDdepth_t
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDdepth_t, rdLocaldepth_t)
+ write(6,*) ' nferr ', nferr, rdLocaldepth_t
+
+ nx = rdLocalt_lon
+ ny = rdLocalt_lat
+ nz = rdLocaldepth_t
+
+ write(6,*) nx, ny, nz
+
+ !
+ ! Get IDs for variables
+ !
+ nferr = nf_inq_varid(rd_ncid, 't_lon', rdVarIDt_lon)
+ write(6,*) ' nferr ', nferr, rdVarIDt_lon
+ nferr = nf_inq_varid(rd_ncid, 't_lat', rdVarIDt_lat)
+ write(6,*) ' nferr ', nferr, rdVarIDt_lat
+ nferr = nf_inq_varid(rd_ncid, 'depth_t', rdVarIDdepth_t)
+ write(6,*) ' nferr ', nferr, rdVarIDdepth_t
+ nferr = nf_inq_varid(rd_ncid, 'TEMP', rdVarIDTEMP)
+ write(6,*) ' nferr ', nferr, rdVarIDTEMP
+
+ nferrs = nf_open('TS/woce_s_ann.3600x2431x42interp.r4.nc', NF_SHARE, rd_ncids)
+ nferrs = nf_inq_varid(rd_ncids, 'SALT', rdVarIDSALT)
+ write(6,*) ' nferrs ', nferrs, rdVarIDSALT
+
+ nferru = nf_open('TS/ws.old_ncep_1958-2000avg.interp3600x2431.nc', NF_SHARE, rd_ncidu)
+ nferru = nf_inq_varid(rd_ncidu, 'TAUX', rdVarIDTAUX)
+ nferru = nf_inq_varid(rd_ncidu, 'TAUY', rdVarIDTAUY)
+ write(6,*) ' nferru ', nferru, rdVarIDTAUX, rdVarIDTAUY
+
+ end subroutine read_TS_init
+
+ subroutine read_TS_fields(t_lon, t_lat, depth_t, TEMP, SALT, TAUX, TAUY)
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ real (kind=4), dimension(:), intent(out) :: t_lon, t_lat, depth_t
+ real (kind=4), dimension(:,:,:), intent(out) :: TEMP, SALT
+ real (kind=4), dimension(:,:), intent(out) :: TAUX, TAUY
+
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ integer :: nferr, nferrs, nferru
+
+ start1(1) = 1
+ count1(1) = rdLocalt_lon
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lon, start1, count1, t_lon)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lon
+
+ start1(1) = 1
+ count1(1) = rdLocalt_lat
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDt_lat, start1, count1, t_lat)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDt_lat
+
+ start1(1) = 1
+ count1(1) = rdLocaldepth_t
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDdepth_t, start1, count1, depth_t)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDdepth_t
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdLocalt_lon
+ count3(2) = rdLocalt_lat
+ count3(3) = rdLocaldepth_t
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDTEMP, start3, count3, TEMP)
+ write(6,*) ' nferr ', nferr, rd_ncid, rdVarIDTEMP
+ write(6,*) ' temperature' , minval(TEMP), maxval(TEMP)
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdLocalt_lon
+ count3(2) = rdLocalt_lat
+ count3(3) = rdLocaldepth_t
+ nferrs = nf_get_vara_real(rd_ncids, rdVarIDSALT, start3, count3, SALT)
+ write(6,*) ' nferrs ', nferrs, rd_ncids, rdVarIDSALT
+ write(6,*) ' salinity' , minval(SALT), maxval(SALT)
+
+ start2(1) = 1
+ start2(2) = 1
+ count2(1) = rdLocalt_lon
+ count2(2) = rdLocalt_lat
+ nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUX, start2, count2, TAUX)
+ nferru = nf_get_vara_real(rd_ncidu, rdVarIDTAUY, start2, count2, TAUY)
+ write(6,*) ' nferru ', nferru, rd_ncidu, rdVarIDTAUX, rdVarIDTAUY
+ write(6,*) ' TAUX' , minval(TAUX), maxval(TAUX)
+ write(6,*) ' TAUY' , minval(TAUY), maxval(TAUY)
+
+
+ end subroutine read_TS_fields
+
+
+ subroutine read_TS_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr, nferrs, nferru
+
+ nferr = nf_close(rd_ncid)
+ write(6,*) ' nferr ', nferr
+
+ nferrs = nf_close(rd_ncids)
+ write(6,*) ' nferrs ', nferrs
+
+ nferru = nf_close(rd_ncidu)
+ write(6,*) ' nferru ', nferru
+
+
+ end subroutine read_TS_finalize
+
+end module read_TS
Copied: branches/land_ice/icesheet/src/module_read_netcdf.F (from rev 1322, branches/land_ice/icesheet/src/module_read_netcdf.F)
===================================================================
--- branches/land_ice/icesheet/src/module_read_netcdf.F (rev 0)
+++ branches/land_ice/icesheet/src/module_read_netcdf.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,465 @@
+module read_netcdf
+
+ integer :: rd_ncid
+ integer :: rdDimIDTime
+ integer :: rdDimIDnCells
+ integer :: rdDimIDnEdges
+ integer :: rdDimIDnVertices
+ integer :: rdDimIDmaxEdges
+ integer :: rdDimIDmaxEdges2
+ integer :: rdDimIDnVertLevels
+ integer :: rdDimIDTWO
+ integer :: rdDimIDvertexDegree
+ integer :: rdVarIDlatCell
+ integer :: rdVarIDlonCell
+ integer :: rdVarIDxCell
+ integer :: rdVarIDyCell
+ integer :: rdVarIDzCell
+ integer :: rdVarIDindexToCellID
+ integer :: rdVarIDlatEdge
+ integer :: rdVarIDlonEdge
+ integer :: rdVarIDxEdge
+ integer :: rdVarIDyEdge
+ integer :: rdVarIDzEdge
+ integer :: rdVarIDindexToEdgeID
+ integer :: rdVarIDlatVertex
+ integer :: rdVarIDlonVertex
+ integer :: rdVarIDxVertex
+ integer :: rdVarIDyVertex
+ integer :: rdVarIDzVertex
+ integer :: rdVarIDindexToVertexID
+ integer :: rdVarIDcellsOnEdge
+ integer :: rdVarIDnEdgesOnCell
+ integer :: rdVarIDnEdgesOnEdge
+ integer :: rdVarIDedgesOnCell
+ integer :: rdVarIDedgesOnEdge
+ integer :: rdVarIDweightsOnEdge
+ integer :: rdVarIDdvEdge
+ integer :: rdVarIDdcEdge
+ integer :: rdVarIDangleEdge
+ integer :: rdVarIDareaCell
+ integer :: rdVarIDareaTriangle
+ integer :: rdVarIDcellsOnCell
+ integer :: rdVarIDverticesOnCell
+ integer :: rdVarIDverticesOnEdge
+ integer :: rdVarIDedgesOnVertex
+ integer :: rdVarIDcellsOnVertex
+ integer :: rdVarIDkiteAreasOnVertex
+ integer :: rdVarIDfEdge
+ integer :: rdVarIDfVertex
+
+ integer :: rdLocalnCells
+ integer :: rdLocalnEdges
+ integer :: rdLocalnVertices
+ integer :: rdLocalmaxEdges
+ integer :: rdLocalmaxEdges2
+ integer :: rdLocalnVertLevels
+ integer :: rdLocalTWO
+ integer :: rdLocalvertexDegree
+
+ contains
+
+ subroutine read_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ maxEdges2, &
+ nVertLevels, &
+ TWO, &
+ vertexDegree &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(out) :: nCells
+ integer, intent(out) :: nEdges
+ integer, intent(out) :: nVertices
+ integer, intent(out) :: maxEdges
+ integer, intent(out) :: maxEdges2
+ integer, intent(out) :: nVertLevels
+ integer, intent(out) :: TWO
+ integer, intent(out) :: vertexDegree
+
+ integer :: nferr
+
+
+ nferr = nf_open('grid.nc', NF_SHARE, rd_ncid)
+
+ !
+ ! Get IDs for variable dimensions
+ !
+ nferr = nf_inq_unlimdim(rd_ncid, rdDimIDTime)
+ nferr = nf_inq_dimid(rd_ncid, 'nCells', rdDimIDnCells)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnCells, rdLocalnCells)
+ nferr = nf_inq_dimid(rd_ncid, 'nEdges', rdDimIDnEdges)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnEdges, rdLocalnEdges)
+ nferr = nf_inq_dimid(rd_ncid, 'nVertices', rdDimIDnVertices)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertices, rdLocalnVertices)
+ nferr = nf_inq_dimid(rd_ncid, 'maxEdges', rdDimIDmaxEdges)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges, rdLocalmaxEdges)
+ nferr = nf_inq_dimid(rd_ncid, 'maxEdges2', rdDimIDmaxEdges2)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDmaxEdges2, rdLocalmaxEdges2)
+ nferr = nf_inq_dimid(rd_ncid, 'nVertLevels', rdDimIDnVertLevels)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnVertLevels, rdLocalnVertLevels)
+ nferr = nf_inq_dimid(rd_ncid, 'vertexDegree', rdDimIDvertexDegree)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDvertexDegree, rdLocalvertexDegree)
+ nferr = nf_inq_dimid(rd_ncid, 'TWO', rdDimIDTWO)
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDTWO, rdLocalTWO)
+
+
+ nCells = rdLocalnCells
+ nEdges = rdLocalnEdges
+ nVertices = rdLocalnVertices
+ maxEdges = rdLocalmaxEdges
+ maxEdges2 = rdLocalmaxEdges2
+ nVertLevels = rdLocalnVertLevels
+ vertexDegree = rdLocalvertexDegree
+ TWO = rdLocalTWO
+
+ !
+ ! Get IDs for variables
+ !
+ nferr = nf_inq_varid(rd_ncid, 'latCell', rdVarIDlatCell)
+ nferr = nf_inq_varid(rd_ncid, 'lonCell', rdVarIDlonCell)
+ nferr = nf_inq_varid(rd_ncid, 'xCell', rdVarIDxCell)
+ nferr = nf_inq_varid(rd_ncid, 'yCell', rdVarIDyCell)
+ nferr = nf_inq_varid(rd_ncid, 'zCell', rdVarIDzCell)
+ nferr = nf_inq_varid(rd_ncid, 'indexToCellID', rdVarIDindexToCellID)
+ nferr = nf_inq_varid(rd_ncid, 'latEdge', rdVarIDlatEdge)
+ nferr = nf_inq_varid(rd_ncid, 'lonEdge', rdVarIDlonEdge)
+ nferr = nf_inq_varid(rd_ncid, 'xEdge', rdVarIDxEdge)
+ nferr = nf_inq_varid(rd_ncid, 'yEdge', rdVarIDyEdge)
+ nferr = nf_inq_varid(rd_ncid, 'zEdge', rdVarIDzEdge)
+ nferr = nf_inq_varid(rd_ncid, 'indexToEdgeID', rdVarIDindexToEdgeID)
+ nferr = nf_inq_varid(rd_ncid, 'latVertex', rdVarIDlatVertex)
+ nferr = nf_inq_varid(rd_ncid, 'lonVertex', rdVarIDlonVertex)
+ nferr = nf_inq_varid(rd_ncid, 'xVertex', rdVarIDxVertex)
+ nferr = nf_inq_varid(rd_ncid, 'yVertex', rdVarIDyVertex)
+ nferr = nf_inq_varid(rd_ncid, 'zVertex', rdVarIDzVertex)
+ nferr = nf_inq_varid(rd_ncid, 'indexToVertexID', rdVarIDindexToVertexID)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnEdge', rdVarIDcellsOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'nEdgesOnCell', rdVarIDnEdgesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'nEdgesOnEdge', rdVarIDnEdgesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnCell', rdVarIDedgesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnEdge', rdVarIDedgesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'weightsOnEdge', rdVarIDweightsOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'dvEdge', rdVarIDdvEdge)
+ nferr = nf_inq_varid(rd_ncid, 'dcEdge', rdVarIDdcEdge)
+ nferr = nf_inq_varid(rd_ncid, 'angleEdge', rdVarIDangleEdge)
+ nferr = nf_inq_varid(rd_ncid, 'areaCell', rdVarIDareaCell)
+ nferr = nf_inq_varid(rd_ncid, 'areaTriangle', rdVarIDareaTriangle)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnCell', rdVarIDcellsOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'verticesOnCell', rdVarIDverticesOnCell)
+ nferr = nf_inq_varid(rd_ncid, 'verticesOnEdge', rdVarIDverticesOnEdge)
+ nferr = nf_inq_varid(rd_ncid, 'edgesOnVertex', rdVarIDedgesOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'cellsOnVertex', rdVarIDcellsOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'kiteAreasOnVertex', rdVarIDkiteAreasOnVertex)
+ nferr = nf_inq_varid(rd_ncid, 'fEdge', rdVarIDfEdge)
+ nferr = nf_inq_varid(rd_ncid, 'fVertex', rdVarIDfVertex)
+
+ end subroutine read_netcdf_init
+
+
+ subroutine read_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ xCell, &
+ yCell, &
+ zCell, &
+ indexToCellID, &
+ latEdge, &
+ lonEdge, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ indexToEdgeID, &
+ latVertex, &
+ lonVertex, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: time
+ real (kind=8), dimension(:), intent(out) :: latCell
+ real (kind=8), dimension(:), intent(out) :: lonCell
+ real (kind=8), dimension(:), intent(out) :: xCell
+ real (kind=8), dimension(:), intent(out) :: yCell
+ real (kind=8), dimension(:), intent(out) :: zCell
+ integer, dimension(:), intent(out) :: indexToCellID
+ real (kind=8), dimension(:), intent(out) :: latEdge
+ real (kind=8), dimension(:), intent(out) :: lonEdge
+ real (kind=8), dimension(:), intent(out) :: xEdge
+ real (kind=8), dimension(:), intent(out) :: yEdge
+ real (kind=8), dimension(:), intent(out) :: zEdge
+ integer, dimension(:), intent(out) :: indexToEdgeID
+ real (kind=8), dimension(:), intent(out) :: latVertex
+ real (kind=8), dimension(:), intent(out) :: lonVertex
+ real (kind=8), dimension(:), intent(out) :: xVertex
+ real (kind=8), dimension(:), intent(out) :: yVertex
+ real (kind=8), dimension(:), intent(out) :: zVertex
+ integer, dimension(:), intent(out) :: indexToVertexID
+ integer, dimension(:,:), intent(out) :: cellsOnEdge
+ integer, dimension(:), intent(out) :: nEdgesOnCell
+ integer, dimension(:), intent(out) :: nEdgesOnEdge
+ integer, dimension(:,:), intent(out) :: edgesOnCell
+ integer, dimension(:,:), intent(out) :: edgesOnEdge
+ real (kind=8), dimension(:,:), intent(out) :: weightsOnEdge
+ real (kind=8), dimension(:), intent(out) :: dvEdge
+ real (kind=8), dimension(:), intent(out) :: dcEdge
+ real (kind=8), dimension(:), intent(out) :: angleEdge
+ real (kind=8), dimension(:), intent(out) :: areaCell
+ real (kind=8), dimension(:), intent(out) :: areaTriangle
+ integer, dimension(:,:), intent(out) :: cellsOnCell
+ integer, dimension(:,:), intent(out) :: verticesOnCell
+ integer, dimension(:,:), intent(out) :: verticesOnEdge
+ integer, dimension(:,:), intent(out) :: edgesOnVertex
+ integer, dimension(:,:), intent(out) :: cellsOnVertex
+ real (kind=8), dimension(:,:), intent(out) :: kiteAreasOnVertex
+ real (kind=8), dimension(:), intent(out) :: fEdge
+ real (kind=8), dimension(:), intent(out) :: fVertex
+
+ integer :: nferr
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ start1(1) = 1
+
+ start2(1) = 1
+ start2(2) = 1
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatCell, start1, count1, latCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonCell, start1, count1, lonCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxCell, start1, count1, xCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyCell, start1, count1, yCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzCell, start1, count1, zCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToCellID, start1, count1, indexToCellID)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatEdge, start1, count1, latEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonEdge, start1, count1, lonEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxEdge, start1, count1, xEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyEdge, start1, count1, yEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzEdge, start1, count1, zEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToEdgeID, start1, count1, indexToEdgeID)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlatVertex, start1, count1, latVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDlonVertex, start1, count1, lonVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDxVertex, start1, count1, xVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDyVertex, start1, count1, yVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDzVertex, start1, count1, zVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDindexToVertexID, start1, count1, indexToVertexID)
+
+ start2(2) = 1
+ count2( 1) = rdLocalTWO
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnEdge, start2, count2, cellsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnCell, start2, count2, edgesOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges2
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnEdge, start2, count2, edgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges2
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDweightsOnEdge, start2, count2, weightsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDdvEdge, start1, count1, dvEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDdcEdge, start1, count1, dcEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDangleEdge, start1, count1, angleEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnCells
+ count1( 1) = rdLocalnCells
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDareaCell, start1, count1, areaCell)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDareaTriangle, start1, count1, areaTriangle)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnCell, start2, count2, cellsOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalmaxEdges
+ count2( 2) = rdLocalnCells
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnCell, start2, count2, verticesOnCell)
+
+ start2(2) = 1
+ count2( 1) = rdLocalTWO
+ count2( 2) = rdLocalnEdges
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDverticesOnEdge, start2, count2, verticesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_int(rd_ncid, rdVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = rdLocalvertexDegree
+ count2( 2) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnEdges
+ count1( 1) = rdLocalnEdges
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDfEdge, start1, count1, fEdge)
+
+ start1(1) = 1
+ count1( 1) = rdLocalnVertices
+ count1( 1) = rdLocalnVertices
+ nferr = nf_get_vara_double(rd_ncid, rdVarIDfVertex, start1, count1, fVertex)
+
+ end subroutine read_netcdf_fields
+
+
+ subroutine read_netcdf_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(rd_ncid)
+
+ end subroutine read_netcdf_finalize
+
+end module read_netcdf
Copied: branches/land_ice/icesheet/src/module_read_topo.F (from rev 1322, branches/land_ice/icesheet/src/module_read_topo.F)
===================================================================
--- branches/land_ice/icesheet/src/module_read_topo.F (rev 0)
+++ branches/land_ice/icesheet/src/module_read_topo.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,138 @@
+module read_topo
+
+ integer :: rd_ncid
+ integer :: rdDimIDnx
+ integer :: rdDimIDny
+ integer :: rdVarIDx
+ integer :: rdVarIDy
+ integer :: rdVarIDthk
+ integer :: rdVarIDtopg
+
+
+ integer :: rdLocalnx
+ integer :: rdLocalny
+
+ contains
+
+ subroutine read_topo_init( nx, ny )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(out) :: nx, ny
+
+ integer :: nferr
+
+ integer, dimension(3) :: dimids
+
+
+! nferr = nf_open('topo/Greenland_5km_v1.1.nc', NF_SHARE, rd_ncid)
+! nferr = nf_open('topo/gis_5km.180511.nc', NF_SHARE, rd_ncid)
+ nferr = nf_open('topo/gis_20km.180511.nc', NF_SHARE, rd_ncid)
+ write(6,*) ' nferr ncid ', nferr, rd_ncid
+
+ !
+ ! Get IDs for variable dimensions
+ !
+
+ nferr = nf_inq_dimid(rd_ncid, 'x1', rdDimIDnx)
+ write(6,*) ' nferr xdim-id ', nferr, rdDimIDnx
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDnx, rdLocalnx)
+ write(6,*) ' nferr xdim-length ', nferr, rdLocalnx
+
+ nferr = nf_inq_dimid(rd_ncid, 'y1', rdDimIDny)
+ write(6,*) ' nferr ydim-id ', nferr, rdDimIDny
+ nferr = nf_inq_dimlen(rd_ncid, rdDimIDny, rdLocalny)
+ write(6,*) ' nferr ydim-length ', nferr, rdLocalny
+
+ nx = rdLocalnx
+ ny = rdLocalny
+ print *, 'x,y dims of .nc var = ', nx, ny
+
+ !
+ ! Get IDs for variables
+ !
+ 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,topg)
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ real (kind=4), dimension(:), intent(out) :: x,y
+ real (kind=4), dimension(:,:), intent(out) :: thk, topg
+
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ integer :: nferr
+
+ start1(1) = 1
+ count1(1) = rdLocalnx
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDx, start1, count1, x)
+ write(6,*) ' nferr ncid xvar-id (get x-val) ', nferr, rd_ncid, rdVarIDx
+
+ start1(1) = 1
+ count1(1) = rdLocalny
+ nferr = nf_get_vara_real(rd_ncid, rdVarIDy, start1, count1, y)
+ write(6,*) ' nferr ncid yvar-id (get y-val) ', nferr, rd_ncid, rdVarIDy
+
+ start2(1) = 1
+ start2(2) = 1
+ count2(1) = rdLocalnx
+ count2(2) = rdLocalny
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+ count3(1) = rdlocalnx
+ count3(2) = rdlocalny
+ count3(3) = 1
+
+ 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
+
+
+ subroutine read_topo_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(rd_ncid)
+ write(6,*) ' nferr ', nferr
+
+
+ end subroutine read_topo_finalize
+
+end module read_topo
Copied: branches/land_ice/icesheet/src/module_write_netcdf.F (from rev 1322, branches/land_ice/icesheet/src/module_write_netcdf.F)
===================================================================
--- branches/land_ice/icesheet/src/module_write_netcdf.F (rev 0)
+++ branches/land_ice/icesheet/src/module_write_netcdf.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,566 @@
+module write_netcdf
+
+ integer :: wr_ncid
+ integer :: wrDimIDTime
+ integer :: wrDimIDnCells
+ integer :: wrDimIDnEdges
+ integer :: wrDimIDnVertices
+ integer :: wrDimIDmaxEdges
+ integer :: wrDimIDmaxEdges2
+ integer :: wrDimIDTWO
+ integer :: wrDimIDvertexDegree
+ integer :: wrDimIDnVertLevels
+ integer :: wrVarIDlatCell
+ integer :: wrVarIDlonCell
+ integer :: wrVarIDxCell
+ integer :: wrVarIDyCell
+ integer :: wrVarIDzCell
+ integer :: wrVarIDindexToCellID
+ integer :: wrVarIDlatEdge
+ integer :: wrVarIDlonEdge
+ integer :: wrVarIDxEdge
+ integer :: wrVarIDyEdge
+ integer :: wrVarIDzEdge
+ integer :: wrVarIDindexToEdgeID
+ integer :: wrVarIDlatVertex
+ integer :: wrVarIDlonVertex
+ integer :: wrVarIDxVertex
+ integer :: wrVarIDyVertex
+ integer :: wrVarIDzVertex
+ integer :: wrVarIDindexToVertexID
+ integer :: wrVarIDmaxLevelCell
+ integer :: wrVarIDcellsOnEdge
+ integer :: wrVarIDnEdgesOnCell
+ integer :: wrVarIDnEdgesOnEdge
+ integer :: wrVarIDedgesOnCell
+ integer :: wrVarIDedgesOnEdge
+ integer :: wrVarIDweightsOnEdge
+ integer :: wrVarIDdvEdge
+ integer :: wrVarIDdcEdge
+ integer :: wrVarIDangleEdge
+ integer :: wrVarIDareaCell
+ integer :: wrVarIDareaTriangle
+ integer :: wrVarIDcellsOnCell
+ integer :: wrVarIDverticesOnCell
+ integer :: wrVarIDverticesOnEdge
+ integer :: wrVarIDedgesOnVertex
+ integer :: wrVarIDcellsOnVertex
+ integer :: wrVarIDkiteAreasOnVertex
+ integer :: wrVarIDfEdge
+ integer :: wrVarIDfVertex
+ integer :: wrVarIDboundaryEdge
+ integer :: wrVarIDboundaryVertex
+
+ integer :: wrVarIDh
+
+ ! ice sheet variables
+ integer :: wrVarIDthck
+ integer :: wrVarIDtopg
+
+ integer :: wrLocalnCells
+ integer :: wrLocalnEdges
+ integer :: wrLocalnVertices
+ integer :: wrLocalmaxEdges
+ integer :: wrLocalnVertLevels
+ integer :: wrLocalvertexDegree
+
+ contains
+
+ subroutine write_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ nVertLevels, &
+ vertexDegree, &
+ sphere_radius, &
+ on_a_sphere &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertices
+ integer, intent(in) :: maxEdges
+ integer, intent(in) :: nVertLevels
+ integer, intent(in) :: vertexDegree
+ character (len=16) :: on_a_sphere
+ double precision :: sphere_radius
+
+
+ integer :: nferr
+ integer, dimension(10) :: dimlist
+
+
+ wrLocalnCells = nCells
+ wrLocalnEdges = nEdges
+ wrLocalnVertices = nVertices
+ wrLocalmaxEdges = maxEdges
+ wrLocalnVertLevels = nVertLevels
+ wrLocalvertexDegree = vertexDegree
+
+ nferr = nf_create('icesheet.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
+
+ !
+ ! Define dimensions
+ !
+ nferr = nf_def_dim(wr_ncid, 'nCells', nCells, wrDimIDnCells)
+ nferr = nf_def_dim(wr_ncid, 'nEdges', nEdges, wrDimIDnEdges)
+ nferr = nf_def_dim(wr_ncid, 'nVertices', nVertices, wrDimIDnVertices)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges', maxEdges, wrDimIDmaxEdges)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges2', 2*maxEdges, wrDimIDmaxEdges2)
+ nferr = nf_def_dim(wr_ncid, 'TWO', 2, wrDimIDTWO)
+ nferr = nf_def_dim(wr_ncid, 'vertexDegree', vertexDegree, wrDimIDvertexDegree)
+ nferr = nf_def_dim(wr_ncid, 'nVertLevels', nVertLevels, wrDimIDnVertLevels)
+ nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
+
+ !
+ ! Define variables
+ !
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'latCell', NF_DOUBLE, 1, dimlist, wrVarIDlatCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'lonCell', NF_DOUBLE, 1, dimlist, wrVarIDlonCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'xCell', NF_DOUBLE, 1, dimlist, wrVarIDxCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'yCell', NF_DOUBLE, 1, dimlist, wrVarIDyCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'zCell', NF_DOUBLE, 1, dimlist, wrVarIDzCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'indexToCellID', NF_INT, 1, dimlist, wrVarIDindexToCellID)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'latEdge', NF_DOUBLE, 1, dimlist, wrVarIDlatEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'lonEdge', NF_DOUBLE, 1, dimlist, wrVarIDlonEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'xEdge', NF_DOUBLE, 1, dimlist, wrVarIDxEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'yEdge', NF_DOUBLE, 1, dimlist, wrVarIDyEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'zEdge', NF_DOUBLE, 1, dimlist, wrVarIDzEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'indexToEdgeID', NF_INT, 1, dimlist, wrVarIDindexToEdgeID)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'latVertex', NF_DOUBLE, 1, dimlist, wrVarIDlatVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'lonVertex', NF_DOUBLE, 1, dimlist, wrVarIDlonVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'xVertex', NF_DOUBLE, 1, dimlist, wrVarIDxVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'yVertex', NF_DOUBLE, 1, dimlist, wrVarIDyVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'zVertex', NF_DOUBLE, 1, dimlist, wrVarIDzVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'indexToVertexID', NF_INT, 1, dimlist, wrVarIDindexToVertexID)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'maxLevelCell', NF_INT, 1, dimlist, wrVarIDmaxLevelCell)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'cellsOnEdge', NF_INT, 2, dimlist, wrVarIDcellsOnEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnCell', NF_INT, 1, dimlist, wrVarIDnEdgesOnCell)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnEdge', NF_INT, 1, dimlist, wrVarIDnEdgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'edgesOnCell', NF_INT, 2, dimlist, wrVarIDedgesOnCell)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'edgesOnEdge', NF_INT, 2, dimlist, wrVarIDedgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'weightsOnEdge', NF_DOUBLE, 2, dimlist, wrVarIDweightsOnEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dvEdge', NF_DOUBLE, 1, dimlist, wrVarIDdvEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dcEdge', NF_DOUBLE, 1, dimlist, wrVarIDdcEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'angleEdge', NF_DOUBLE, 1, dimlist, wrVarIDangleEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'areaCell', NF_DOUBLE, 1, dimlist, wrVarIDareaCell)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'areaTriangle', NF_DOUBLE, 1, dimlist, wrVarIDareaTriangle)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'cellsOnCell', NF_INT, 2, dimlist, wrVarIDcellsOnCell)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'verticesOnCell', NF_INT, 2, dimlist, wrVarIDverticesOnCell)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'verticesOnEdge', NF_INT, 2, dimlist, wrVarIDverticesOnEdge)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'edgesOnVertex', NF_INT, 2, dimlist, wrVarIDedgesOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'cellsOnVertex', NF_INT, 2, dimlist, wrVarIDcellsOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'kiteAreasOnVertex', NF_DOUBLE, 2, dimlist, wrVarIDkiteAreasOnVertex)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'fEdge', NF_DOUBLE, 1, dimlist, wrVarIDfEdge)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'fVertex', NF_DOUBLE, 1, dimlist, wrVarIDfVertex)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'boundaryEdge', NF_INT, 2, dimlist, wrVarIDboundaryEdge)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'boundaryVertex', NF_INT, 2, dimlist, wrVarIDboundaryVertex)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+
+ nferr = nf_def_var(wr_ncid, 'h', NF_DOUBLE, 3, dimlist, wrVarIDh)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnCells
+ dimlist( 3) = wrDimIDTime
+
+ ! ice sheet variables (NOTE: for now, we are assuming vert dim of 1, so this works). If nVertLevels > 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)
+
+ nferr = nf_enddef(wr_ncid)
+
+ end subroutine write_netcdf_init
+
+
+ subroutine write_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ xCell, &
+ yCell, &
+ zCell, &
+ indexToCellID, &
+ latEdge, &
+ lonEdge, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ indexToEdgeID, &
+ latVertex, &
+ lonVertex, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ indexToVertexID, &
+ maxLevelCell, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex, &
+ boundaryEdge, &
+ boundaryVertex, &
+ h, &
+ topg, &
+ thck &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: time
+ real (kind=8), dimension(:), intent(in) :: latCell
+ real (kind=8), dimension(:), intent(in) :: lonCell
+ real (kind=8), dimension(:), intent(in) :: xCell
+ real (kind=8), dimension(:), intent(in) :: yCell
+ real (kind=8), dimension(:), intent(in) :: zCell
+ integer, dimension(:), intent(in) :: indexToCellID
+ real (kind=8), dimension(:), intent(in) :: latEdge
+ real (kind=8), dimension(:), intent(in) :: lonEdge
+ real (kind=8), dimension(:), intent(in) :: xEdge
+ real (kind=8), dimension(:), intent(in) :: yEdge
+ real (kind=8), dimension(:), intent(in) :: zEdge
+ integer, dimension(:), intent(in) :: indexToEdgeID
+ real (kind=8), dimension(:), intent(in) :: latVertex
+ real (kind=8), dimension(:), intent(in) :: lonVertex
+ real (kind=8), dimension(:), intent(in) :: xVertex
+ real (kind=8), dimension(:), intent(in) :: yVertex
+ real (kind=8), dimension(:), intent(in) :: zVertex
+ integer, dimension(:), intent(in) :: indexToVertexID
+ integer, dimension(:), intent(in) :: maxLevelCell
+ integer, dimension(:,:), intent(in) :: cellsOnEdge
+ integer, dimension(:), intent(in) :: nEdgesOnCell
+ integer, dimension(:), intent(in) :: nEdgesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnCell
+ integer, dimension(:,:), intent(in) :: edgesOnEdge
+ real (kind=8), dimension(:,:), intent(in) :: weightsOnEdge
+ real (kind=8), dimension(:), intent(in) :: dvEdge
+ real (kind=8), dimension(:), intent(in) :: dcEdge
+ real (kind=8), dimension(:), intent(in) :: angleEdge
+ real (kind=8), dimension(:), intent(in) :: areaCell
+ real (kind=8), dimension(:), intent(in) :: areaTriangle
+ integer, dimension(:,:), intent(in) :: cellsOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnVertex
+ integer, dimension(:,:), intent(in) :: cellsOnVertex
+ real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
+ real (kind=8), dimension(:), intent(in) :: fEdge
+ real (kind=8), dimension(:), intent(in) :: fVertex
+ integer, dimension(:,:), intent(in) :: boundaryEdge
+ integer, dimension(:,:), intent(in) :: boundaryVertex
+
+ real (kind=8), dimension(:,:,:), intent(in) :: h
+
+ ! ice sheet fields
+ real (kind=8), dimension(:,:,:), intent(in) :: thck
+ real (kind=8), dimension(:,:,:), intent(in) :: topg
+
+ integer :: nferr
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ start1(1) = 1
+
+ start2(1) = 1
+ start2(2) = 1
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+
+ start4(1) = 1
+ start4(2) = 1
+ start4(3) = 1
+ start4(4) = 1
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatCell, start1, count1, latCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonCell, start1, count1, lonCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxCell, start1, count1, xCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyCell, start1, count1, yCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzCell, start1, count1, zCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToCellID, start1, count1, indexToCellID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatEdge, start1, count1, latEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonEdge, start1, count1, lonEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxEdge, start1, count1, xEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyEdge, start1, count1, yEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzEdge, start1, count1, zEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToEdgeID, start1, count1, indexToEdgeID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatVertex, start1, count1, latVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonVertex, start1, count1, lonVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxVertex, start1, count1, xVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyVertex, start1, count1, yVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzVertex, start1, count1, zVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToVertexID, start1, count1, indexToVertexID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDmaxLevelCell, start1, count1, maxLevelCell)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnEdge, start2, count2, cellsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnCell, start2, count2, edgesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnEdge, start2, count2, edgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDweightsOnEdge, start2, count2, weightsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdvEdge, start1, count1, dvEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdcEdge, start1, count1, dcEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDangleEdge, start1, count1, angleEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaCell, start1, count1, areaCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaTriangle, start1, count1, areaTriangle)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnCell, start2, count2, cellsOnCell)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnCell, start2, count2, verticesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnEdge, start2, count2, verticesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalvertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfEdge, start1, count1, fEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfVertex, start1, count1, fVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryEdge, start2, count2, boundaryEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryVertex, start2, count2, boundaryVertex)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnCells
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDh, start3, count3, h)
+
+ ! 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
+
+
+ subroutine write_netcdf_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(wr_ncid)
+
+ end subroutine write_netcdf_finalize
+
+end module write_netcdf
Copied: branches/land_ice/icesheet/src/utilities.F (from rev 1322, branches/land_ice/icesheet/src/utilities.F)
===================================================================
--- branches/land_ice/icesheet/src/utilities.F (rev 0)
+++ branches/land_ice/icesheet/src/utilities.F 2012-01-09 22:19:16 UTC (rev 1329)
@@ -0,0 +1,771 @@
+module utilities
+
+contains
+
+subroutine write_OpenDX( on_a_sphere, &
+ nCells, &
+ nVertices, &
+ nEdges, &
+ vertexDegree, &
+ maxEdges, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ cellsOnVertex, &
+ edgesOnCell, &
+ areaCell, &
+ maxLevelCell, &
+ depthCell, &
+ kiteAreasOnVertex )
+
+ implicit none
+
+ character (len=16), intent(in) :: on_a_sphere
+ integer, intent(in) :: nCells, nVertices, vertexDegree, nEdges, maxEdges
+ real (kind=8), dimension(nCells), intent(inout) :: xCell
+ real (kind=8), dimension(nCells), intent(inout) :: yCell
+ real (kind=8), dimension(nCells), intent(inout) :: zCell
+ real (kind=8), dimension(nVertices), intent(inout) :: xVertex
+ real (kind=8), dimension(nVertices), intent(inout) :: yVertex
+ real (kind=8), dimension(nVertices), intent(inout) :: zVertex
+ real (kind=8), dimension(nEdges), intent(inout) :: xEdge
+ real (kind=8), dimension(nEdges), intent(inout) :: yEdge
+ real (kind=8), dimension(nEdges), intent(inout) :: zEdge
+ integer, dimension(nCells), intent(in) :: nEdgesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: edgesOnCell
+ integer, dimension(2,nEdges), intent(in) :: verticesOnEdge
+ integer, dimension(vertexDegree, nVertices), intent(in) :: cellsOnVertex
+ integer, dimension(nCells), intent(in) :: maxLevelCell
+ real (kind=8), dimension(nCells), intent(in) :: areaCell
+ real (kind=8), dimension(nCells), intent(in) :: depthCell
+ real (kind=8), dimension(vertexDegree,nVertices), intent(in) :: kiteAreasOnVertex
+
+ character(len=80) :: a, b, c, d, e, f
+ integer :: i, j, k, nVerticesTotal, iEdge, iLoop, iFace, Vert(4), Edge(4), iVertex, i1, i2, jp1
+ integer :: nKitesTotal, iCell, iEdge1, iEdge2, iVertex11, iVertex12, iVertex21, iVertex22, ksave
+ real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xscale, work(nCells)
+ real (kind=8) :: xv, yv, zv, xc, yc, zc, dist
+ logical (kind=8) :: eflag
+
+ if(on_a_sphere.eq.'NO ') then
+ write(6,*) ' write_dx, not on a sphere '
+ endif
+
+ xscale = 1.00
+ xCell = xCell*xscale
+ yCell = yCell*xscale
+ zCell = zCell*xscale
+ xVertex = xVertex*xscale
+ yVertex = yVertex*xscale
+ zVertex = zVertex*xscale
+ xEdge = xEdge*xscale
+ yEdge = yEdge*xscale
+ zEdge = zEdge*xscale
+
+ write(6,*) 'xCell', minval(xCell), maxval(xCell)
+ write(6,*) ' nCells', nCells
+ write(6,*) ' nEdges', nEdges
+ write(6,*) ' nVertices', nVertices
+ write(6,*) ' nEdgesOnCell',minval(nEdgesOnCell), maxval(nEdgesOnCell)
+
+ open(unit=1,file='dx/vector.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file vector.position.data')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,*)
+
+ a = trim('object 0 class array type float rank 1 shape 3 items')
+ b = trim('ascii data file vector.data')
+ c = trim('attribute "dep" string "positions"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "vector" class field')
+ b = trim('component "positions" "positions list"')
+ c = trim('component "data" 0')
+ write(1,10) a
+ write(1,10) b
+ write(1,10) c
+
+ close(1)
+
+ open(unit=14,file='dx/vector.position.data',form='formatted',status='unknown')
+ do i=1,nCells
+ write(14,22) xCell(i), yCell(i), zCell(i)
+ enddo
+ close(14)
+
+
+
+ nVerticesTotal = 0
+ do i=1,nCells
+ nVerticesTotal = nVerticesTotal + nEdgesOnCell(i)
+ enddo
+ write(6,*) 'total number of vertices', nVerticesTotal
+
+ open(unit=1,file='dx/ocean.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file ocean.position.data')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,*)
+ 10 format(a70,i10)
+
+ a = trim('object "edge list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.edge.data')
+ c = trim('attribute "ref" string "positions"')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "loops list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.loop.data')
+ c = trim('attribute "ref" string "edges"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "face list" class array type int rank 0 items')
+ b = trim('ascii data file ocean.face.data')
+ c = trim('attribute "ref" string "loops"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object 0 class array type float rank 0 items')
+ b = trim('data file ocean.area.data')
+ c = trim('attribute "dep" string "faces"')
+ write(1,10) a, nCells
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "area" class field')
+ b = trim('component "positions" "positions list"')
+ c = trim('component "edges" "edge list"')
+ d = trim('component "loops" "loops list"')
+ e = trim('component "faces" "face list"')
+ f = trim('component "data" 0')
+ write(1,10) a
+ write(1,10) b
+ write(1,10) c
+ write(1,10) d
+ write(1,10) e
+ write(1,10) f
+
+ close(1)
+
+ work = depthCell
+
+ open(unit=10,file='dx/ocean.area.data',form='formatted',status='unknown')
+ open(unit=11,file='dx/ocean.face.data',form='formatted',status='unknown')
+ open(unit=12,file='dx/ocean.loop.data',form='formatted',status='unknown')
+ open(unit=13,file='dx/ocean.edge.data',form='formatted',status='unknown')
+ open(unit=14,file='dx/ocean.position.data',form='formatted',status='unknown')
+
+ iLoop = 0
+ iEdge = 0
+ do i=1,nCells
+ write(10,20) work(i)
+ write(11,21) i-1
+ write(12,21) iLoop
+ iLoop = iLoop + nEdgesOnCell(i)
+
+ eflag = .false.
+ do j=1,nEdgesOnCell(i)
+ k = verticesOnCell(j,i)
+ xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
+ xc = xCell(i); yc = yCell(i); zc = zCell(i)
+ dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
+ if(dist.gt.5.0e5.and.on_a_sphere.eq.'NO ') then
+ eflag = .true.
+ endif
+ enddo
+
+ if(eflag) then
+
+ do j=1,nEdgesOnCell(i)
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ k = verticesOnCell(j,i)
+ xv = xVertex(k); yv = yVertex(k); zv = zVertex(k)
+ xc = xCell(i); yc = yCell(i); zc = zCell(i)
+ dist = sqrt( (xc-xv)**2 + (yc-yv)**2 + (zc-zv)**2 )
+ if(dist.gt.5.0e5) then
+ write(14,22) xc, yc, zc
+ else
+ write(14,22) xv, yv, zv
+ endif
+ enddo
+
+ else
+
+ do j=1,nEdgesOnCell(i)
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ k = verticesOnCell(j,i)
+ if(k.le.0) write(6,*) ' vert1 ',k, verticesOnCell(:,i)
+ write(14,22) xVertex(k), yVertex(k), zVertex(k)
+ write(15,23) j,i,k,xVertex(k), yVertex(k), zVertex(k)
+ enddo
+ endif
+ enddo
+
+ 20 format(e20.10)
+ 21 format(i20)
+ 22 format(3e20.10)
+ 23 format(3i6, 3e20.10)
+
+ close(10)
+ close(11)
+ close(12)
+ close(13)
+ close(14)
+
+ ! nVerticesTotal = 0
+ ! nKitesTotal = 0
+ ! do i=1,nCells
+ ! nKitesTotal = nKitesTotal + nEdgesOnCell(i)
+ ! enddo
+ ! nVerticesTotal = nKitesTotal*4
+ ! write(6,*) nKitesTotal, nVerticesTotal
+
+ ! open(unit=1,file='dx/kite.dx',form='formatted',status='unknown')
+
+ ! a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ ! b = trim('ascii data file kite.position.data')
+ ! write(1,10) a, nVerticesTotal
+ ! write(1,10) b
+ ! write(1,*)
+
+ ! a = trim('object "edge list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.edge.data')
+ ! c = trim('attribute "ref" string "positions"')
+ ! write(1,10) a, nVerticesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "loops list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.loop.data')
+ ! c = trim('attribute "ref" string "edges"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "face list" class array type int rank 0 items')
+ ! b = trim('ascii data file kite.face.data')
+ ! c = trim('attribute "ref" string "loops"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object 0 class array type float rank 0 items')
+ ! b = trim('data file kite.area.data')
+ ! c = trim('attribute "dep" string "faces"')
+ ! write(1,10) a, nKitesTotal
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,*)
+
+ ! a = trim('object "area" class field')
+ ! b = trim('component "positions" "positions list"')
+ ! c = trim('component "edges" "edge list"')
+ ! d = trim('component "loops" "loops list"')
+ ! e = trim('component "faces" "face list"')
+ ! f = trim('component "data" 0')
+ ! write(1,10) a
+ ! write(1,10) b
+ ! write(1,10) c
+ ! write(1,10) d
+ ! write(1,10) e
+ ! write(1,10) f
+
+ ! close(1)
+
+ ! open(unit=10,file='dx/kite.area.data',form='formatted',status='unknown')
+ ! open(unit=11,file='dx/kite.face.data',form='formatted',status='unknown')
+ ! open(unit=12,file='dx/kite.loop.data',form='formatted',status='unknown')
+ ! open(unit=13,file='dx/kite.edge.data',form='formatted',status='unknown')
+ ! open(unit=14,file='dx/kite.position.data',form='formatted',status='unknown')
+
+ ! iLoop = 0
+ ! iEdge = 0
+ ! iFace = 0
+
+ ! do iCell=1,nCells
+ ! do j=1,nEdgesOnCell(iCell)
+ ! iEdge1 = edgesOnCell(j,iCell)
+ ! jp1 = j+1
+ ! if(j.eq.nEdgesOnCell(iCell)) jp1=1
+ ! iEdge2 = edgesOnCell(jp1,iCell)
+
+ ! iVertex11 = verticesOnEdge(1,iEdge1)
+ ! iVertex21 = verticesOnEdge(2,iEdge1)
+ ! iVertex12 = verticesOnEdge(1,iEdge2)
+ ! ivertex22 = verticesOnEdge(2,iEdge2)
+
+ ! if(iVertex11.eq.iVertex12.or.iVertex11.eq.iVertex22) then
+ ! iVertex = iVertex11
+ ! elseif(iVertex21.eq.iVertex12.or.iVertex21.eq.iVertex22) then
+ ! iVertex = iVertex21
+ ! else
+ ! write(6,*) iVertex11, iVertex21, iVertex12, iVertex22
+ ! stop
+ ! endif
+
+ ! ksave = 0
+ ! do k=1,vertexDegree
+ ! if(cellsOnVertex(k,iVertex).eq.iCell) ksave=k
+ ! enddo
+ ! if(ksave.eq.0) then
+ ! write(6,*) ' can not find iCell'
+ ! write(6,*) cellsOnVertex(:,iVertex)
+ ! write(6,*) iCell
+ ! write(6,*) iEdge1, iEdge2
+ ! write(6,*) iVertex11, iVertex21, iVertex21, iVertex22
+ ! write(6,*) iVertex
+ ! stop
+ ! endif
+
+ ! write(11,21) iFace
+ ! write(12,21) iLoop
+ ! iFace = iFace + 1
+ ! iLoop = iLoop + 4
+ ! do k=1,4
+ ! write(13,21) iEdge
+ ! iEdge = iEdge + 1
+ ! enddo
+ !
+ ! x1 = xCell(iCell) ; y1 = yCell(iCell) ; z1 = zCell(iCell)
+ ! x2 = xEdge(iEdge1) ; y2 = yEdge(iEdge1) ; z2 = zEdge(iEdge1)
+ ! x3 = xVertex(iVertex); y3 = yVertex(iVertex); z3 = zVertex(iVertex)
+ ! x4 = xEdge(iEdge2) ; y4 = yEdge(iEdge2) ; z4 = zEdge(iEdge2)
+ !
+ ! write(14,22) x1, y1, z1
+ ! write(14,22) x2, y2, z2
+ ! write(14,22) x3, y3, z3
+ ! write(14,22) x4, y4, z4
+ ! write(10,22) kiteAreasOnVertex(ksave,iVertex)
+
+ ! enddo
+ ! enddo
+
+end subroutine write_OpenDX
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE CONVERT_LX
+!
+! Convert (lat,lon) to an (x, y, z) location on a sphere with specified radius.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine convert_lx(x, y, z, radius, lat, lon)
+
+ implicit none
+
+ real, intent(in) :: radius
+ real, intent(in) :: lat, lon
+ real, intent(out) :: x, y, z
+
+ z = radius * sin(lat)
+ x = radius * cos(lon) * cos(lat)
+ y = radius * sin(lon) * cos(lat)
+
+end subroutine convert_lx
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE CONVERT_XL
+!
+! Convert (x, y, z) to a (lat, lon) location on a sphere with
+! radius sqrt(x^2 + y^2 + z^2).
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine convert_xl(x, y, z, lat,lon)
+
+ implicit none
+
+ real, intent(in) :: x, y, z
+ real, intent(out) :: lat, lon
+
+ real :: dl, clat, pii, rtod
+ real :: eps
+ parameter (eps=1.e-10)
+
+ pii = 2.*asin(1.0)
+ rtod=180./pii
+ dl = sqrt(x*x + y*y + z*z)
+
+ lat = asin(z/dl)
+
+! check for being close to either pole
+
+ if (abs(x) > eps) then
+
+ if (abs(y) > eps) then
+
+ lon = atan(abs(y/x))
+
+ if ((x <= 0.) .and. (y >= 0.)) then
+ lon = pii-lon
+ else if ((x <= 0.) .and. (y < 0.)) then
+ lon = lon+pii
+ else if ((x >= 0.) .and. (y <= 0.)) then
+ lon = 2*pii-lon
+ end if
+
+ else ! we're either on longitude 0 or 180
+
+ if (x > 0) then
+ lon = 0.
+ else
+ lon = pii
+ end if
+
+ end if
+
+ else if (abs(y) > eps) then
+
+ if (y > 0) then
+ lon = pii/2.
+ else
+ lon = 3.*pii/2.
+ end if
+
+ else ! we are at a pole
+
+ lon = 0.
+
+ end if
+
+end subroutine convert_xl
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine transform_from_lonlat_to_xyz(xin, yin, zin, ulon, ulat, ux, uy, uz)
+!
+! transform vector measured in latitude/longitude space to a vector measured in x,y,z
+!
+! INTENT(IN)
+! xin = x position
+! yin = y position
+! zin = z position
+! ulon = east component of vector
+! ulat = north component of vector
+!
+! INTENT(OUT)
+! ux = x component of vector
+! uy = y component of vector
+! uz = z component of vector
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+real, intent(in) :: xin, yin, zin, ulon, ulat
+real, intent(out) :: ux, uy, uz
+real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X2_to_X1(3,3), r
+integer :: i,j,k
+logical :: l_Pole
+real, parameter :: epsvt = 1.0e-10
+
+!-----------------------------------------------------------------------
+! define the e1, e2, and e3 directions
+!-----------------------------------------------------------------------
+ X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
+ X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
+ X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
+
+!-----------------------------------------------------------------------
+! find the vectors (measured in X1) that point in the local
+! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
+!-----------------------------------------------------------------------
+ h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
+ call unit_vector_in_3space(h(3,:))
+
+!-----------------------------------------------------------------------
+! g(:) is a work array and holds the vector pointing to the North Pole.
+! measured in X1
+!-----------------------------------------------------------------------
+ g(:) = X1(3,:)
+
+!-----------------------------------------------------------------------
+! determine if the local vertical hits a pole
+!-----------------------------------------------------------------------
+ l_Pole = .false.
+ r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
+ r = abs(r) + epsvt
+ if(r.gt.1.0) then
+ l_Pole = .true.
+ h(3,:) = h(3,:) + epsvt
+ call unit_vector_in_3space(h(3,:))
+ endif
+
+!-----------------------------------------------------------------------
+! find the vector that is perpendicular to the local vertical vector
+! and points in the direction of of the North pole, this defines the local
+! north direction. measured in X1
+!-----------------------------------------------------------------------
+ call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
+
+!-----------------------------------------------------------------------
+! take the cross product of the local North direction and the local vertical
+! to find the local east vector. still in X1
+!-----------------------------------------------------------------------
+ call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
+
+!-----------------------------------------------------------------------
+! put these 3 vectors into a matrix X2
+!-----------------------------------------------------------------------
+ X2(1,:) = h(1,:) ! local east (measured in X1)
+ X2(2,:) = h(2,:) ! local north (measured in X1)
+ X2(3,:) = h(3,:) ! local vertical (measured in X1)
+
+!-----------------------------------------------------------------------
+! compute the transformation matrix
+!-----------------------------------------------------------------------
+ trans_X2_to_X1(:,:) = matmul(X1,transpose(X2))
+
+!-----------------------------------------------------------------------
+! transform (ulon, ulat) into (x,y,z)
+!-----------------------------------------------------------------------
+ p(1) = ulon; p(2) = ulat; p(3) = 0
+ g(:) = matmul(trans_X2_to_X1(:, :), p(:))
+ ux = g(1); uy = g(2); uz = g(3)
+
+end subroutine transform_from_lonlat_to_xyz
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine transform_from_xyz_to_lonlat(xin, yin, zin, ux, uy, uz, ulon, ulat)
+!
+! transform vector measured in x,y,z space to a vector measured in latitude/longitude space
+!
+! INTENT(IN)
+! xin = x position
+! yin = y position
+! zin = z position
+! ux = x component of vector
+! uy = y component of vector
+! uz = z component of vector
+!
+! INTENT(OUT)
+! ulon = east component of vector
+! ulat = north component of vector
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+real, intent(in) :: xin, yin, zin, ux, uy, uz
+real, intent(out) :: ulon, ulat
+real :: h(3,3), p(3), q(3), g(3), X1(3,3), X2(3,3), trans_X1_to_X2(3,3), r
+integer :: i,j,k
+logical :: l_Pole
+real, parameter :: epsvt = 1.0e-10
+
+!-----------------------------------------------------------------------
+! define the e1, e2, and e3 directions
+!-----------------------------------------------------------------------
+ X1(1,1) = 1.0; X1(1,2) = 0.0; X1(1,3) = 0.0
+ X1(2,1) = 0.0; X1(2,2) = 1.0; X1(2,3) = 0.0
+ X1(3,1) = 0.0; X1(3,2) = 0.0; X1(3,3) = 1.0
+
+!-----------------------------------------------------------------------
+! find the vectors (measured in X1) that point in the local
+! east (h(1,:)), north (h(2,:)), and vertical (h(3,:)) direction
+!-----------------------------------------------------------------------
+ h(3,1) = xin; h(3,2) = yin; h(3,3) = zin
+ call unit_vector_in_3space(h(3,:))
+
+!-----------------------------------------------------------------------
+! g(:) is a work array and holds the vector pointing to the North Pole.
+! measured in X1
+!-----------------------------------------------------------------------
+ g(:) = X1(3,:)
+
+!-----------------------------------------------------------------------
+! determine if the local vertical hits a pole
+!-----------------------------------------------------------------------
+ l_Pole = .false.
+ r = g(1)*h(3,1) + g(2)*h(3,2) + g(3)*h(3,3)
+ r = abs(r) + epsvt
+ if(r.gt.1.0) then
+ l_Pole = .true.
+ h(3,:) = h(3,:) + epsvt
+ call unit_vector_in_3space(h(3,:))
+ endif
+
+!-----------------------------------------------------------------------
+! find the vector that is perpendicular to the local vertical vector
+! and points in the direction of of the North pole, this defines the local
+! north direction. measured in X1
+!-----------------------------------------------------------------------
+ call vector_on_tangent_plane ( h(3,:), g(:), h(2,:) )
+
+!-----------------------------------------------------------------------
+! take the cross product of the local North direction and the local vertical
+! to find the local east vector. still in X1
+!-----------------------------------------------------------------------
+ call cross_product_in_3space ( h(2,:), h(3,:), h(1,:) )
+
+!-----------------------------------------------------------------------
+! put these 3 vectors into a matrix X2
+!-----------------------------------------------------------------------
+ X2(1,:) = h(1,:) ! local east (measured in X1)
+ X2(2,:) = h(2,:) ! local north (measured in X1)
+ X2(3,:) = h(3,:) ! local vertical (measured in X1)
+
+!-----------------------------------------------------------------------
+! compute the transformation matrix
+!-----------------------------------------------------------------------
+ trans_X1_to_X2(:,:) = matmul(X2,transpose(X1))
+
+!-----------------------------------------------------------------------
+! transform (ulon, ulat) into (x,y,z)
+!-----------------------------------------------------------------------
+ p(1) = ux; p(2) = uy; p(3) = uz
+ g(:) = matmul(trans_X1_to_X2(:, :), p(:))
+ ulon = g(1); ulat= g(2);
+
+end subroutine transform_from_xyz_to_lonlat
+
+!======================================================================
+! BEGINNING OF UNIT_VECTOR_IN_3SPACE
+!======================================================================
+ subroutine unit_vector_in_3space (p_1)
+
+!-----------------------------------------------------------------------
+! PURPOSE : normalize p_1 to unit length and overwrite p_1
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(inout)
+!-----------------------------------------------------------------------
+ real , intent(inout) :: &
+ p_1 (:)
+
+!-----------------------------------------------------------------------
+! local
+!-----------------------------------------------------------------------
+ real :: length
+
+ length = SQRT (p_1(1)**2 + p_1(2)**2 + p_1(3)**2 )
+ length = 1.0/length
+ p_1(1) = p_1(1)*length
+ p_1(2) = p_1(2)*length
+ p_1(3) = p_1(3)*length
+
+ end subroutine unit_vector_in_3space
+!======================================================================
+! END OF UNIT_VECTOR_IN_3SPACE
+!======================================================================
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_3SPACE
+!======================================================================
+ subroutine cross_product_in_3space(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real , intent(in) :: &
+ p_1 (:), &
+ p_2 (:)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real , intent(out) :: &
+ p_out (:)
+
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+ end subroutine cross_product_in_3space
+!======================================================================
+! END OF CROSS_PRODUCT_IN_3SPACE
+!======================================================================
+
+!======================================================================
+! BEGINNING OF VECTOR_ON_TANGENT_PLANE
+!======================================================================
+ subroutine vector_on_tangent_plane(p_1, p_2, p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE : given two points measured in (x,y,z) and lying on
+! the unit sphere, find the vector (p_out) that lies on the plane
+! perpendicular to the p_1 vector and points in the direction of
+! the projection of p_2 onto the tangent plane.
+!
+! NOTE : p_1 and p_2 are assumed to be of unit length
+! NOTE : p_out is normalized to unit length
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real , intent(in) :: &
+ p_1 (:), &
+ p_2 (:)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real , intent(out) :: &
+ p_out (:)
+
+!-----------------------------------------------------------------------
+! local
+!-----------------------------------------------------------------------
+ real :: &
+ work (3), t1(3), t2(3)
+
+! work (1) = - p_1(2) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
+! + p_1(3) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) )
+
+! work (2) = + p_1(1) * ( -p_1(2) * p_2(1) + p_1(1) * p_2(2) ) &
+! - p_1(3) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
+
+! work (3) = - p_1(1) * ( p_1(3) * p_2(1) - p_1(1) * p_2(3) ) &
+! + p_1(2) * ( -p_1(3) * p_2(2) + p_1(2) * p_2(3) )
+
+
+ t1(:) = p_2(:) - p_1(:)
+ t2(:) = p_1
+
+ call unit_vector_in_3space (t1)
+ call unit_vector_in_3space (t2)
+
+ call cross_product_in_3space(t1(:), t2(:), work(:))
+ call unit_vector_in_3space (work)
+ call cross_product_in_3space(t2(:),work(:),p_out(:))
+ call unit_vector_in_3space (p_out)
+
+ end subroutine vector_on_tangent_plane
+!======================================================================
+! END OF VECTOR_ON_TANGENT_PLANE
+!======================================================================
+
+end module utilities
</font>
</pre>