<p><b>ringler@lanl.gov</b> 2010-04-06 10:25:23 -0600 (Tue, 06 Apr 2010)</p><p><br>
this directory contains code to convert POP input files to a<br>
MPAS grid.nc file. the code is not particularly clean, but has<br>
been tested using the POP x3 mesh.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/Makefile
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/Makefile (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/Makefile 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,67 @@
+# 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
+NETCDF = /usr/local
+
+# 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 = pop_to_mpas.o \
+ module_cell_indexing.o \
+ module_sphere_utilities.o \
+ module_data_types.o \
+ module_write_netcdf.o
+
+all: ptm
+
+module_sphere_utilities.o: module_data_types.o
+
+pop_to_mpas.o: module_cell_indexing.o module_write_netcdf.o module_sphere_utilities.o module_data_types.o
+
+ptm: $(OBJS)
+ $(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
+
+clean:
+ $(RM) *.o *.mod pop *.f90
+
+.F.o:
+ $(RM) $@ $*.mod
+ $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
+ $(FC) $(FFLAGS) -c $*.f90 $(INCLUDES)
+ #$(RM) $*.f90
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_cell_indexing.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_cell_indexing.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_cell_indexing.F 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,152 @@
+module cell_indexing
+
+! this subroutine provide index mapping for quad meshes dimensioned (nx, ny)
+
+ integer, parameter :: maxEdges = 4
+
+ integer :: nx, ny, nVertLevels, nTracers, nVertexDegree
+ integer, dimension(20) :: nproc
+ character (len = 256) :: POP_GridFileName, POP_KmtFileName, POP_DepthFileName
+
+
+ contains
+
+
+ subroutine cell_indexing_read_nl()
+
+ implicit none
+
+ namelist /periodic_grid/ nx, ny, nVertLevels, nTracers, nproc, nVertexDegree &
+ , POP_GridFileName, POP_KmtFileName, POP_DepthFileName
+
+ nx = 200
+ ny = 200
+ nVertLevels = 1
+ nTracers = 2
+ nproc(:) = -1
+ POP_GridFileName = 'unknown_POP_GridFileName'
+ POP_KmtFileName = 'unknown_POP_KmtFileName'
+ POP_DepthFileName = 'unknown_POP_DepthFileName'
+
+ open(20,file='namelist.input',status='old')
+ read(20,periodic_grid)
+ close(20)
+
+ write(*,periodic_grid)
+
+ end subroutine cell_indexing_read_nl
+
+
+ subroutine cellColRow(idx, iCol, iRow)
+
+ implicit none
+
+ integer, intent(in) :: idx
+ integer, intent(out) :: iCol, iRow
+
+ iRow = ((idx-1) / nx) + 1
+ iCol = mod((idx-1), nx) + 1
+
+ end subroutine cellColRow
+
+
+ integer function cellIdx(iCol, iRow)
+
+ implicit none
+
+ integer, intent(in) :: iCol, iRow
+
+ cellIdx = (iRow-1)*nx + iCol
+
+ end function cellIdx
+
+
+ integer function cellOnCell(iCol, iRow, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCol, iRow, neighborNumber
+
+ integer :: mx, px, my, py
+
+ mx = iCol - 1
+ if (mx == 0) mx = nx
+ my = iRow - 1
+ if (my == 0) my = ny
+ px = iCol + 1
+ if (px == nx + 1) px = 1
+ py = iRow + 1
+ if (py == ny + 1) py = 1
+
+ if (neighborNumber == 1) then
+ cellOnCell = cellIdx(mx, iRow)
+ else if (neighborNumber == 2) then
+ cellOnCell = cellIdx(iCol, my)
+ else if (neighborNumber == 3) then
+ cellOnCell = cellIdx(px, iRow)
+ else if (neighborNumber == 4) then
+ cellOnCell = cellIdx(iCol, py)
+ endif
+
+ end function cellOnCell
+
+
+ integer function edgeOnCell(iCell, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCell, neighborNumber
+
+ integer :: myRow, myCol
+
+ call cellColRow(iCell, myCol, myRow)
+
+ if (neighborNumber == 1) then
+ edgeOnCell = 2*(iCell - 1) + 1
+ else if (neighborNumber == 2) then
+ edgeOnCell = 2*(iCell - 1) + 2
+ else if (neighborNumber == 3) then
+ edgeOnCell = 2*(cellOnCell(myCol, myRow, 3) - 1) + 1
+ else if (neighborNumber == 4) then
+ edgeOnCell = 2*(cellOnCell(myCol, myRow, 4) - 1) + 2
+ end if
+
+ end function edgeOnCell
+
+
+ integer function vertexOnCell(iCell, neighborNumber)
+
+ implicit none
+
+ integer, intent(in) :: iCell, neighborNumber
+
+ integer :: myRow, myCol, mx, my, px, py
+
+ call cellColRow(iCell, myCol, myRow)
+
+!
+! unlike original Todd version, MM put vertex in lower left instead of upper left
+!
+ mx = myCol - 1
+ if (mx == 0) mx = nx
+ my = myRow - 1
+ if (my == 0) my = ny
+ px = myCol + 1
+ if (px == nx + 1) px = 1
+ py = myRow + 1
+ if (py == ny + 1) py = 1
+
+ if (neighborNumber == 1) then
+ vertexOnCell = cellIdx(myCol,myRow)
+ else if (neighborNumber == 2) then
+ vertexOnCell = cellIdx(px,myRow)
+ else if (neighborNumber == 3) then
+ vertexOnCell = cellIdx(px,py)
+ else if (neighborNumber == 4) then
+ vertexOnCell = cellIdx(myCol,py)
+ end if
+
+ end function vertexOnCell
+
+
+end module cell_indexing
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_data_types.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_data_types.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_data_types.F 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,206 @@
+module data_types
+
+ integer, parameter :: LESS = -1, EQUAL = 0, GREATER = 1
+
+ type geo_point
+ real :: lat, lon
+ end type geo_point
+
+ type send_list_ptr
+ integer :: nodeID
+ integer :: nNodeList
+ integer, pointer, dimension(:) :: nodeList
+ type (send_list_ptr), pointer :: next
+ end type send_list_ptr
+
+ type recv_list_ptr
+ integer :: nodeID
+ integer :: nNodeList
+ integer, pointer, dimension(:) :: nodeList
+ type (recv_list_ptr), pointer :: next
+ end type recv_list_ptr
+
+ type adjacency_list
+ integer :: nNodes
+ integer :: nNeighbors
+ integer, pointer, dimension(:) :: neighbor, start, len
+ end type adjacency_list
+
+ type binary_tree
+ integer :: node1, node2
+ integer :: vertex1, vertex2
+ real :: lat1, lon1, lat2, lon2
+ type (binary_tree), pointer :: left, right, parent
+ end type binary_tree
+
+ contains
+
+
+ integer function cmp_points(a, b)
+
+ implicit none
+
+ type (geo_point), intent(in) :: a, b
+
+ if (a%lat > b%lat) then
+ cmp_points = GREATER
+ else if (a%lat == b%lat) then
+ if (a%lon > b%lon) then
+ cmp_points = GREATER
+ else if (a%lon == b%lon) then
+ cmp_points = EQUAL
+ else
+ cmp_points = LESS
+ end if
+ else
+ cmp_points = LESS
+ end if
+
+ end function cmp_points
+
+
+ subroutine swap_points(a, b)
+
+ implicit none
+
+ type (geo_point), intent(inout) :: a, b
+
+ type (geo_point) :: temp
+
+ temp = a
+ a = b
+ b = temp
+
+ end subroutine swap_points
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE CONVERT_ADJACENCY_LIST
+ !
+ ! Convert adjacency list from format provided by STRIPACK to format used in
+ ! our code.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine convert_adjacency_list(n, lend, nvc, list, lptr, alist)
+
+ implicit none
+
+ integer, intent(in) :: n, nvc
+ integer, dimension(n), intent(in) :: lend
+ integer, dimension(nvc), intent(in) :: lptr
+ integer, dimension(nvc), intent(in) :: list
+ type (adjacency_list), intent(inout) :: alist
+
+ integer :: i, j, k, len, ipos
+
+ len = 0
+
+ ! Count total number of nodes
+ do i=1,n
+
+ ! Scan neighbors of i
+ k = lend(i)
+ k = lptr(lend(i))
+ len = len + 1
+
+ do while (k /= lend(i))
+ k = lptr(k)
+ len = len + 1
+ end do
+
+ end do
+
+ alist % nNodes = n
+ alist % nNeighbors = len
+ allocate(alist % neighbor(len))
+ allocate(alist % start(n))
+ allocate(alist % len(n))
+
+ ipos = 0
+ do i=1,n
+
+ ! Scan neighbors of i
+ k = lend(i)
+ k = lptr(lend(i))
+ ipos = ipos + 1
+
+ alist % start(i) = ipos
+ alist % neighbor(ipos) = list(k)
+ alist % len(i) = 1
+
+ do while (k /= lend(i))
+ k = lptr(k)
+ ipos = ipos + 1
+
+ alist % neighbor(ipos) = list(k)
+ alist % len(i) = alist % len(i) + 1
+ end do
+ end do
+
+ end subroutine convert_adjacency_list
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE CONVERT_CORNER_LIST
+ !
+ ! Convert VC list from format provided by STRIPACK to format used in
+ ! our code.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine convert_corner_list(n, lend, nvc, listc, lptr, clist)
+
+ implicit none
+
+ integer, intent(in) :: n, nvc
+ integer, dimension(n), intent(in) :: lend
+ integer, dimension(nvc), intent(in) :: lptr
+ integer, dimension(nvc), intent(in) :: listc
+ type (adjacency_list), intent(inout) :: clist
+
+ integer :: i, j, k, len, ipos
+
+ len = 0
+
+ ! Count total number of nodes
+ do i=1,n
+
+ ! Scan neighbors of i
+ k = lend(i)
+ k = lptr(lend(i))
+ len = len + 1
+
+ do while (k /= lend(i))
+ k = lptr(k)
+ len = len + 1
+ end do
+
+ end do
+
+ clist % nNodes = n
+ clist % nNeighbors = len
+ allocate(clist % neighbor(len))
+ allocate(clist % start(n))
+ allocate(clist % len(n))
+
+ ipos = 0
+ do i=1,n
+
+ ! Scan neighbors of i
+ k = lend(i)
+ k = lptr(lend(i))
+ ipos = ipos + 1
+
+ clist % start(i) = ipos
+ clist % neighbor(ipos) = listc(k)
+ clist % len(i) = 1
+
+ do while (k /= lend(i))
+ k = lptr(k)
+ ipos = ipos + 1
+
+ clist % neighbor(ipos) = listc(k)
+ clist % len(i) = clist % len(i) + 1
+ end do
+ end do
+
+ end subroutine convert_corner_list
+
+end module data_types
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_sphere_utilities.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_sphere_utilities.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_sphere_utilities.F 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,959 @@
+module sphere_utilities
+
+ contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION TRIANGLE_AREA
+!
+! Given the (latitude, longitude) coordinates of the corners of a triangle,
+! plus the radius of the sphere, compute the area of the triangle.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function triangle_area(p1, p2, p3, radius)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2, p3
+ real, intent(in) :: radius
+
+ real :: a, b, c, s, e, pii, tanqe
+
+ pii = 2.*asin(1.0)
+
+ a = sphere_distance(p1,p2,radius)
+ b = sphere_distance(p2,p3,radius)
+ c = sphere_distance(p3,p1,radius)
+ s = 0.5*(a+b+c)
+
+ tanqe = sqrt(tan(0.5*s)*tan(0.5*(s-a))*tan(0.5*(s-b))*tan(0.5*(s-c)))
+ e = 4.*atan(tanqe)
+ triangle_area = radius*radius*e
+
+end function triangle_area
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION OBTUSE
+!
+! Given the (latitude, longitude) coordinates of the corners of a triangle,
+! determine if the triangle is obtuse
+!
+! obtuse.ne.0 then the triangle is obtuse
+! value of 1,2,3 means that angle associated with p1,p2,p3 is > 90
+! obtuse = 0 then the triangle is not obtuse
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+integer function obtuse(p1, p2, p3)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2, p3
+
+ real :: x1(3), x2(3), x3(3), dot, r(3), s(3), rmag, smag
+
+ obtuse = 0
+
+ call convert_lx(x1(1), x1(2), x1(3), 1.0, p1)
+ call convert_lx(x2(1), x2(2), x2(3), 1.0, p2)
+ call convert_lx(x3(1), x3(2), x3(3), 1.0, p3)
+
+ ! test angle formed by x3,x1,x2
+ r(:) = x3(:) - x1(:)
+ s(:) = x2(:) - x1(:)
+ rmag = sqrt(r(1)**2+r(2)**2+r(3)**2)
+ smag = sqrt(s(1)**2+s(2)**2+s(3)**2)
+ r(:) = r(:) / rmag
+ s(:) = s(:) / smag
+ dot = r(1)*s(1) + r(2)*s(2) + r(3)*s(3)
+ if(dot.lt.0) obtuse = 1
+
+ ! test angle formed by x1,x2,x3
+ r(:) = x1(:) - x2(:)
+ s(:) = x3(:) - x2(:)
+ rmag = sqrt(r(1)**2+r(2)**2+r(3)**2)
+ smag = sqrt(s(1)**2+s(2)**2+s(3)**2)
+ r(:) = r(:) / rmag
+ s(:) = s(:) / smag
+ dot = r(1)*s(1) + r(2)*s(2) + r(3)*s(3)
+ if(dot.lt.0) obtuse = 2
+
+ ! test angle formed by x2,x3,x1
+ r(:) = x2(:) - x3(:)
+ s(:) = x1(:) - x3(:)
+ rmag = sqrt(r(1)**2+r(2)**2+r(3)**2)
+ smag = sqrt(s(1)**2+s(2)**2+s(3)**2)
+ r(:) = r(:) / rmag
+ s(:) = s(:) / smag
+ dot = r(1)*s(1) + r(2)*s(2) + r(3)*s(3)
+ if(dot.lt.0) obtuse = 3
+
+end function obtuse
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION SPHERE_DISTANCE
+!
+! Given two (latitude, longitude) coordinates on the surface of a sphere,
+! plus the radius of the sphere, compute the great circle distance between
+! the points.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function sphere_distance(p1, p2, radius)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2
+ real, intent(in) :: radius
+
+ real :: arg1
+
+ arg1 = sqrt( sin(0.5*(p2%lat-p1%lat))**2 + &
+ cos(p1%lat)*cos(p2%lat)*sin(0.5*(p2%lon-p1%lon))**2 )
+ sphere_distance = 2.*radius*asin(arg1)
+
+end function sphere_distance
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION PLANE_DISTANCE
+!
+! Given two (latitude, longitude) coordinates on the surface of a sphere,
+! plus the radius of the sphere, compute the secant distance between
+! the points.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function plane_distance(p1, p2, radius)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2
+ real, intent(in) :: radius
+
+ real :: x1, x2, y1, y2, z1, z2
+
+ z1 = sin(p1%lat)
+ z2 = sin(p2%lat)
+ x1 = cos(p1%lon)*cos(p1%lat)
+ x2 = cos(p2%lon)*cos(p2%lat)
+ y1 = sin(p1%lon)*cos(p1%lat)
+ y2 = sin(p2%lon)*cos(p2%lat)
+
+ plane_distance = radius*sqrt((z1-z2)**2+(x1-x2)**2+(y1-y2)**2)
+
+end function plane_distance
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION ARC_ANGLE
+!
+! Given two (latitude, longitude) coordinates on the surface of a sphere,
+! compute the angle between the points as measured from the origin of the
+! sphere.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function arc_angle(p1, p2)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2
+
+ real :: arg1
+
+ arg1 = sqrt( sin(0.5*(p2%lat-p1%lat))**2 + &
+ cos(p1%lat)*cos(p2%lat)*sin(0.5*(p2%lon-p1%lon))**2 )
+ arc_angle = 2.*asin(arg1)
+
+end function arc_angle
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION GREAT_CIRCLE_POINTS
+!
+! Return n points equally spaced along the great circle arc between (lat1,lon1)
+! and (lat2,lon2). These points include the end points of the arc.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine great_circle_points(p1, p2, pl, n)
+
+ use data_types
+
+ implicit none
+
+ integer, intent(in) :: n
+ type (geo_point), intent(in) :: p1, p2
+ type (geo_point), dimension(n), intent(inout) :: pl
+
+ real :: x1, x2, y1, y2, z1, z2
+ real :: dx, dl
+ real :: x, y, z
+ integer :: i
+ real :: dtheta, dinc, dt
+
+ real :: pii, rtod
+
+ pii = 2.*asin(1.0)
+ rtod = 180./pii
+
+! write(6,*) ' in gcp ',rtod*lat1,rtod*lon1,rtod*lat2,rtod*lon2
+
+ if (n < 2) then
+ write(6,*) ' n less than 2 in great_circle_points '
+ stop
+ end if
+
+ if (n == 2) then
+ pl(1) = p1
+ pl(2) = p2
+ end if
+
+ dtheta = arc_angle(p1, p2)
+ dinc = dtheta/float(n-1)
+
+ call convert_lx(x1,y1,z1,1.,p1)
+ call convert_lx(x2,y2,z2,1.,p2)
+
+! set the end points
+
+ pl(1) = p1
+ pl(n) = p2
+
+! write(6,*) ' x1,y1,z1 ',x1,y1,z1
+! write(6,*) ' x2,y2,z2 ',x2,y2,z2
+
+! compute the interior points. see notes for derivation
+
+ do i=2,n-1
+ dt = float(i-1)*dinc
+
+ if (dt <= 0.5*dtheta) then
+ dx = 1.-tan(0.5*dtheta-dt)/tan(0.5*dtheta)
+! write(6,*) ' case 1 ',dx
+ x = x1+0.5*dx*(x2-x1)
+ y = y1+0.5*dx*(y2-y1)
+ z = z1+0.5*dx*(z2-z1)
+ else
+ dt = dtheta-dt
+ dx = 1.-tan(0.5*dtheta-dt)/tan(0.5*dtheta)
+! write(6,*) ' case 2 ',dx
+ x = x2+0.5*dx*(x1-x2)
+ y = y2+0.5*dx*(y1-y2)
+ z = z2+0.5*dx*(z1-z2)
+ end if
+
+! write(6,*) ' x,y,z ',x,y,z
+
+ call convert_xl(x,y,z,pl(i))
+ enddo
+
+end subroutine great_circle_points
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE DIVIDE_TRIANGLE
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!subroutine divide_triangle( p1, p2, p3, pnew)
+!
+! use data_types
+!
+! implicit none
+!
+! type (geo_point), intent(in) :: p1, p2, p3
+! type (geo_point), dimension(6), intent(inout) :: pnew
+!
+! real :: t_area, area_total, radius
+! type (geo_point), dimension(3) :: pts
+! type (geo_point) :: c
+!
+! radius = 1.
+! pnew(1) = p1
+! pnew(4) = p2
+! pnew(6) = p3
+!
+! call great_circle_points(p1,p2,pts,3)
+! pnew(2) = pts(2)
+!
+! call great_circle_points(p1,p3,pts,3)
+! pnew(3) = pts(2)
+!
+! call great_circle_points(p2,p3,pts,3)
+! pnew(5) = pts(2)
+!
+!
+! write(6,*) ' '
+! write(6,*) ' original triangle '
+! write(6,*) p1%lat, p1%lon
+! write(6,*) p2%lat, p2%lon
+! write(6,*) p3%lat, p3%lon
+!
+! t_area = triangle_area(p1,p2,p3,radius)
+! write(6,*) ' area ',t_area
+! call compute_voronoi_corner(p1,p2,p3,c)
+! write(6,*) ' voronoi corner ',c%lat,c%lon
+!
+! area_total = 0.
+!
+! write(6,*) ' '
+! write(6,*) ' new triangles '
+!
+! write(6,*) ' triangle 1 '
+! write(6,*) pnew(1)%lat,pnew(1)%lon
+! write(6,*) pnew(1)%lat,pnew(2)%lon
+! write(6,*) pnew(1)%lat,pnew(3)%lon
+! t_area = triangle_area( pnew(1),pnew(2),pnew(3),radius)
+! area_total = area_total + t_area
+! write(6,*) ' area ',t_area
+!
+! write(6,*) ' triangle 2 '
+! write(6,*) pnew(2)%lat,pnew(2)%lon
+! write(6,*) pnew(4)%lat,pnew(4)%lon
+! write(6,*) pnew(5)%lat,pnew(5)%lon
+! t_area = triangle_area( pnew(2),pnew(4),pnew(5),radius)
+! area_total = area_total + t_area
+! write(6,*) ' area ',t_area
+!
+! write(6,*) ' triangle 3 '
+! write(6,*) pnew(2)%lat,pnew(2)%lon
+! write(6,*) pnew(5)%lat,pnew(5)%lon
+! write(6,*) pnew(3)%lat,pnew(3)%lon
+! t_area = triangle_area( pnew(2),pnew(5),pnew(3),radius)
+! area_total = area_total + t_area
+! write(6,*) ' area ',t_area
+!
+! write(6,*) ' triangle 4 '
+! write(6,*) pnew(3)%lat,pnew(3)%lon
+! write(6,*) pnew(5)%lat,pnew(5)%lon
+! write(6,*) pnew(6)%lat,pnew(6)%lon
+! t_area = triangle_area( pnew(3),pnew(5),pnew(6),radius)
+! area_total = area_total + t_area
+! write(6,*) ' area ',t_area
+! write(6,*) ' total area is ',area_total
+!
+!end subroutine divide_triangle
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE COMPUTE_VORONOI_CORNER
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!subroutine compute_voronoi_corner( p0, p1, p2, vc )
+!
+! use data_types
+!
+! implicit none
+!
+! type (geo_point), intent(in) :: p0, p1, p2
+! type (geo_point), intent(out) :: vc
+!
+! real :: x0, y0, z0, x1, y1, z1, x2, y2, z2, xc, yc, zc, cabs
+! real :: a1, a2, a3, b1, b2, b3
+! real :: dot0
+!
+! z0 = sin(p0%lat)
+! z1 = sin(p1%lat)
+! z2 = sin(p2%lat)
+!
+! x0 = cos(p0%lon)*cos(p0%lat)
+! x1 = cos(p1%lon)*cos(p1%lat)
+! x2 = cos(p2%lon)*cos(p2%lat)
+!
+! y0 = sin(p0%lon)*cos(p0%lat)
+! y1 = sin(p1%lon)*cos(p1%lat)
+! y2 = sin(p2%lon)*cos(p2%lat)
+!
+! a1 = x2-x0
+! a2 = y2-y0
+! a3 = z2-z0
+!
+! b1 = x1-x0
+! b2 = y1-y0
+! b3 = z1-z0
+!
+!
+! xc = a2*b3-a3*b2
+! yc = a3*b1-a1*b3
+! zc = a1*b2-a2*b1
+! cabs = sqrt(xc*xc+yc*yc+zc*zc)
+!
+!! write(6,*) ' cabs = ',cabs
+!! write(6,*) ' xc, yc, zc = ',xc,yc,zc
+!! write(6,*) ' x0, y0, z0 = ',x0,y0,z0
+!! write(6,*) ' x1, y1, z1 = ',x1,y1,z1
+!! write(6,*) ' x2, y2, z2 = ',x2,y2,z2
+! dot0 = x0*xc+y0*yc+z0*zc
+!! write(6,*) ' dot is ',dot0
+!
+! if( dot0 < 0.) then ! flip p1 with p2
+!
+! z2 = sin(p1%lat)
+! z1 = sin(p2%lat)
+!
+! x2 = cos(p1%lon)*cos(p1%lat)
+! x1 = cos(p2%lon)*cos(p2%lat)
+!
+! y2 = sin(p1%lon)*cos(p1%lat)
+! y1 = sin(p2%lon)*cos(p2%lat)
+!
+! a1 = x2-x0
+! a2 = y2-y0
+! a3 = z2-z0
+!
+! b1 = x1-x0
+! b2 = y1-y0
+! b3 = z1-z0
+!
+!
+! xc = a2*b3-a3*b2
+! yc = a3*b1-a1*b3
+! zc = a1*b2-a2*b1
+! cabs = sqrt(xc*xc+yc*yc+zc*zc)
+!
+!! write(6,*) ' flipping '
+!! write(6,*) ' cabs = ',cabs
+!! write(6,*) ' xc, yc, zc = ',xc,yc,zc
+! dot0 = x0*xc+y0*yc+z0*zc
+!! write(6,*) ' dot is ',dot0
+!
+! end if
+!
+!
+! xc = xc/cabs
+! yc = yc/cabs
+! zc = zc/cabs
+!
+! call convert_xl(xc,yc,zc,vc)
+!
+!end subroutine compute_voronoi_corner
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! 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, latlon)
+
+ use data_types
+
+ implicit none
+
+ real, intent(in) :: radius
+ type (geo_point), intent(in) :: latlon
+ real, intent(out) :: x, y, z
+
+ z = radius * sin(latlon%lat)
+ x = radius * cos(latlon%lon) * cos(latlon%lat)
+ y = radius * sin(latlon%lon) * cos(latlon%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, latlon)
+
+ use data_types
+
+ implicit none
+
+ real, intent(in) :: x, y, z
+ type (geo_point), intent(out) :: latlon
+
+ 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)
+
+ latlon%lat = asin(z/dl)
+
+! check for being close to either pole
+
+ if (abs(x) > eps) then
+
+ if (abs(y) > eps) then
+
+ latlon%lon = atan(abs(y/x))
+
+ if ((x <= 0.) .and. (y >= 0.)) then
+ latlon%lon = pii-latlon%lon
+ else if ((x <= 0.) .and. (y < 0.)) then
+ latlon%lon = latlon%lon+pii
+ else if ((x >= 0.) .and. (y <= 0.)) then
+ latlon%lon = 2*pii-latlon%lon
+ end if
+
+ else ! we're either on longitude 0 or 180
+
+ if (x > 0) then
+ latlon%lon = 0.
+ else
+ latlon%lon = pii
+ end if
+
+ end if
+
+ else if (abs(y) > eps) then
+
+ if (y > 0) then
+ latlon%lon = pii/2.
+ else
+ latlon%lon = 3.*pii/2.
+ end if
+
+ else ! we are at a pole
+
+ latlon%lon = 0.
+
+ end if
+
+end subroutine convert_xl
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE GC_INTERSECT
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine gc_intersect(p0, p1, p2, p3, pc)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p0, p1, p2, p3
+ type (geo_point), intent(out) :: pc
+
+ real :: x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3
+ real :: n1, n2, n3, m1, m2, m3
+ real :: xc, yc, zc, dot
+ real, parameter :: radius=1.0
+
+ call convert_lx(x0,y0,z0,radius,p0)
+ call convert_lx(x1,y1,z1,radius,p1)
+ call convert_lx(x2,y2,z2,radius,p2)
+ call convert_lx(x3,y3,z3,radius,p3)
+
+ n1 = (y0 * z1 - y1 * z0)
+ n2 = -(x0 * z1 - x1 * z0)
+ n3 = (x0 * y1 - x1 * y0)
+
+ m1 = (y2 * z3 - y3 * z2)
+ m2 = -(x2 * z3 - x3 * z2)
+ m3 = (x2 * y3 - x3 * y2)
+
+ xc = (n2 * m3 - n3 * m2)
+ yc = -(n1 * m3 - n3 * m1)
+ zc = (n1 * m2 - n2 * m1)
+
+ dot = x0*xc + y0*yc + z0*zc
+
+ if (dot < 0.0) then
+ xc = -xc
+ yc = -yc
+ zc = -zc
+ end if
+
+ call convert_xl(xc,yc,zc,pc)
+
+end subroutine gc_intersect
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION POS_ANG
+!
+! Normalize an angle, given in radians, to lie in the interval [0,2*PI].
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function pos_ang(angle)
+
+ implicit none
+
+ real, intent(in) :: angle
+
+ real :: pii
+
+ pii = 2.*asin(1.0)
+ pos_ang = angle
+
+ if(angle > 2.*pii) then
+ pos_ang = angle - 2.*pii
+ else if(angle < 0.) then
+ pos_ang = angle + 2.*pii
+ end if
+
+end function pos_ang
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION MERIDIAN_ANGLE
+!
+! Find the angle between the meridian that intersects point (lat1,lon1)
+! and the great circle passing through points (lat1,lon1) (lat2,lon2).
+! (lat1,lon1) is the vertex of the angle.
+!
+! Convention: zero points north, 90 points west, -90 point east,
+! points south 180, -180
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function meridian_angle(p1, p2)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2
+
+type (geo_point) :: np
+
+ real :: pii, da, db, dc
+ type (geo_point) :: p3
+ real :: cosa
+ real :: eps
+ parameter (eps = 1.e-04)
+real :: ax, ay, az
+real :: bx, by, bz
+real :: cx, cy, cz
+
+np = p1
+np%lat = np%lat + 0.05
+
+call convert_lx(ax, ay, az, 1.0, p1)
+call convert_lx(bx, by, bz, 1.0, np)
+call convert_lx(cx, cy, cz, 1.0, p2)
+
+meridian_angle = plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, ax, ay, az)
+return
+
+ if (p1%lon == p2%lon) then
+
+ meridian_angle = 0.0
+
+ else
+
+ pii = 2.*asin(1.0)
+ dc = arc_angle(p1,p2)
+
+ p3%lon = p1%lon
+ if (p1%lat + dc <= pii/2.0) then
+ p3%lat = p1%lat+dc
+ else
+ p3%lat = p1%lat-dc
+ end if
+ db = arc_angle(p1,p3)
+ da = arc_angle(p2,p3)
+
+! see spherical trig section on online wolfram pages - eq(11) ->
+
+ cosa = max(-1.,min(1.,(cos(da)-cos(db)*cos(dc))/(sin(db)*sin(dc))))
+ meridian_angle = acos(cosa)
+
+
+ if (((p2%lon > p1%lon) .and. (p2%lon - p1%lon <= pii)) .or. &
+ ((p2%lon < p1%lon) .and. (p1%lon - p2%lon >= pii))) then
+ meridian_angle = -abs(meridian_angle)
+ else
+ meridian_angle = abs(meridian_angle)
+ end if
+
+ end if
+
+end function meridian_angle
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE CENTER_OF_MASS
+!
+! Find centriod of the triangle whose corners are at (lat1,lon1), (lat2,lon2),
+! and (lat3,lon3).
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine center_of_mass(p1, p2, p3, pc)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2, p3
+ type (geo_point), intent(out) :: pc
+
+ real :: x1, x2, x3, xc
+ real :: y1, y2, y3, yc
+ real :: z1, z2, z3, zc
+
+ call convert_lx(x1,y1,z1,1.,p1)
+ call convert_lx(x2,y2,z2,1.,p2)
+ call convert_lx(x3,y3,z3,1.,p3)
+
+ xc = (x1+x2+x3)/3.
+ yc = (y1+y2+y3)/3.
+ zc = (z1+z2+z3)/3.
+
+ call convert_xl(xc,yc,zc,pc)
+
+end subroutine center_of_mass
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE DIVIDE_TRIANGLE
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine divide_triangle(p1, p2, p3, n, p)
+
+ use data_types
+
+ implicit none
+
+ type (geo_point), intent(in) :: p1, p2, p3
+ integer, intent(in) :: n
+ type (geo_point), dimension(3,n), intent(out) :: p
+
+ integer :: i, j, k
+ integer :: glevel ! Level of decomposition
+ type (geo_point), allocatable, dimension(:) :: p1p2, p1p3
+ type (geo_point), allocatable, dimension(:,:) :: line
+
+ glevel = nint(log(real(n)) / log(4.0)) ! Each subdivision gives four times the number of
+ ! triangles, so log4(n) gives the level decomposition
+
+ glevel = (2 ** glevel) + 1
+ allocate(line(glevel, glevel))
+ allocate(p1p2(glevel))
+ allocate(p1p3(glevel))
+
+ call great_circle_points(p1, p2, p1p2, glevel)
+ call great_circle_points(p1, p3, p1p3, glevel)
+
+ line(1,1) = p1
+ line(1,2) = p1p2(2)
+ line(2,2) = p1p3(2)
+
+ do i = 3,glevel
+ call great_circle_points(p1p2(i), p1p3(i), line(:,i), i)
+!do j=1,i
+!write(0,*) j,i,' P ',line(j,i)%lat*180./3.14159, line(j,i)%lon*180./3.14159
+!end do
+ end do
+
+ k = 1
+ do i = 1,glevel-1
+ do j = 1,i
+ p(1,k) = line(j,i)
+ p(2,k) = line(j,i+1)
+ p(3,k) = line(j+1,i+1)
+!write(0,*) j,i, ' - ',p(1,k)%lat*180./3.14159,p(1,k)%lon*180./3.14159
+!write(0,*) j,i+1, ' - ',p(2,k)%lat*180./3.14159,p(2,k)%lon*180./3.14159
+!write(0,*) j+1,i+1, ' - ',p(3,k)%lat*180./3.14159,p(3,k)%lon*180./3.14159
+ k = k + 1
+ end do
+ end do
+
+!write(0,*) '-----------'
+ do i = glevel,3,-1
+ do j = 2,i-1
+ p(1,k) = line(j,i)
+ p(2,k) = line(j,i-1)
+ p(3,k) = line(j-1,i-1)
+!write(0,*) j,i, ' - ',p(1,k)%lat*180./3.14159,p(1,k)%lon*180./3.14159
+!write(0,*) j,i-1, ' - ',p(2,k)%lat*180./3.14159,p(2,k)%lon*180./3.14159
+!write(0,*) j-1,i-1, ' - ',p(3,k)%lat*180./3.14159,p(3,k)%lon*180./3.14159
+ k = k + 1
+ end do
+ end do
+
+ deallocate(line)
+ deallocate(p1p2)
+ deallocate(p1p3)
+
+end subroutine divide_triangle
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE POINT_TO_PLANE
+!
+! Find projection (xp, yp, zp) of a point (Qx,Qy,Qz) onto the plane defined by
+! the equation ax+by+cz+d=0
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine point_to_plane(a, b, c, d, Qx, Qy, Qz, xp, yp, zp)
+
+ implicit none
+
+ real, intent(in) :: a, b, c, d ! The coefficients in the equation of the plane
+ real, intent(in) :: Qx, Qy, Qz ! The coordinates of the point Q to be projected to the plane
+ real, intent(out) :: xp, yp, zp ! The coordinates of the point projected in the plane
+
+ real :: Px, Py, Pz ! A point P in the plane ax + by + cz + d = 0
+ real :: PQx, PQy, PQz ! Components of the vector from P to Q
+ real :: PQn ! The dot product of PQ and the vector normal to the plane
+ real :: m2 ! The magnitude and squared magnitude of the vector n normal to the plane
+
+ m2 = (a**2.0 + b**2.0 + c**2.0)
+
+ Px = -d*a/m2
+ Py = -d*b/m2
+ Pz = -d*c/m2
+
+ PQx = Qx - Px
+ PQy = Qy - Py
+ PQz = Qz - Pz
+
+ PQn = PQx * a + PQy * b + PQz * c
+
+ ! . Q
+ ! n ^ /
+ ! | /
+ ! |/
+ ! ----------.-------------------
+ ! P
+
+ xp = Qx - PQn * a / m2
+ yp = Qy - PQn * b / m2
+ zp = Qz - PQn * c / m2
+
+end subroutine point_to_plane
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE POINT_TO_SPHERE
+!
+! Find projection (xp, yp, zp) of a point (Qx,Qy,Qz) in the plane defined by
+! the equation ax+by+cz+d=0 onto the surface of the sphere with radius r
+! centered at the origin.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine point_to_sphere(a, b, c, d, r, Qx, Qy, Qz, xp, yp, zp)
+
+ implicit none
+
+ real, intent(in) :: a, b, c, d ! The coefficients in the equation of the plane
+ real, intent(in) :: r ! The radius of the sphere
+ real, intent(in) :: Qx, Qy, Qz ! The coordinates of the point Q to be projected to the sphere
+ real, intent(out) :: xp, yp, zp ! The coordinates of the point projected to the sphere
+
+ real :: aa, bb, cc ! Coefficients of quadratic equation
+ real :: disc, t1, t2
+
+ ! Solve for the interesection of the line (Qx - at, Qy - bt, Qz - ct) and the
+ ! sphere x^2 + y^2 + z^2 - r^2 = 0
+ aa = a**2.0 + b**2.0 + c**2.0
+ bb = -2.0*(Qx*a + Qy*b + Qz*c)
+ cc = Qx**2.0 + Qy**2.0 + Qz**2.0 - r**2.0
+
+ disc = bb**2.0 - 4.0*aa*cc
+
+ if (disc < 0.0) then ! Point has no projection on the surface of the sphere
+ xp = 0.0
+ yp = 0.0
+ zp = 0.0
+ else if (disc == 0.0) then ! Point has exactly one projection (line through point and
+ t1 = -bb / (2.0*aa)
+ xp = Qx - a*t1 ! and normal to plane is tangent to sphere
+ yp = Qy - b*t1
+ zp = Qz - c*t1
+ else ! Point has two projections; choose the one that is closest
+ t1 = (-bb + sqrt(disc)) / (2.0*aa)
+ t2 = (-bb - sqrt(disc)) / (2.0*aa)
+ if (abs(t1) <= abs(t2)) then
+ xp = Qx - a*t1
+ yp = Qy - b*t1
+ zp = Qz - c*t1
+ else
+ xp = Qx - a*t2
+ yp = Qy - b*t2
+ zp = Qz - c*t2
+ end if
+ end if
+
+end subroutine point_to_sphere
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! SUBROUTINE ROTATE_ABOUT_VECTOR
+!
+! Rotates the point (x,y,z) through an angle theta about the vector
+! originating at (a, b, c) and having direction (u, v, w).
+!
+! Reference: http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/ArbitraryAxisRotation.html
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine rotate_about_vector(x, y, z, theta, a, b, c, u, v, w, xp, yp, zp)
+
+ implicit none
+
+ real, intent(in) :: x, y, z, theta, a, b, c, u, v, w
+ real, intent(out) :: xp, yp, zp
+
+ real :: vw2, uw2, uv2
+ real :: m
+
+ vw2 = v**2.0 + w**2.0
+ uw2 = u**2.0 + w**2.0
+ uv2 = u**2.0 + v**2.0
+ m = sqrt(u**2.0 + v**2.0 + w**2.0)
+
+ xp = (a*vw2 + u*(-b*v-c*w+u*x+v*y+w*z) + ((x-a)*vw2+u*(b*v+c*w-v*y-w*z))*cos(theta) + m*(-c*v+b*w-w*y+v*z)*sin(theta))/m**2.0
+ yp = (b*uw2 + v*(-a*u-c*w+u*x+v*y+w*z) + ((y-b)*uw2+v*(a*u+c*w-u*x-w*z))*cos(theta) + m*( c*u-a*w+w*x-u*z)*sin(theta))/m**2.0
+ zp = (c*uv2 + w*(-a*u-b*v+u*x+v*y+w*z) + ((z-c)*uv2+w*(a*u+b*v-u*x-v*y))*cos(theta) + m*(-b*u+a*v-v*x+u*y)*sin(theta))/m**2.0
+
+end subroutine rotate_about_vector
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FUNCTION PLANE_ANGLE
+!
+! Computes the angle between vectors AB and AC, given points A, B, and C, and
+! a vector (u,v,w) normal to the plane.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+
+ implicit none
+
+ real, intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w
+
+ real :: ABx, ABy, ABz ! The components of the vector AB
+ real :: mAB ! The magnitude of AB
+ real :: ACx, ACy, ACz ! The components of the vector AC
+ real :: mAC ! The magnitude of AC
+
+ real :: Dx ! The i-components of the cross product AB x AC
+ real :: Dy ! The j-components of the cross product AB x AC
+ real :: Dz ! The k-components of the cross product AB x AC
+
+ real :: cos_angle
+
+ ABx = bx - ax
+ ABy = by - ay
+ ABz = bz - az
+ mAB = sqrt(ABx**2.0 + ABy**2.0 + ABz**2.0)
+
+ ACx = cx - ax
+ ACy = cy - ay
+ ACz = cz - az
+ mAC = sqrt(ACx**2.0 + ACy**2.0 + ACz**2.0)
+
+
+ Dx = (ABy * ACz) - (ABz * ACy)
+ Dy = -((ABx * ACz) - (ABz * ACx))
+ Dz = (ABx * ACy) - (ABy * ACx)
+
+ cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
+
+ if (cos_angle < -1.0) then
+ cos_angle = -1.0
+ else if (cos_angle > 1.0) then
+ cos_angle = 1.0
+ end if
+
+ if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
+ plane_angle = acos(cos_angle)
+ else
+ plane_angle = -acos(cos_angle)
+ end if
+
+end function plane_angle
+
+end module sphere_utilities
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_write_netcdf.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_write_netcdf.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/module_write_netcdf.F 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,631 @@
+module write_netcdf
+
+ integer :: wr_ncid
+ integer :: wrDimIDTime
+ integer :: wrDimIDnCells
+ integer :: wrDimIDnEdges
+ integer :: wrDimIDnVertices
+ integer :: wrDimIDmaxEdges
+ integer :: wrDimIDmaxEdges2
+ integer :: wrDimIDTWO
+ integer :: wrDimIDvertexDegree
+ integer :: wrDimIDnVertLevels
+ integer :: wrDimIDnTracers
+ integer :: wrVarIDlatCell
+ integer :: wrVarIDlonCell
+ integer :: wrVarIDxCell
+ integer :: wrVarIDyCell
+ integer :: wrVarIDzCell
+ integer :: wrVarIDindexToCellID
+ integer :: wrVarIDlatEdge
+ integer :: wrVarIDlonEdge
+ integer :: wrVarIDxEdge
+ integer :: wrVarIDyEdge
+ integer :: wrVarIDzEdge
+ integer :: wrVarIDindexToEdgeID
+ integer :: wrVarIDlatVertex
+ integer :: wrVarIDlonVertex
+ integer :: wrVarIDxVertex
+ integer :: wrVarIDyVertex
+ integer :: wrVarIDzVertex
+ integer :: wrVarIDindexToVertexID
+ integer :: wrVarIDcellsOnEdge
+ integer :: wrVarIDnEdgesOnCell
+ integer :: wrVarIDnEdgesOnEdge
+ integer :: wrVarIDedgesOnCell
+ integer :: wrVarIDedgesOnEdge
+ integer :: wrVarIDweightsOnEdge
+ integer :: wrVarIDdvEdge
+ integer :: wrVarIDdcEdge
+ integer :: wrVarIDangleEdge
+ integer :: wrVarIDareaCell
+ integer :: wrVarIDareaTriangle
+ integer :: wrVarIDcellsOnCell
+ integer :: wrVarIDverticesOnCell
+ integer :: wrVarIDverticesOnEdge
+ integer :: wrVarIDedgesOnVertex
+ integer :: wrVarIDcellsOnVertex
+ integer :: wrVarIDkiteAreasOnVertex
+ integer :: wrVarIDfEdge
+ integer :: wrVarIDfVertex
+ integer :: wrVarIDh_s
+ integer :: wrVarIDu
+ integer :: wrVarIDboundaryEdge
+ integer :: wrVarIDboundaryVertex
+ integer :: wrVarIDu_src
+ integer :: wrVarIDv
+ integer :: wrVarIDh
+ integer :: wrVarIDvh
+ integer :: wrVarIDcirculation
+ integer :: wrVarIDvorticity
+ integer :: wrVarIDke
+ integer :: wrVarIDtracers
+
+ integer :: wrLocalnCells
+ integer :: wrLocalnEdges
+ integer :: wrLocalnVertices
+ integer :: wrLocalmaxEdges
+ integer :: wrLocalnVertLevels
+ integer :: wrLocalnTracers
+ integer :: wrLocalnVertexDegree
+
+ contains
+
+ subroutine write_netcdf_init( &
+ nCells, &
+ nEdges, &
+ nVertices, &
+ maxEdges, &
+ nVertLevels, &
+ nTracers, &
+ nVertexDegree &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertices
+ integer, intent(in) :: maxEdges
+ integer, intent(in) :: nVertLevels
+ integer, intent(in) :: nTracers
+ integer, intent(in) :: nVertexDegree
+
+ integer :: nferr
+ integer, dimension(10) :: dimlist
+
+
+ wrLocalnCells = nCells
+ wrLocalnEdges = nEdges
+ wrLocalnVertices = nVertices
+ wrLocalmaxEdges = maxEdges
+ wrLocalnVertLevels = nVertLevels
+ wrLocalnTracers = nTracers
+ wrLocalnVertexDegree = nVertexDegree
+
+ nferr = nf_create('grid.nc', IOR(NF_CLOBBER,NF_64BIT_OFFSET), wr_ncid)
+
+ !
+ ! Define dimensions
+ !
+ nferr = nf_def_dim(wr_ncid, 'nCells', nCells, wrDimIDnCells)
+ nferr = nf_def_dim(wr_ncid, 'nEdges', nEdges, wrDimIDnEdges)
+ nferr = nf_def_dim(wr_ncid, 'nVertices', nVertices, wrDimIDnVertices)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges', maxEdges, wrDimIDmaxEdges)
+ nferr = nf_def_dim(wr_ncid, 'maxEdges2', 2*maxEdges, wrDimIDmaxEdges2)
+ nferr = nf_def_dim(wr_ncid, 'TWO', 2, wrDimIDTWO)
+! nferr = nf_def_dim(wr_ncid, 'THREE', 3, wrDimIDTHREE)
+ nferr = nf_def_dim(wr_ncid, 'vertexDegree', nVertexDegree, wrDimIDvertexDegree)
+ nferr = nf_def_dim(wr_ncid, 'nVertLevels', nVertLevels, wrDimIDnVertLevels)
+ nferr = nf_def_dim(wr_ncid, 'nTracers', nTracers, wrDimIDnTracers)
+ nferr = nf_def_dim(wr_ncid, 'Time', NF_UNLIMITED, wrDimIDTime)
+
+ !
+ ! Define variables
+ !
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'latCell', NF_DOUBLE, 1, dimlist, wrVarIDlatCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'lonCell', NF_DOUBLE, 1, dimlist, wrVarIDlonCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'xCell', NF_DOUBLE, 1, dimlist, wrVarIDxCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'yCell', NF_DOUBLE, 1, dimlist, wrVarIDyCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'zCell', NF_DOUBLE, 1, dimlist, wrVarIDzCell)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'indexToCellID', NF_INT, 1, dimlist, wrVarIDindexToCellID)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'latEdge', NF_DOUBLE, 1, dimlist, wrVarIDlatEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'lonEdge', NF_DOUBLE, 1, dimlist, wrVarIDlonEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'xEdge', NF_DOUBLE, 1, dimlist, wrVarIDxEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'yEdge', NF_DOUBLE, 1, dimlist, wrVarIDyEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'zEdge', NF_DOUBLE, 1, dimlist, wrVarIDzEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'indexToEdgeID', NF_INT, 1, dimlist, wrVarIDindexToEdgeID)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'latVertex', NF_DOUBLE, 1, dimlist, wrVarIDlatVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'lonVertex', NF_DOUBLE, 1, dimlist, wrVarIDlonVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'xVertex', NF_DOUBLE, 1, dimlist, wrVarIDxVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'yVertex', NF_DOUBLE, 1, dimlist, wrVarIDyVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'zVertex', NF_DOUBLE, 1, dimlist, wrVarIDzVertex)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'indexToVertexID', NF_INT, 1, dimlist, wrVarIDindexToVertexID)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'cellsOnEdge', NF_INT, 2, dimlist, wrVarIDcellsOnEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnCell', NF_INT, 1, dimlist, wrVarIDnEdgesOnCell)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'nEdgesOnEdge', NF_INT, 1, dimlist, wrVarIDnEdgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'edgesOnCell', NF_INT, 2, dimlist, wrVarIDedgesOnCell)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'edgesOnEdge', NF_INT, 2, dimlist, wrVarIDedgesOnEdge)
+ dimlist( 1) = wrDimIDmaxEdges2
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'weightsOnEdge', NF_DOUBLE, 2, dimlist, wrVarIDweightsOnEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dvEdge', NF_DOUBLE, 1, dimlist, wrVarIDdvEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'dcEdge', NF_DOUBLE, 1, dimlist, wrVarIDdcEdge)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'angleEdge', NF_DOUBLE, 1, dimlist, wrVarIDangleEdge)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'areaCell', NF_DOUBLE, 1, dimlist, wrVarIDareaCell)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'areaTriangle', NF_DOUBLE, 1, dimlist, wrVarIDareaTriangle)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'cellsOnCell', NF_INT, 2, dimlist, wrVarIDcellsOnCell)
+ dimlist( 1) = wrDimIDmaxEdges
+ dimlist( 2) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'verticesOnCell', NF_INT, 2, dimlist, wrVarIDverticesOnCell)
+ dimlist( 1) = wrDimIDTWO
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'verticesOnEdge', NF_INT, 2, dimlist, wrVarIDverticesOnEdge)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'edgesOnVertex', NF_INT, 2, dimlist, wrVarIDedgesOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'cellsOnVertex', NF_INT, 2, dimlist, wrVarIDcellsOnVertex)
+ dimlist( 1) = wrDimIDvertexDegree
+ dimlist( 2) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'kiteAreasOnVertex', NF_DOUBLE, 2, dimlist, wrVarIDkiteAreasOnVertex)
+ dimlist( 1) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, 'fEdge', NF_DOUBLE, 1, dimlist, wrVarIDfEdge)
+ dimlist( 1) = wrDimIDnVertices
+ nferr = nf_def_var(wr_ncid, 'fVertex', NF_DOUBLE, 1, dimlist, wrVarIDfVertex)
+ dimlist( 1) = wrDimIDnCells
+ nferr = nf_def_var(wr_ncid, 'h_s', NF_DOUBLE, 1, dimlist, wrVarIDh_s)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'u', NF_DOUBLE, 3, dimlist, wrVarIDu)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ nferr = nf_def_var(wr_ncid, '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, 'u_src', NF_DOUBLE, 2, dimlist, wrVarIDu_src)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'v', NF_DOUBLE, 3, dimlist, wrVarIDv)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnCells
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'h', NF_DOUBLE, 3, dimlist, wrVarIDh)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnEdges
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'vh', NF_DOUBLE, 3, dimlist, wrVarIDvh)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnVertices
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'circulation', NF_DOUBLE, 3, dimlist, wrVarIDcirculation)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnVertices
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'vorticity', NF_DOUBLE, 3, dimlist, wrVarIDvorticity)
+ dimlist( 1) = wrDimIDnVertLevels
+ dimlist( 2) = wrDimIDnCells
+ dimlist( 3) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'ke', NF_DOUBLE, 3, dimlist, wrVarIDke)
+ dimlist( 1) = wrDimIDnTracers
+ dimlist( 2) = wrDimIDnVertLevels
+ dimlist( 3) = wrDimIDnCells
+ dimlist( 4) = wrDimIDTime
+ nferr = nf_def_var(wr_ncid, 'tracers', NF_DOUBLE, 4, dimlist, wrVarIDtracers)
+
+ nferr = nf_enddef(wr_ncid)
+
+ end subroutine write_netcdf_init
+
+
+ subroutine write_netcdf_fields( &
+ time, &
+ latCell, &
+ lonCell, &
+ xCell, &
+ yCell, &
+ zCell, &
+ indexToCellID, &
+ latEdge, &
+ lonEdge, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ indexToEdgeID, &
+ latVertex, &
+ lonVertex, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex, &
+ h_s, &
+ boundaryEdge, &
+ boundaryVertex, &
+ u_src, &
+ u, &
+ v, &
+ h, &
+ vh, &
+ circulation, &
+ vorticity, &
+ ke, &
+ tracers &
+ )
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer, intent(in) :: time
+ real (kind=8), dimension(:), intent(in) :: latCell
+ real (kind=8), dimension(:), intent(in) :: lonCell
+ real (kind=8), dimension(:), intent(in) :: xCell
+ real (kind=8), dimension(:), intent(in) :: yCell
+ real (kind=8), dimension(:), intent(in) :: zCell
+ integer, dimension(:), intent(in) :: indexToCellID
+ real (kind=8), dimension(:), intent(in) :: latEdge
+ real (kind=8), dimension(:), intent(in) :: lonEdge
+ real (kind=8), dimension(:), intent(in) :: xEdge
+ real (kind=8), dimension(:), intent(in) :: yEdge
+ real (kind=8), dimension(:), intent(in) :: zEdge
+ integer, dimension(:), intent(in) :: indexToEdgeID
+ real (kind=8), dimension(:), intent(in) :: latVertex
+ real (kind=8), dimension(:), intent(in) :: lonVertex
+ real (kind=8), dimension(:), intent(in) :: xVertex
+ real (kind=8), dimension(:), intent(in) :: yVertex
+ real (kind=8), dimension(:), intent(in) :: zVertex
+ integer, dimension(:), intent(in) :: indexToVertexID
+ integer, dimension(:,:), intent(in) :: cellsOnEdge
+ integer, dimension(:), intent(in) :: nEdgesOnCell
+ integer, dimension(:), intent(in) :: nEdgesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnCell
+ integer, dimension(:,:), intent(in) :: edgesOnEdge
+ real (kind=8), dimension(:,:), intent(in) :: weightsOnEdge
+ real (kind=8), dimension(:), intent(in) :: dvEdge
+ real (kind=8), dimension(:), intent(in) :: dcEdge
+ real (kind=8), dimension(:), intent(in) :: angleEdge
+ real (kind=8), dimension(:), intent(in) :: areaCell
+ real (kind=8), dimension(:), intent(in) :: areaTriangle
+ integer, dimension(:,:), intent(in) :: cellsOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnCell
+ integer, dimension(:,:), intent(in) :: verticesOnEdge
+ integer, dimension(:,:), intent(in) :: edgesOnVertex
+ integer, dimension(:,:), intent(in) :: cellsOnVertex
+ real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
+ real (kind=8), dimension(:), intent(in) :: fEdge
+ real (kind=8), dimension(:), intent(in) :: fVertex
+ real (kind=8), dimension(:), intent(in) :: h_s
+ integer, dimension(:,:), intent(in) :: boundaryEdge
+ integer, dimension(:,:), intent(in) :: boundaryVertex
+ real (kind=8), dimension(:,:), intent(in) :: u_src
+ real (kind=8), dimension(:,:,:), intent(in) :: u
+ real (kind=8), dimension(:,:,:), intent(in) :: v
+ real (kind=8), dimension(:,:,:), intent(in) :: h
+ real (kind=8), dimension(:,:,:), intent(in) :: vh
+ real (kind=8), dimension(:,:,:), intent(in) :: circulation
+ real (kind=8), dimension(:,:,:), intent(in) :: vorticity
+ real (kind=8), dimension(:,:,:), intent(in) :: ke
+ real (kind=8), dimension(:,:,:,:), intent(in) :: tracers
+
+ integer :: nferr
+ integer, dimension(1) :: start1, count1
+ integer, dimension(2) :: start2, count2
+ integer, dimension(3) :: start3, count3
+ integer, dimension(4) :: start4, count4
+
+ start1(1) = 1
+
+ start2(1) = 1
+ start2(2) = 1
+
+ start3(1) = 1
+ start3(2) = 1
+ start3(3) = 1
+
+ start4(1) = 1
+ start4(2) = 1
+ start4(3) = 1
+ start4(4) = 1
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatCell, start1, count1, latCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonCell, start1, count1, lonCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxCell, start1, count1, xCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyCell, start1, count1, yCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzCell, start1, count1, zCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToCellID, start1, count1, indexToCellID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatEdge, start1, count1, latEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonEdge, start1, count1, lonEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxEdge, start1, count1, xEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyEdge, start1, count1, yEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzEdge, start1, count1, zEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToEdgeID, start1, count1, indexToEdgeID)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlatVertex, start1, count1, latVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlonVertex, start1, count1, lonVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDxVertex, start1, count1, xVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDyVertex, start1, count1, yVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDzVertex, start1, count1, zVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDindexToVertexID, start1, count1, indexToVertexID)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnEdge, start2, count2, cellsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnCell, start1, count1, nEdgesOnCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDnEdgesOnEdge, start1, count1, nEdgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnCell, start2, count2, edgesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnEdge, start2, count2, edgesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = 2*wrLocalmaxEdges
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDweightsOnEdge, start2, count2, weightsOnEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdvEdge, start1, count1, dvEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDdcEdge, start1, count1, dcEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDangleEdge, start1, count1, angleEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaCell, start1, count1, areaCell)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDareaTriangle, start1, count1, areaTriangle)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnCell, start2, count2, cellsOnCell)
+
+ start2(2) = 1
+ count2( 1) = wrLocalmaxEdges
+ count2( 2) = wrLocalnCells
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnCell, start2, count2, verticesOnCell)
+
+ start2(2) = 1
+ count2( 1) = 2
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDverticesOnEdge, start2, count2, verticesOnEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDedgesOnVertex, start2, count2, edgesOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDcellsOnVertex, start2, count2, cellsOnVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertexDegree
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfEdge, start1, count1, fEdge)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnVertices
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDfVertex, start1, count1, fVertex)
+
+ start1(1) = 1
+ count1( 1) = wrLocalnCells
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDh_s, start1, count1, h_s)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryEdge, start2, count2, boundaryEdge)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnVertices
+ nferr = nf_put_vara_int(wr_ncid, wrVarIDboundaryVertex, start2, count2, boundaryVertex)
+
+ start2(2) = 1
+ count2( 1) = wrLocalnVertLevels
+ count2( 2) = wrLocalnEdges
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDu_src, start2, count2, u_src)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDu, start3, count3, u)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDv, start3, count3, v)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnCells
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDh, start3, count3, h)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnEdges
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDvh, start3, count3, vh)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnVertices
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDcirculation, start3, count3, circulation)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnVertices
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDvorticity, start3, count3, vorticity)
+
+ start3(3) = time
+ count3( 1) = wrLocalnVertLevels
+ count3( 2) = wrLocalnCells
+ count3( 3) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDke, start3, count3, ke)
+
+ start4(4) = time
+ count4( 1) = wrLocalnTracers
+ count4( 2) = wrLocalnVertLevels
+ count4( 3) = wrLocalnCells
+ count4( 4) = 1
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDtracers, start4, count4, tracers)
+
+ end subroutine write_netcdf_fields
+
+
+ subroutine write_netcdf_finalize()
+
+ implicit none
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_close(wr_ncid)
+
+ end subroutine write_netcdf_finalize
+
+end module write_netcdf
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/namelist.input
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/namelist.input (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/namelist.input 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,11 @@
+&periodic_grid
+ nx = 100,
+ ny = 116,
+ nVertLevels = 1,
+ nTracers = 2,
+ nVertexDegree = 4,
+ nproc = 2, 4, 8,
+ POP_GridFileName = 'grid.final.da'
+ POP_KmtFileName = 'topo_file_no_marginal_seas'
+ POP_DepthFileName = 'deptht.r8'
+/
Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas.F 2010-04-06 16:25:23 UTC (rev 180)
@@ -0,0 +1,1295 @@
+program pop_to_mpas
+
+ use cell_indexing
+ use write_netcdf
+ use sphere_utilities
+ use data_types
+
+ implicit none
+
+! specify grid resolution with variable dc (meters)
+ real (kind=8), parameter :: dc = 10000.0 ! Distance between cell centers
+
+ real (kind=8), parameter :: pi = 3.141592653589793
+ real (kind=8), parameter :: ONE = 1.0_8
+ real (kind=8), parameter :: TWO = 2.0_8
+ real (kind=8), parameter :: FOUR = 4.0_8
+ real (kind=8), parameter :: omega = 7.292123625e-5_8 ! angular vel. of Earth 1/s
+!maltrud set radius
+! NOTE-- did not work with real value of a--kite areas were bad
+! so apply at the end
+ real (kind=8), parameter :: a = 6371220.0 ! Radius of Earth (m)
+ real, parameter :: a_unit = 1.0 ! Radius of unit sphere
+ real, parameter :: xscale = 1.01
+
+ integer, allocatable, dimension(:) :: indexToCellID, indexToEdgeID, indexToVertexID
+ integer, allocatable, dimension(:) :: nEdgesOnCell, nEdgesOnEdge
+ integer, allocatable, dimension(:,:) :: cellsOnCell, edgesOnCell, verticesOnCell
+ integer, allocatable, dimension(:,:) :: cellsOnEdge, edgesOnEdge, verticesOnEdge
+ integer, allocatable, dimension(:,:) :: edgesOnVertex, cellsOnVertex, boundaryEdge, boundaryVertex
+ real (kind=8), allocatable, dimension(:) :: areaTriangle, areaCell, angleEdge
+ real (kind=8), allocatable, dimension(:) :: dcEdge, dvEdge
+ real (kind=8), allocatable, dimension(:) :: latCell, lonCell, xCell, yCell, zCell
+ real (kind=8), allocatable, dimension(:) :: latEdge, lonEdge, xEdge, yEdge, zEdge
+ real (kind=8), allocatable, dimension(:) :: latVertex, lonVertex, xVertex, yVertex, zVertex
+ real (kind=8), allocatable, dimension(:,:) :: weightsOnEdge, kiteAreasOnVertex, u_src
+ real (kind=8), allocatable, dimension(:) :: fEdge, fVertex, h_s
+ real (kind=8), allocatable, dimension(:,:,:) :: u, v, h, vh, circulation, vorticity, ke
+ real (kind=8), allocatable, dimension(:,:,:,:) :: tracers
+
+ integer :: i, j, np, iCell, iVertex, iEdge, iEdge_s, iVertex_s
+ integer :: nCells, nEdges, nVertices
+ integer :: iRow, iCol, ii, jj, jm1, jp1, k
+ integer :: nprocx, nprocy, itmp(20), m
+ real (kind=8) :: r
+ character (len=32) :: decomp_fname
+
+ integer :: iVert1, iVert2, iEdge1, iEdge2, iCell1, iCell2
+ integer :: iVert3, iVert4, iEdge3, iEdge4, iCell3, iCell4
+ real (kind=8) :: ymid, f0, beta, ytmp, ymax, stress, distance &
+ ,AreaLeft, AreaRight, lonPerturb0, latPerturb0 &
+ ,scalePerturb, heightPerturb, MaxDepth
+ real, dimension(5,3) :: AvgXYZ
+ real :: RadiusTest
+
+ integer, allocatable, dimension(:,:) :: POP_KMT, cellID_ij
+ integer, allocatable, dimension(:,:,:) :: edgeID_ij, vertexID_ij
+ real (kind=8), allocatable, dimension(:,:) :: &
+ POP_uLat, POP_uLon, POP_Angle, POP_Depth
+ type (geo_point) :: latlon1, latlon2, latlon3, latlon4, latlon5
+ integer :: recLength, ip1, im1
+
+!
+! begin
+!
+ call cell_indexing_read_nl()
+
+ allocate(POP_uLat(nx,ny), POP_uLon(nx,ny), POP_Angle(nx,ny))
+ allocate(POP_Depth(nx,ny))
+ allocate(POP_KMT(nx,ny))
+ allocate(cellID_ij(nx,ny))
+ allocate(vertexID_ij(nx,ny,nVertexDegree))
+ allocate(edgeID_ij(nx,ny,nVertexDegree))
+!
+! open and read ULAT, ULON, ANGLE from POP grid file
+!
+ inquire (iolength=recLength) POP_uLat
+ open(57, file = trim(POP_GridFileName), access = 'direct', &
+ recl = recLength, form = 'unformatted', status = 'old')
+ read(57,rec=1) POP_uLat
+ read(57,rec=2) POP_uLon
+ read(57,rec=7) POP_Angle
+ close(57)
+ write(6,*) 'lat/lon',minval(POP_uLat),maxval(POP_uLat),minval(POP_uLon),maxval(POP_uLon)
+!
+! open and read POP depth file
+!
+ if (trim(POP_DepthFileName) /= 'unknown_POP_DepthFileName') then
+ inquire (iolength=recLength) POP_Depth
+ open(57, file = trim(POP_DepthFileName), access = 'direct', &
+ recl = recLength, form = 'unformatted', status = 'old')
+ read(57,rec=1) POP_Depth
+ close(57)
+ endif
+ write(6,*) ' POP_Depth', minval(POP_Depth), maxval(POP_Depth)
+
+!
+! open and read POP KMT file
+!
+ inquire (iolength=recLength) POP_KMT
+ open(57, file = trim(POP_KmtFileName), access = 'direct', &
+ recl = recLength, form = 'unformatted', status = 'old')
+ read(57,rec=1) POP_KMT
+ close(57)
+ write(6,*) ' POP_KMT',minval(POP_KMT), maxval(POP_KMT)
+
+ iCell = 1
+ iEdge = 1
+ iVertex = 1
+ cellID_ij = 0
+ edgeID_ij = 0
+ vertexID_ij = 0
+
+ iEdge_s = 0
+ iVertex_s = 0
+!
+! skip bottom and top rows since they are land
+! NOTE: tripole will be different for top row
+!
+ do iRow = 2, ny - 1
+ do iCol = 1, nx
+ if (POP_KMT(iCol,iRow) /= 0) then
+ cellID_ij(iCol,iRow) = iCell
+ iCell = iCell + 1
+ endif
+ enddo
+ enddo
+
+ nCells = iCell - 1
+ write(*,*)' done labelling cells; nCells = ', nCells
+ write(*,*)' max value of cellID_ij ', maxval(cellID_ij)
+
+ do iRow = 2, ny - 1
+ do iCol = 1, nx
+ iCell = cellID_ij(iCol,iRow)
+ im1 = merge(nx,iCol-1,iCol==1)
+ ip1 = merge(1 ,iCol+1,iCol==nx)
+ jm1 = iRow - 1
+ jp1 = iRow + 1
+
+ if(iCell.eq.0) cycle
+
+ ! non-zero iCell always counts edge to the right
+ j=3
+ edgeID_ij(iCol,iRow,j) = iEdge; iEdge = iEdge + 1
+
+ ! non-zero iCell always counts edge to the top
+ j=4
+ edgeID_ij(iCol,iRow,j) = iEdge; iEdge = iEdge + 1
+
+ ! check cell to the left
+ j=1
+ if(cellID_ij(im1, iRow).eq.0) then
+ edgeID_ij(iCol, iRow, j) = iEdge; iEdge = iEdge + 1
+ else
+ edgeID_ij(iCol, iRow, j) = edgeID_ij(im1, iRow, 3)
+ endif
+
+ ! check cell to the bottom
+ j=2
+ if(cellID_ij(iCol, jm1).eq.0) then
+ edgeID_ij(iCol, iRow, j) = iEdge; iEdge = iEdge + 1
+ else
+ edgeID_ij(iCol, iRow, j) = edgeID_ij(iCol, jm1, 4)
+ endif
+
+ ! non-zero iCells always owns vertex up and to the right
+ k=3
+ vertexID_ij(iCol,iRow,k) = iVertex; iVertex = iVertex + 1
+
+ ! check cell to the left
+ k=4
+ if(cellID_ij(im1,iRow).eq.0) then
+ vertexID_ij(iCol,iRow,k) = iVertex; iVertex = iVertex + 1
+ else
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,iRow,3)
+ endif
+
+ ! check cell to the bottom and left
+ k=1
+ if(cellID_ij(im1,iRow).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,iRow,2)
+ elseif(cellID_ij(im1,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,jm1,3)
+ elseif(cellID_ij(iCol,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(iCol,jm1,4)
+ else
+ vertexID_ij(iCol,iRow,k) = iVertex; iVertex = iVertex + 1
+ endif
+
+ ! check cell to the bottom
+ k=2
+ if(cellID_ij(iCol,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(iCol,jm1,3)
+ elseif(cellID_ij(ip1,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(ip1,jm1,4)
+ else
+ vertexID_ij(iCol,iRow,k) = iVertex; iVertex = iVertex + 1
+ endif
+
+ enddo
+ enddo
+
+ do iRow = 2, ny - 1
+ do iCol = 1, nx
+ iCell = cellID_ij(iCol,iRow)
+ im1 = merge(nx,iCol-1,iCol==1)
+ ip1 = merge(1 ,iCol+1,iCol==nx)
+ jm1 = iRow - 1
+ jp1 = iRow + 1
+
+ if(iCell.eq.0) cycle
+
+ ! check cell to the left
+ j=1
+ if(cellID_ij(im1, iRow).ne.0) then
+ if(edgeID_ij(iCol, iRow, j).eq.0) then
+ edgeID_ij(iCol, iRow, j) = edgeID_ij(im1, iRow, 3)
+ endif
+ endif
+
+ ! check cell to the bottom
+ j=2
+ if(cellID_ij(iCol, jm1).ne.0) then
+ if(edgeID_ij(iCol, iRow, j).eq.0) then
+ edgeID_ij(iCol, iRow, j) = edgeID_ij(iCol, jm1, 4)
+ endif
+ endif
+
+ ! check cell to the left
+ k=4
+ if(vertexID_ij(iCol,iRow,k).eq.0) then
+ if(cellID_ij(im1,iRow).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,iRow,3)
+ endif
+ endif
+
+ ! check cell to the bottom and left
+ k=1
+ if(vertexID_ij(iCol,iRow,k).eq.0) then
+ if(cellID_ij(im1,iRow).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,iRow,2)
+ elseif(cellID_ij(im1,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(im1,jm1,3)
+ elseif(cellID_ij(iCol,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(iCol,jm1,4)
+ endif
+ endif
+
+ ! check cell to the bottom
+ k=2
+ if(vertexID_ij(iCol,iRow,k).eq.0) then
+ if(cellID_ij(iCol,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(iCol,jm1,3)
+ elseif(cellID_ij(ip1,jm1).ne.0) then
+ vertexID_ij(iCol,iRow,k) = vertexID_ij(ip1,jm1,4)
+ endif
+ endif
+
+ enddo
+ enddo
+
+
+ do iRow=2,ny-1
+ do iCol=1,nx
+ iCell = cellID_ij(iCol,iRow)
+ if(iCell.eq.0) cycle
+ if(minval(edgeID_ij(iCol,iRow,:)).le.0) then
+ write(6,*) edgeID_ij(iCol,iRow,:)
+ endif
+ enddo
+ enddo
+
+ nEdges = iEdge - 1
+ nVertices = iVertex - 1
+ write(*,*) ' number of cells ', nCells
+ write(*,*)' done labelling edges, vertices; nEdges, nVertices = ', &
+ nEdges, nVertices
+
+ ! make sure each ID is only counted once
+ do k=1,nCells
+ m = 0
+ do iRow=2,ny-1
+ do iCol=1,nx
+ if(cellID_ij(iCol,iRow).eq.k) m=m+1
+ enddo
+ enddo
+ if(m.ne.1) then
+ write(6,*) ' cellID_ij invalid', iCell, m
+ stop
+ endif
+ enddo
+
+ ! make sure each ID is only counted once
+ do k=1,nEdges
+ m = 0
+ do iRow=2,ny-1
+ do iCol=1,nx
+ do j=1,nVertexDegree
+ if(edgeID_ij(iCol,iRow,j).eq.k) m=m+1
+ enddo
+ enddo
+ enddo
+ if(m.gt.2.or.m.lt.1) then
+ write(6,*) ' edgeID_ij invalid', iEdge, m
+ stop
+ endif
+ enddo
+
+ ! make sure each ID is only counted once
+ do k=1,nVertices
+ m = 0
+ do iRow=2,ny-1
+ do iCol=1,nx
+ do j=1,nVertexDegree
+ if(vertexID_ij(iCol,iRow,j).eq.k) m=m+1
+ enddo
+ enddo
+ enddo
+ if(m.gt.4.or.m.lt.1) then
+ write(6,*) ' vertexID_ij invalid', iVertex, m
+ stop
+ endif
+ enddo
+
+
+ allocate(indexToCellID(nCells))
+ allocate(indexToEdgeID(nEdges))
+ allocate(indexToVertexID(nVertices))
+
+ allocate(nEdgesOnCell(nCells))
+ allocate(cellsOnCell(maxEdges, nCells))
+ allocate(edgesOnCell(maxEdges, nCells))
+ allocate(verticesOnCell(maxEdges, nCells))
+
+ allocate(nEdgesOnEdge(nEdges))
+ allocate(cellsOnEdge(2,nEdges))
+ allocate(verticesOnEdge(2,nEdges))
+ allocate(edgesOnEdge(2*maxEdges,nEdges))
+ allocate(weightsOnEdge(2*maxEdges,nEdges))
+
+ allocate(edgesOnVertex(nVertexDegree,nVertices))
+ allocate(cellsOnVertex(nVertexDegree,nVertices))
+ allocate(kiteAreasOnVertex(nVertexDegree,nVertices))
+
+ allocate(areaTriangle(nVertices))
+ allocate(areaCell(nCells))
+
+ allocate(dcEdge(nEdges))
+ allocate(dvEdge(nEdges))
+ allocate(angleEdge(nEdges))
+
+ allocate(latCell(nCells))
+ allocate(lonCell(nCells))
+ allocate(xCell(nCells))
+ allocate(yCell(nCells))
+ allocate(zCell(nCells))
+ allocate(latEdge(nEdges))
+ allocate(lonEdge(nEdges))
+ allocate(xEdge(nEdges))
+ allocate(yEdge(nEdges))
+ allocate(zEdge(nEdges))
+ allocate(latVertex(nVertices))
+ allocate(lonVertex(nVertices))
+ allocate(xVertex(nVertices))
+ allocate(yVertex(nVertices))
+ allocate(zVertex(nVertices))
+
+ allocate(fEdge(nEdges))
+ allocate(fVertex(nVertices))
+ allocate(h_s(nCells))
+ allocate(boundaryEdge(nVertLevels, nEdges))
+ allocate(boundaryVertex(nVertLevels, nEdges))
+ allocate(u_src(nVertLevels, nEdges))
+
+ allocate(u(nVertLevels,nEdges,1))
+ allocate(v(nVertLevels,nEdges,1))
+ allocate(vh(nVertLevels,nEdges,1))
+ allocate(h(nVertLevels,nCells,1))
+ allocate(circulation(nVertLevels,nVertices,1))
+ allocate(vorticity(nVertLevels,nVertices,1))
+ allocate(ke(nVertLevels,nCells,1))
+ allocate(tracers(nTracers,nVertLevels,nCells,1))
+
+ latCell = 0.0
+ lonCell = 0.0
+ latEdge = 0.0
+ lonEdge = 0.0
+ latVertex = 0.0
+ lonVertex = 0.0
+!
+! now that all cells, edges and vertices have been labelled, we can set
+! up connectivity, etc.
+!
+ nEdgesOnEdge = 0
+ nEdgesOnCell = 0
+
+ do iRow = 2, ny - 1
+ do iCol = 1, nx
+ iCell = cellID_ij(iCol,iRow)
+
+ if(iCell.eq.0) cycle
+
+ iEdge1 = edgeID_ij(iCol,iRow,1)
+ iEdge2 = edgeID_ij(iCol,iRow,2)
+ iEdge3 = edgeID_ij(iCol,iRow,3)
+ iEdge4 = edgeID_ij(iCol,iRow,4)
+
+ iVert1 = vertexID_ij(iCol,iRow,1)
+ iVert2 = vertexID_ij(iCol,iRow,2)
+ iVert3 = vertexID_ij(iCol,iRow,3)
+ iVert4 = vertexID_ij(iCol,iRow,4)
+
+ do j=1,nVertexDegree
+ do k=1,nVertexDegree
+ if(j.ne.k) then
+ if(edgeID_ij(iCol,iRow,j).eq.edgeID_ij(iCol,iRow,k)) then
+ write(6,*) 'edgeID invalid'
+ write(6,*) edgeID_ij(iCol,iRow,:)
+ stop
+ endif
+ if(vertexID_ij(iCol,iRow,j).eq.vertexID_ij(iCol,iRow,k)) then
+ write(6,*) 'vertexID invalid'
+ write(6,*) vertexID_ij(iCol,iRow,:)
+ stop
+ endif
+ endif
+ enddo
+ enddo
+
+ nEdgesOnCell(iCell) = nVertexDegree
+
+ do j = 1, nVertexDegree
+ edgesOnCell(j,iCell) = edgeID_ij(iCol,iRow,j)
+ verticesOnCell(j,iCell) = vertexID_ij(iCol,iRow,j)
+ enddo
+
+ im1 = merge(nx,iCol-1,iCol==1)
+ ip1 = merge(1 ,iCol+1,iCol==nx)
+ jm1 = iRow - 1
+ jp1 = iRow + 1
+
+ ! iCell is to the right of iEdge1 and to the top of iEdge2
+ cellsOnEdge(2,iEdge1) = iCell
+ cellsOnEdge(2,iEdge2) = iCell
+
+ ! iCell is to the left of iEdge3 and to the bottom of iEdge 4
+ cellsOnEdge(1,iEdge3) = iCell
+ cellsOnEdge(1,iEdge4) = iCell
+
+ ! iCell is down and to the left of iVert3
+ cellsOnVertex(1,iVert3) = cellID_ij(iCol,iRow)
+
+ ! iCell is down and to the right of iVert4
+ cellsOnVertex(2,iVert4) = cellID_ij(iCol,iRow)
+
+ ! iCell is up and the to the right of iVert1
+ cellsOnVertex(3,iVert1) = cellID_ij(iCol,iRow)
+
+ ! iCell is up and to the left of iVert2
+ cellsOnVertex(4,iVert2) = cellID_ij(iCol,iRow)
+
+ ! fill in cellsOnCell
+ cellsOnCell(1,iCell) = cellID_ij(im1,iRow)
+ cellsOnCell(2,iCell) = cellID_ij(iCol,jm1)
+ cellsOnCell(3,iCell) = cellID_ij(ip1,iRow)
+ cellsOnCell(4,iCell) = cellID_ij(iCol,jp1)
+
+ latVertex(iVert1) = POP_uLat(im1, iRow-1)
+ lonVertex(iVert1) = POP_uLon(im1, iRow-1)
+ latVertex(iVert2) = POP_uLat(iCol,iRow-1)
+ lonVertex(iVert2) = POP_uLon(iCol,iRow-1)
+ latVertex(iVert3) = POP_uLat(iCol,iRow)
+ lonVertex(iVert3) = POP_uLon(iCol,iRow)
+ latVertex(iVert4) = POP_uLat(im1, iRow)
+ lonVertex(iVert4) = POP_uLon(im1, iRow)
+
+ edgesOnVertex(1,iVert1) = iEdge1
+ edgesOnVertex(4,iVert1) = iEdge2
+ edgesOnVertex(2,iVert2) = iEdge2
+ edgesOnVertex(1,iVert2) = iEdge3
+ edgesOnVertex(2,iVert3) = iEdge4
+ edgesOnVertex(3,iVert3) = iEdge3
+ edgesOnVertex(3,iVert4) = iEdge1
+ edgesOnVertex(4,iVert4) = iEdge4
+
+ verticesOnEdge(1,edgesOnCell(1,iCell)) = verticesOnCell(1,iCell)
+ verticesOnEdge(2,edgesOnCell(1,iCell)) = verticesOnCell(4,iCell)
+ verticesOnEdge(1,edgesOnCell(2,iCell)) = verticesOnCell(2,iCell)
+ verticesOnEdge(2,edgesOnCell(2,iCell)) = verticesOnCell(1,iCell)
+
+ verticesOnEdge(1,edgesOnCell(3,iCell)) = verticesOnCell(2,iCell)
+ verticesOnEdge(2,edgesOnCell(3,iCell)) = verticesOnCell(3,iCell)
+ verticesOnEdge(1,edgesOnCell(4,iCell)) = verticesOnCell(3,iCell)
+ verticesOnEdge(2,edgesOnCell(4,iCell)) = verticesOnCell(4,iCell)
+
+ edgesOnEdge(1,edgesOnCell(3,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(2,edgesOnCell(3,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(3,edgesOnCell(3,iCell)) = edgesOnCell(2,iCell)
+
+ edgesOnEdge(1,edgesOnCell(4,iCell)) = edgesOnCell(1,iCell)
+ edgesOnEdge(2,edgesOnCell(4,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(3,edgesOnCell(4,iCell)) = edgesOnCell(3,iCell)
+
+ edgesOnEdge(4,edgesOnCell(1,iCell)) = edgesOnCell(2,iCell)
+ edgesOnEdge(5,edgesOnCell(1,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(6,edgesOnCell(1,iCell)) = edgesOnCell(4,iCell)
+
+ edgesOnEdge(4,edgesOnCell(2,iCell)) = edgesOnCell(3,iCell)
+ edgesOnEdge(5,edgesOnCell(2,iCell)) = edgesOnCell(4,iCell)
+ edgesOnEdge(6,edgesOnCell(2,iCell)) = edgesOnCell(1,iCell)
+
+ weightsOnEdge(1,edgesOnCell(3,iCell)) = ONE / FOUR
+ weightsOnEdge(2,edgesOnCell(3,iCell)) = 0.0
+ weightsOnEdge(3,edgesOnCell(3,iCell)) = ONE / FOUR
+
+ weightsOnEdge(1,edgesOnCell(4,iCell)) = -ONE / FOUR
+ weightsOnEdge(2,edgesOnCell(4,iCell)) = 0.0
+ weightsOnEdge(3,edgesOnCell(4,iCell)) = -ONE / FOUR
+
+ weightsOnEdge(4,edgesOnCell(1,iCell)) = ONE / FOUR
+ weightsOnEdge(5,edgesOnCell(1,iCell)) = 0.0
+ weightsOnEdge(6,edgesOnCell(1,iCell)) = ONE / FOUR
+
+ weightsOnEdge(4,edgesOnCell(2,iCell)) = -ONE / FOUR
+ weightsOnEdge(5,edgesOnCell(2,iCell)) = 0.0
+ weightsOnEdge(6,edgesOnCell(2,iCell)) = -ONE / FOUR
+
+ end do
+ end do
+
+ do iEdge=1,nEdges
+ m = 0
+ do iVertex=1,nVertices
+ do j=1,nVertexDegree
+ if(edgesOnVertex(j,iVertex).eq.iEdge) m=m+1
+ enddo
+ enddo
+
+ if(m.ne.2) then
+ write(6,*) m,iEdge
+ do iVertex=1,nVertices
+ do j=1,nVertexDegree
+ if(edgesOnVertex(j,iVertex).eq.iEdge) then
+ write(6,*) iVertex
+ write(6,*) ' ', edgesOnVertex(:,iVertex)
+ endif
+ enddo
+ enddo
+ write(6,*)
+ stop
+ endif
+ enddo
+ write(6,*) ' edgesOnVertex passed sanity check'
+ write(6,*)
+
+
+! convert lat/lon to XYZ at each vertex
+ do iVertex = 1,nVertices
+ latlon1%lat = latVertex(iVertex)
+ latlon1%lon = lonVertex(iVertex)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+ xVertex(iVertex) = AvgXYZ(1,1)
+ yVertex(iVertex) = AvgXYZ(1,2)
+ zVertex(iVertex) = AvgXYZ(1,3)
+ enddo
+
+!
+! average vertex coordinates to create edge coordinates
+!
+ do iCell = 1, nCells
+!
+! edge 1
+!
+ iVert1 = verticesOnEdge(1,edgesOnCell(1,iCell))
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+
+ iVert2 = verticesOnEdge(2,edgesOnCell(1,iCell))
+ latlon2%lat = latVertex(iVert2)
+ latlon2%lon = lonVertex(iVert2)
+ call convert_lx(AvgXYZ(2,1), AvgXYZ(2,2), AvgXYZ(2,3), a_unit, latlon2)
+
+ AvgXYZ(3,:) = (AvgXYZ(1,:) + AvgXYZ(2,:))/TWO
+ RadiusTest = sqrt(AvgXYZ(3,1)**2 + AvgXYZ(3,2)**2 + AvgXYZ(3,3)**2)
+ AvgXYZ(3,:) = AvgXYZ(3,:)*a_unit/RadiusTest
+
+ xEdge(edgesOnCell(1,iCell)) = AvgXYZ(3,1)
+ yEdge(edgesOnCell(1,iCell)) = AvgXYZ(3,2)
+ zEdge(edgesOnCell(1,iCell)) = AvgXYZ(3,3)
+
+! calculate distance between vertices
+ distance = sphere_distance(latlon1, latlon2, a_unit)
+ dvEdge(edgesOnCell(1,iCell)) = distance
+
+! convert back to lat,lon
+ call convert_xl(AvgXYZ(3,1), AvgXYZ(3,2), AvgXYZ(3,3), latlon2)
+ latEdge(edgesOnCell(1,iCell)) = latlon2%lat
+ lonEdge(edgesOnCell(1,iCell)) = latlon2%lon
+
+!
+! edge 2
+!
+ iVert1 = verticesOnEdge(1,edgesOnCell(2,iCell))
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+
+ iVert2 = verticesOnEdge(2,edgesOnCell(2,iCell))
+ latlon2%lat = latVertex(iVert2)
+ latlon2%lon = lonVertex(iVert2)
+ call convert_lx(AvgXYZ(2,1), AvgXYZ(2,2), AvgXYZ(2,3), a_unit, latlon2)
+
+ AvgXYZ(3,:) = (AvgXYZ(1,:) + AvgXYZ(2,:))/TWO
+ RadiusTest = sqrt(AvgXYZ(3,1)**2 + AvgXYZ(3,2)**2 + AvgXYZ(3,3)**2)
+ AvgXYZ(3,:) = AvgXYZ(3,:)*a_unit/RadiusTest
+
+ xEdge(edgesOnCell(2,iCell)) = AvgXYZ(3,1)
+ yEdge(edgesOnCell(2,iCell)) = AvgXYZ(3,2)
+ zEdge(edgesOnCell(2,iCell)) = AvgXYZ(3,3)
+
+! calculate distance between vertices
+ distance = sphere_distance(latlon1, latlon2, a_unit)
+ dvEdge(edgesOnCell(2,iCell)) = distance
+
+! convert back to lat,lon
+ call convert_xl(AvgXYZ(3,1), AvgXYZ(3,2), AvgXYZ(3,3), latlon2)
+ latEdge(edgesOnCell(2,iCell)) = latlon2%lat
+ lonEdge(edgesOnCell(2,iCell)) = latlon2%lon
+
+!
+! edge 3
+!
+ iVert1 = verticesOnEdge(1,edgesOnCell(3,iCell))
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+
+ iVert2 = verticesOnEdge(2,edgesOnCell(3,iCell))
+ latlon2%lat = latVertex(iVert2)
+ latlon2%lon = lonVertex(iVert2)
+ call convert_lx(AvgXYZ(2,1), AvgXYZ(2,2), AvgXYZ(2,3), a_unit, latlon2)
+
+ AvgXYZ(3,:) = (AvgXYZ(1,:) + AvgXYZ(2,:))/TWO
+ RadiusTest = sqrt(AvgXYZ(3,1)**2 + AvgXYZ(3,2)**2 + AvgXYZ(3,3)**2)
+ AvgXYZ(3,:) = AvgXYZ(3,:)*a_unit/RadiusTest
+
+ xEdge(edgesOnCell(3,iCell)) = AvgXYZ(3,1)
+ yEdge(edgesOnCell(3,iCell)) = AvgXYZ(3,2)
+ zEdge(edgesOnCell(3,iCell)) = AvgXYZ(3,3)
+
+! calculate distance between vertices
+ distance = sphere_distance(latlon1, latlon2, a_unit)
+ dvEdge(edgesOnCell(3,iCell)) = distance
+
+! convert back to lat,lon
+ call convert_xl(AvgXYZ(3,1), AvgXYZ(3,2), AvgXYZ(3,3), latlon2)
+ latEdge(edgesOnCell(3,iCell)) = latlon2%lat
+ lonEdge(edgesOnCell(3,iCell)) = latlon2%lon
+
+!
+! edge 4
+!
+ iVert1 = verticesOnEdge(1,edgesOnCell(4,iCell))
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+
+ iVert2 = verticesOnEdge(2,edgesOnCell(4,iCell))
+ latlon2%lat = latVertex(iVert2)
+ latlon2%lon = lonVertex(iVert2)
+ call convert_lx(AvgXYZ(2,1), AvgXYZ(2,2), AvgXYZ(2,3), a_unit, latlon2)
+
+ AvgXYZ(3,:) = (AvgXYZ(1,:) + AvgXYZ(2,:))/TWO
+ RadiusTest = sqrt(AvgXYZ(3,1)**2 + AvgXYZ(3,2)**2 + AvgXYZ(3,3)**2)
+ AvgXYZ(3,:) = AvgXYZ(3,:)*a_unit/RadiusTest
+
+ xEdge(edgesOnCell(4,iCell)) = AvgXYZ(3,1)
+ yEdge(edgesOnCell(4,iCell)) = AvgXYZ(3,2)
+ zEdge(edgesOnCell(4,iCell)) = AvgXYZ(3,3)
+
+! calculate distance between vertices
+ distance = sphere_distance(latlon1, latlon2, a_unit)
+ dvEdge(edgesOnCell(4,iCell)) = distance
+
+! convert back to lat,lon
+ call convert_xl(AvgXYZ(3,1), AvgXYZ(3,2), AvgXYZ(3,3), latlon2)
+ latEdge(edgesOnCell(4,iCell)) = latlon2%lat
+ lonEdge(edgesOnCell(4,iCell)) = latlon2%lon
+
+ enddo
+
+!
+! average vertex coordinates to create centered coordinates
+!
+ do iCell = 1, nCells
+
+ iVert1 = verticesOnCell(1,iCell)
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ call convert_lx(AvgXYZ(1,1), AvgXYZ(1,2), AvgXYZ(1,3), a_unit, latlon1)
+
+ iVert2 = verticesOnCell(2,iCell)
+ latlon2%lat = latVertex(iVert2)
+ latlon2%lon = lonVertex(iVert2)
+ call convert_lx(AvgXYZ(2,1), AvgXYZ(2,2), AvgXYZ(2,3), a_unit, latlon2)
+
+ iVert3 = verticesOnCell(3,iCell)
+ latlon3%lat = latVertex(iVert3)
+ latlon3%lon = lonVertex(iVert3)
+ call convert_lx(AvgXYZ(3,1), AvgXYZ(3,2), AvgXYZ(3,3), a_unit, latlon3)
+
+ iVert4 = verticesOnCell(4,iCell)
+ latlon4%lat = latVertex(iVert4)
+ latlon4%lon = lonVertex(iVert4)
+ call convert_lx(AvgXYZ(4,1), AvgXYZ(4,2), AvgXYZ(4,3), a_unit, latlon4)
+
+ AvgXYZ(5,:) = (AvgXYZ(1,:) + AvgXYZ(2,:) + AvgXYZ(3,:) + AvgXYZ(4,:))/FOUR
+
+ RadiusTest = sqrt(AvgXYZ(5,1)**2 + AvgXYZ(5,2)**2 + AvgXYZ(5,3)**2)
+ AvgXYZ(5,:) = AvgXYZ(5,:)*a_unit/RadiusTest
+
+ xCell(iCell) = AvgXYZ(5,1)
+ yCell(iCell) = AvgXYZ(5,2)
+ zCell(iCell) = AvgXYZ(5,3)
+
+! convert back to lat,lon
+ call convert_xl(AvgXYZ(5,1), AvgXYZ(5,2), AvgXYZ(5,3), latlon2)
+ latCell(iCell) = latlon2%lat
+ lonCell(iCell) = latlon2%lon
+
+ enddo
+
+ do iRow = 1, ny
+ do iCol = 1, nx
+ iCell = cellID_ij(iCol, iRow)
+ if(iCell.eq.0) cycle
+ indexToCellID(iCell) = iCell
+ angleEdge(edgesOnCell(1,iCell)) = pi / TWO
+ angleEdge(edgesOnCell(2,iCell)) = 0.0
+ end do
+ end do
+
+ write(6,*) maxval(xVertex), minval(xVertex)
+ write(6,*) maxval(yVertex), minval(yVertex)
+
+ do iEdge=1,nEdges
+ indexToEdgeID(iEdge) = iEdge
+ nEdgesOnEdge(iEdge) = 6
+
+!
+! calculate distance across edges connecting cell centers
+!
+
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
+
+ if(iCell1>0.and.iCell2>0) then
+ latlon1%lat = latCell(iCell1)
+ latlon1%lon = lonCell(iCell1)
+ iCell2 = cellsOnEdge(2,iEdge)
+ latlon2%lat = latCell(iCell2)
+ latlon2%lon = lonCell(iCell2)
+ distance = sphere_distance(latlon1, latlon2, a_unit)
+ dcEdge(iEdge) = distance
+ else
+ dcEdge(iEdge) = 1.0
+ endif
+
+ end do
+
+ do iVertex=1,nVertices
+ indexToVertexID(iVertex) = iVertex
+ enddo
+
+
+ areaCell = 0
+ areaTriangle = 0
+ kiteAreasOnVertex = 0
+
+ do iCell = 1, nCells
+
+ iVert1 = verticesOnCell(1,iCell)
+ iVert2 = verticesOnCell(2,iCell)
+ iVert3 = verticesOnCell(3,iCell)
+ iVert4 = verticesOnCell(4,iCell)
+
+ iEdge1 = edgesOnCell(1,iCell)
+ iEdge2 = edgesOnCell(2,iCell)
+ iEdge3 = edgesOnCell(3,iCell)
+ iEdge4 = edgesOnCell(4,iCell)
+
+! lower left of iCell, which is kiteArea(3,iVert1)
+ latlon1%lat = latVertex(iVert1)
+ latlon1%lon = lonVertex(iVert1)
+ latlon2%lat = latCell(iCell)
+ latlon2%lon = lonCell(iCell)
+ latlon3%lat = latEdge(iEdge1)
+ latlon3%lon = lonEdge(iEdge1)
+ AreaLeft = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ latlon3%lat = latEdge(iEdge2)
+ latlon3%lon = lonEdge(iEdge2)
+ AreaRight = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ kiteAreasOnVertex(3,iVert1) = AreaLeft + AreaRight
+
+! lower right of iCell, which is kiteArea(4,iVert2)
+ latlon1%lat = latVertex(iVert2)
+ latlon1%lon = lonVertex(iVert2)
+ latlon3%lat = latEdge(iEdge2)
+ latlon3%lon = lonEdge(iEdge2)
+ AreaLeft = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ latlon3%lat = latEdge(iEdge3)
+ latlon3%lon = lonEdge(iEdge3)
+ AreaRight = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ kiteAreasOnVertex(4,iVert2) = AreaLeft + AreaRight
+
+
+! upper right of iCell, which is kiteArea(1,iVert3)
+ latlon1%lat = latVertex(iVert3)
+ latlon1%lon = lonVertex(iVert3)
+ latlon3%lat = latEdge(iEdge3)
+ latlon3%lon = lonEdge(iEdge3)
+ AreaLeft = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ latlon3%lat = latEdge(iEdge4)
+ latlon3%lon = lonEdge(iEdge4)
+ AreaRight = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ kiteAreasOnVertex(1,iVert3) = AreaLeft + AreaRight
+
+
+! upper left of iCell, which is kiteArea(2,iVert4)
+ latlon1%lat = latVertex(iVert4)
+ latlon1%lon = lonVertex(iVert4)
+ latlon3%lat = latEdge(iEdge4)
+ latlon3%lon = lonEdge(iEdge4)
+ AreaLeft = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ latlon3%lat = latEdge(iEdge1)
+ latlon3%lon = lonEdge(iEdge1)
+ AreaRight = triangle_area(latlon1, latlon2, latlon3, a_unit)
+ kiteAreasOnVertex(2,iVert4) = AreaLeft + AreaRight
+
+
+ areaCell(iCell) = kiteAreasOnVertex(3,iVert1) + &
+ kiteAreasOnVertex(4,iVert2) + &
+ kiteAreasOnVertex(1,iVert3) + &
+ kiteAreasOnVertex(2,iVert4)
+
+ end do
+
+ do iVertex = 1, nVertices
+ areaTriangle(iVertex) = 0.0
+ do j=1,nVertexDegree
+ areaTriangle(iVertex) = areaTriangle(iVertex) + kiteAreasOnVertex(j,iVertex)
+ end do
+ enddo
+
+ write(6,*) ' areas', sum(areaCell), sum(areaTriangle)
+
+ ! fill in boundaryVertex
+ boundaryVertex = -1
+
+ ! fill in boundaryEdge
+ boundaryEdge = 0
+ do iEdge=1,nEdges
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
+ if(iCell1.eq.0.and.iCell2.eq.0) then
+ write(6,*) ' can not be so'
+ stop
+ endif
+
+ if(iCell1.eq.0.or.iCell2.eq.0) then
+ boundaryEdge(:,iEdge) = 1
+ nEdgesOnEdge(iEdge) = 0
+ if(cellsOnEdge(1,iEdge).eq.0) then
+ cellsOnEdge(1,iEdge)=cellsOnEdge(2,iEdge)
+ cellsOnEdge(2,iEdge)=0
+ endif
+ endif
+
+ enddo
+
+ write(6,*) 'rotating cellsOnVertex'
+ do iVertex=1,nVertices
+ do while (cellsOnVertex(1,iVertex).eq.0)
+ do j=2,nVertexDegree
+ cellsOnVertex(j-1,iVertex) = cellsOnVertex(j,iVertex)
+ kiteAreasOnVertex(j-1,iVertex) = kiteAreasOnVertex(j,iVertex)
+ enddo
+ cellsOnVertex(nVertexDegree,iVertex) = 0
+ kiteAreasOnVertex(nVertexDegree,iVertex) = 0.0
+ end do
+ enddo
+
+ do iVertex=1,nVertices
+ if(cellsOnVertex(1,iVertex).eq.0) then
+ write(6,*) cellsOnVertex(:,iVertex)
+ stop
+ endif
+ enddo
+
+ write(6,*) minval(cellsOnEdge(1,:)), maxval(cellsOnEdge(1,:))
+ write(6,*) minval(cellsOnEdge(2,:)), maxval(cellsOnEdge(2,:))
+
+ !
+ ! fill in initial conditions below
+ ! NOTE: these initial conditions will likely be removed
+ ! from the grid.nc files at some point (soon).
+ ! Initialize fields in grid
+ !
+
+ fEdge(:) = 0.0
+ fVertex(:) = 0.0
+ h_s(:) = 0.0
+ u(:,:,:) = 0.0
+ v(:,:,:) = 0.0
+ vh(:,:,:) = 0.0
+ circulation(:,:,:) = 0.0
+ vorticity(:,:,:) = 0.0
+ ke(:,:,:) = 0.0
+ tracers(:,:,:,:) = 0.0
+ h(:,:,:) = 1.0
+! do i=1,nCells
+! r = sqrt((xCell(i) - (nx/2)*dc)**2.0 + (yCell(i) - (ny/2)*dc)**2.0)
+! if (r < 10.0*dc) then
+! tracers(1,1,i,1) = (20.0 / 2.0) * (1.0 + cos(pi*r/(10.0*dc))) + 0.0
+! h(1,i,1) = 1.0 + 0.1*cos(pi*r/(20.0*dc))
+! else
+! tracers(1,1,i,1) = 0.0
+! h(1,i,1) = 1.0
+! end if
+! end do
+
+ MaxDepth = maxval(POP_Depth)
+ do iRow = 2, ny - 1
+ do iCol = 1, nx
+ iCell = cellID_ij(iCol,iRow)
+ if(iCell>0) then
+ if (trim(POP_DepthFileName) /= 'unknown_POP_DepthFileName') then
+ h_s(iCell) = MaxDepth - POP_Depth(iCol,iRow)
+ h(1,iCell,1) = POP_Depth(iCol,iRow)
+ else
+ h_s(iCell) = 1000.0
+ h(1,iCell,1) = h_s(iCell)
+ endif
+ endif
+ enddo
+ enddo
+
+ h_s = 0.0
+ h = 1000.0
+
+ lonPerturb0 = 200.*pi/180. ! note units are radians
+ latPerturb0 = 20.*pi/180.
+ scalePerturb = 0.1
+ heightPerturb = 1.0 ! in meters
+! heightPerturb = 0.0 ! in meters
+ do i=1,nCells
+ r = sqrt((latCell(i) - latPerturb0)**2.0 + (lonCell(i) - lonPerturb0)**2.0)
+ if (r < 10.0*scalePerturb) then
+ distance = heightPerturb*exp(-r*r/(scalePerturb*scalePerturb))
+ else
+ distance = 0.0
+ end if
+ h(1,i,1) = h(1,i,1) + distance
+ end do
+
+ fVertex = 0.0
+ do i = 1,nVertices
+ fVertex(i) = TWO*omega*sin(latVertex(i))
+ enddo
+
+ fEdge = 0.0
+ do i = 1,nEdges
+ fEdge(i) = TWO*omega*sin(latEdge(i))
+ enddo
+
+ j = 0
+ do i=1,nEdges
+ if(boundaryEdge(1,i).gt.0) j=j+1
+ enddo
+
+ xCell = xCell*xscale
+ yCell = yCell*xscale
+ zCell = zCell*xscale
+ xVertex = xVertex*xscale
+ yVertex = yVertex*xscale
+ zVertex = zVertex*xscale
+
+
+ call write_OpenDX( nCells, &
+ nVertices, &
+ nVertexDegree, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ areaCell )
+
+ xCell = xCell/xscale
+ yCell = yCell/xscale
+ zCell = zCell/xscale
+ xVertex = xVertex/xscale
+ yVertex = yVertex/xscale
+ zVertex = zVertex/xscale
+
+
+!
+! multiply by radius to convert from unit sphere to real earth
+!
+ write(6,*) maxval(xCell)
+ xCell = xCell*a
+ write(6,*) maxval(xCell)
+ yCell = yCell*a
+ zCell = zCell*a
+ xEdge = xEdge*a
+ yEdge = yEdge*a
+ zEdge = zEdge*a
+ xVertex = xVertex*a
+ yVertex = yVertex*a
+ zVertex = zVertex*a
+ dvEdge = dvEdge*a
+ dcEdge = dcEdge*a
+ areaCell = areaCell*a*a
+ areaTriangle = areaTriangle*a*a
+ kiteAreasOnVertex = kiteAreasOnVertex*a*a
+
+ write(6,*) ' summary of data '
+ write(6,10) ' grid res ', dc
+ write(6,11) ' domain ', nx, ny
+ write(6,10) ' thickness ', minval(h), maxval(h)
+ write(6,10) ' fVertex ', minval(fVertex), maxval(fVertex)
+ write(6,10) ' fEdge ', minval(fEdge), maxval(fEdge)
+ write(6,10) ' u_src ', minval(u_src), maxval(u_src)
+ write(6,10) ' u ', minval(u), maxval(u)
+ write(6,10) ' xCell ', minval(xCell), maxval(xCell)
+ write(6,10) ' yCell ', minval(yCell), maxval(yCell)
+ write(6,10 ) ' boundaryEdge',float(j), float(j)/float(nx*ny)
+ 10 format(a12,2e15.5)
+ 11 format(a12,2i5)
+
+
+ !
+ ! Write grid to grid.nc file
+ !
+ call write_netcdf_init( nCells, nEdges, nVertices, maxEdges, nVertLevels, nTracers, nVertexDegree)
+
+ call write_netcdf_fields( 1, &
+ latCell, lonCell, xCell, yCell, zCell, indexToCellID, &
+ latEdge, lonEdge, xEdge, yEdge, zEdge, indexToEdgeID, &
+ latVertex, lonVertex, xVertex, yVertex, zVertex, indexToVertexID, &
+ cellsOnEdge, &
+ nEdgesOnCell, &
+ nEdgesOnEdge, &
+ edgesOnCell, &
+ edgesOnEdge, &
+ weightsOnEdge, &
+ dvEdge, &
+ dcEdge, &
+ angleEdge, &
+ areaCell, &
+ areaTriangle, &
+ cellsOnCell, &
+ verticesOnCell, &
+ verticesOnEdge, &
+ edgesOnVertex, &
+ cellsOnVertex, &
+ kiteAreasOnVertex, &
+ fEdge, &
+ fVertex, &
+ h_s, &
+ boundaryEdge, &
+ boundaryVertex, &
+ u_src, &
+ u, &
+ v, &
+ h, &
+ vh, &
+ circulation, &
+ vorticity, &
+ ke, &
+ tracers &
+ )
+
+ call write_netcdf_finalize()
+
+
+ !
+ ! Write a graph.info file to be partitioned by kmetis
+ !
+ !
+ ! Write out a file compatible with metis for block decomposition
+ !
+ m=nEdges
+ do i=1,nCells
+ do j=1,nEdgesOnCell(i)
+ if(cellsOnCell(j,i).eq.0) m=m-1
+ enddo
+ enddo
+
+ open(42,file='graph.info',form='formatted')
+ write(42,*) nCells, m
+ do i=1,nCells
+ itmp = 0; k = 0;
+ do j=1,nEdgesOnCell(i)
+ if(cellsOnCell(j,i).gt.0) then
+ k=k+1; itmp(k)=cellsOnCell(j,i)
+ endif
+ enddo
+ write(42,'(1x,12i8)',advance='no') (itmp(m),m=1,k)
+ write(42,'(1x)')
+ end do
+ close(42)
+
+
+end program pop_to_mpas
+
+
+subroutine decompose_nproc(nproc, nprocx, nprocy)
+
+ implicit none
+
+ integer, intent(in) :: nproc
+ integer, intent(out) :: nprocx, nprocy
+
+ do nprocx=int(sqrt(real(nproc))),1,-1
+ nprocy = nproc / nprocx
+ if (nprocy == ceiling(real(nproc)/real(nprocx))) return
+ end do
+
+end subroutine decompose_nproc
+
+
+ subroutine write_OpenDX( nCells, &
+ nVertices, &
+ nVertexDegree, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ areaCell )
+
+ integer, intent(in) :: nCells, nVertices, nVertexDegree
+ real (kind=8), dimension(nCells), intent(in) :: xCell
+ real (kind=8), dimension(nCells), intent(in) :: yCell
+ real (kind=8), dimension(nCells), intent(in) :: zCell
+ real (kind=8), dimension(nVertices), intent(in) :: xVertex
+ real (kind=8), dimension(nVertices), intent(in) :: yVertex
+ real (kind=8), dimension(nVertices), intent(in) :: zVertex
+ integer, dimension(nCells), intent(in) :: nEdgesOnCell
+ integer, dimension(nVertexDegree,nCells), intent(in) :: verticesOnCell
+ real (kind=8), dimension(nCells), intent(in) :: areaCell
+
+ character(len=80) :: a, b, c, d, e, f
+ integer :: i, j, k, nVerticesTotal, iEdge, iLoop
+
+ write(6,*) 'xCell', minval(xCell), maxval(xCell)
+
+ 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,*) nVerticesTotal, nCells*nVertexDegree
+
+ open(unit=1,file='dx/pop.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file pop.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 pop.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 pop.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 pop.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 pop.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)
+
+ open(unit=10,file='dx/pop.area.data',form='formatted',status='unknown')
+ open(unit=11,file='dx/pop.face.data',form='formatted',status='unknown')
+ open(unit=12,file='dx/pop.loop.data',form='formatted',status='unknown')
+ open(unit=13,file='dx/pop.edge.data',form='formatted',status='unknown')
+ open(unit=14,file='dx/pop.position.data',form='formatted',status='unknown')
+
+ iLoop = 0
+ iEdge = 0
+ do i=1,nCells
+ write(10,20) areaCell(i)
+ write(11,21) i-1
+ write(12,21) iLoop
+ iLoop = iLoop + nEdgesOnCell(i)
+ do j=1,nEdgesOnCell(i)
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ k = verticesOnCell(j,i)
+ write(14,22) xVertex(k), yVertex(k), zVertex(k)
+ write(15,23) j,i,k,xVertex(k), yVertex(k), zVertex(k)
+ enddo
+ 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)
+
+
+ end subroutine write_OpenDX
</font>
</pre>