<p><b>ringler@lanl.gov</b> 2010-05-05 16:22:57 -0600 (Wed, 05 May 2010)</p><p><br>
pop_to_mpas_ChannelQuad.F can be used to create a grid.nc file for a channel, i.e. periodic in x, boundaries in y.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas_ChannelQuad.F
===================================================================
--- branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas_ChannelQuad.F (rev 0)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas_ChannelQuad.F 2010-05-05 22:22:57 UTC (rev 255)
@@ -0,0 +1,1357 @@
+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 = 20000.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 = 1.0 ! Radius of Earth (m)
+ real, parameter :: a_unit = 1.0 ! Radius of unit sphere
+ real, parameter :: xscale = 1.0
+
+ real (kind=8), parameter :: &
+ h_total_max = 2000.0, &
+ u_max = 0.0, &
+ u_src_max = 0.1, & ! max wind stress, N/m2
+ beta = 1.4e-11, &
+ f0 = -1.1e-4
+
+ real (kind=8) :: ymid, ytmp, ymax, xmid, xloc, yloc, pert, ymin
+ real (kind=8), allocatable, dimension(:,:,:) :: utmp
+
+ 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, rho
+ 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) :: stress, distance, n1, n2 &
+ ,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_x, POP_y
+ type (geo_point) :: latlon1, latlon2, latlon3, latlon4, latlon5
+ integer :: recLength, ip1, im1
+
+!
+! begin
+!
+ call cell_indexing_read_nl()
+
+ allocate(utmp(0:nx+1, 0:ny+1, 2))
+ utmp = 0.0
+
+ allocate(POP_uLat(nx,ny), POP_uLon(nx,ny), POP_Angle(nx,ny))
+ allocate(POP_x(nx,ny), POP_y(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)
+ POP_uLat = 0; POP_uLon = 0; POP_Angle = 0
+ write(6,*) 'lat/lon',minval(POP_uLat),maxval(POP_uLat),minval(POP_uLon),maxval(POP_uLon)
+
+ do i=1,nx
+ do j=1,ny
+ POP_x(i,j) = dc*i
+ POP_y(i,j) = dc*j
+ enddo
+ enddo
+ write(6,*) 'POP_x,y', minval(POP_x), maxval(POP_x), minval(POP_y), maxval(POP_y)
+
+!
+! 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)
+
+ POP_KMT = 0
+ POP_KMT(1:nx,2:ny-1) = nVertLevels
+
+ 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(rho(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)
+
+ xVertex(iVert1) = POP_x(im1 , iRow-1)
+ yVertex(iVert1) = POP_y(im1 , iRow-1)
+ zVertex(iVert1) = 0.0
+ xVertex(iVert2) = POP_x(iCol, iRow-1)
+ yVertex(iVert2) = POP_y(iCol, iRow-1)
+ zVertex(iVert1) = 0.0
+ xVertex(iVert3) = POP_x(iCol, iRow )
+ yVertex(iVert3) = POP_y(iCol, iRow )
+ zVertex(iVert1) = 0.0
+ xVertex(iVert4) = POP_x(im1 , iRow )
+ yVertex(iVert4) = POP_y(im1 , iRow )
+ zVertex(iVert1) = 0.0
+
+ 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,*)
+
+
+!
+! average vertex coordinates to create edge coordinates
+!
+ do iCell = 1, nCells
+!
+! edge 1
+!
+ do iEdge = 1,4
+
+ iVert1 = verticesOnEdge(1,edgesOnCell(iEdge,iCell))
+ iVert2 = verticesOnEdge(2,edgesOnCell(iEdge,iCell))
+
+ AvgXYZ(1,1) = xVertex(iVert1)
+ AvgXYZ(1,2) = yVertex(iVert1)
+ AvgXYZ(1,3) = zVertex(iVert1)
+
+ AvgXYZ(2,1) = xVertex(iVert2)
+ AvgXYZ(2,2) = yVertex(iVert2)
+ AvgXYZ(2,3) = zVertex(iVert2)
+
+ AvgXYZ(3,:) = (AvgXYZ(1,:) + AvgXYZ(2,:))/TWO
+
+ xEdge(edgesOnCell(iEdge,iCell)) = AvgXYZ(3,1)
+ yEdge(edgesOnCell(iEdge,iCell)) = AvgXYZ(3,2)
+ zEdge(edgesOnCell(iEdge,iCell)) = AvgXYZ(3,3)
+
+! calculate distance between vertices
+ dvEdge(edgesOnCell(iEdge,iCell)) = dc
+
+! convert back to lat,lon
+ latEdge(edgesOnCell(iEdge,iCell)) = 0.0
+ lonEdge(edgesOnCell(iEdge,iCell)) = 0.0
+
+ enddo
+ enddo
+
+!
+! average vertex coordinates to create centered coordinates
+!
+ do iCell = 1, nCells
+
+ iVert1 = verticesOnCell(1,iCell)
+ iVert2 = verticesOnCell(2,iCell)
+ iVert3 = verticesOnCell(3,iCell)
+ iVert4 = verticesOnCell(4,iCell)
+
+ AvgXYZ(1,1) = xVertex(iVert1)
+ AvgXYZ(1,2) = yVertex(iVert1)
+ AvgXYZ(1,3) = zVertex(iVert1)
+
+ AvgXYZ(2,1) = xVertex(iVert2)
+ AvgXYZ(2,2) = yVertex(iVert2)
+ AvgXYZ(2,3) = zVertex(iVert2)
+
+ AvgXYZ(3,1) = xVertex(iVert3)
+ AvgXYZ(3,2) = yVertex(iVert3)
+ AvgXYZ(3,3) = zVertex(iVert3)
+
+ AvgXYZ(4,1) = xVertex(iVert4)
+ AvgXYZ(4,2) = yVertex(iVert4)
+ AvgXYZ(4,3) = zVertex(iVert4)
+
+ AvgXYZ(5,:) = (AvgXYZ(1,:) + AvgXYZ(2,:) + AvgXYZ(3,:) + AvgXYZ(4,:))/FOUR
+
+ xCell(iCell) = AvgXYZ(5,1)
+ yCell(iCell) = AvgXYZ(5,2)
+ zCell(iCell) = AvgXYZ(5,3)
+
+ latCell(iCell) = 0.0
+ lonCell(iCell) = 0.0
+
+ 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,*) ' vertex positions'
+ 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
+ dcEdge(iEdge) = dc
+ 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)
+ kiteAreasOnVertex(3,iVert1) = dc**2/4
+
+! lower right of iCell, which is kiteArea(4,iVert2)
+ kiteAreasOnVertex(4,iVert2) = dc**2/4
+
+! upper right of iCell, which is kiteArea(1,iVert3)
+ kiteAreasOnVertex(1,iVert3) = dc**2/4
+
+! upper left of iCell, which is kiteArea(2,iVert4)
+ kiteAreasOnVertex(2,iVert4) = dc**2/4
+
+ 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,*) 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
+
+ ! setting for three levels
+ h(1,:,:) = 350.0
+ h(2,:,:) = 750.0
+ h(3,:,:) = 1000.0
+ h_s(:) = -( h(1,:,1) + h(2,:,1) + h(3,:,1) )
+
+ rho(1,:,:) = 1000.0
+ rho(2,:,:) = 1003.0
+ rho(3,:,:) = 1008.0
+
+ ! setting for one level
+ ! h(1,:,:) = 250.0
+ ! h_s(:) = -h(1,:,1)
+ ! rho(1,:,:) = 1000.0
+
+ ymid = dc*real(ny/2)
+
+ fVertex = 0.0
+ do i = 1,nVertices
+ fVertex(i) = f0 + (yVertex(i) - ymid) * beta
+ enddo
+
+ fEdge = 0.0
+ do i = 1,nEdges
+ fEdge(i) = f0 + (yEdge(i) - ymid) * beta
+ enddo
+
+ u_src = 0.0
+ do i = 1,nEdges
+ ymin = minval(yEdge)
+ ymax = maxval(yEdge)
+ ymax = ymax - ymin
+ ytmp = yEdge(i)
+ iCell1 = cellsOnEdge(1,i)
+ iCell2 = cellsOnEdge(2,i)
+ if(iCell1>0.and.iCell2>0) then
+ pert = u_src_max * sin ( pi * ytmp / ymax)
+ n1 = xCell(iCell2) - xCell(iCell1)
+ n2 = yCell(iCell2) - yCell(iCell1)
+ distance = sqrt(n1**2+n2**2)
+ n1 = n1 / distance
+ n2 = n2 / distance
+ u_src(1,i) = pert * n1
+
+! hack fix
+ u_src(1,i) = abs(u_src(1,i))
+
+ write(6,*) i,n1,n2, u_src(1,i)
+ endif
+ 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
+ xEdge = xEdge*xscale
+ yEdge = yEdge*xscale
+ zEdge = zEdge*xscale
+
+ call write_OpenDX( nCells, &
+ nVertices, &
+ nEdges, &
+ nVertexDegree, &
+ maxEdges, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ edgesOnCell, &
+ areaCell, &
+ kiteAreasOnVertex )
+
+ xCell = xCell/xscale
+ yCell = yCell/xscale
+ zCell = zCell/xscale
+ xVertex = xVertex/xscale
+ yVertex = yVertex/xscale
+ zVertex = zVertex/xscale
+ xEdge = xEdge/xscale
+ yEdge = yEdge/xscale
+ zEdge = zEdge/xscale
+
+ write(6,*) 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, &
+ rho, &
+ 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, &
+ nEdges, &
+ nVertexDegree, &
+ maxEdges, &
+ xCell, &
+ yCell, &
+ zCell, &
+ xVertex, &
+ yVertex, &
+ zVertex, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
+ nEdgesOnCell, &
+ verticesOnCell, &
+ edgesOnCell, &
+ areaCell, &
+ kiteAreasOnVertex )
+
+ implicit none
+
+ integer, intent(in) :: nCells, nVertices, nVertexDegree, nEdges, maxEdges
+ 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
+ real (kind=8), dimension(nEdges), intent(in) :: xEdge
+ real (kind=8), dimension(nEdges), intent(in) :: yEdge
+ real (kind=8), dimension(nEdges), intent(in) :: zEdge
+ integer, dimension(nCells), intent(in) :: nEdgesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: edgesOnCell
+ real (kind=8), dimension(nCells), intent(in) :: areaCell
+ real (kind=8), dimension(nVertexDegree,nVertices), intent(in) :: kiteAreasOnVertex
+
+ character(len=80) :: a, b, c, d, e, f
+ integer :: i, j, k, nVerticesTotal, iEdge, iLoop, iFace, Vert(4), Edge(4), iVertex, i1, i2
+ real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4
+
+ 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*4
+
+ 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)
+ if(k.le.0) write(6,*) ' vert1 ',k, verticesOnCell(:,i)
+ write(14,22) xVertex(k), yVertex(k), zVertex(k)
+ write(15,23) j,i,k,xVertex(k), yVertex(k), zVertex(k)
+ enddo
+ enddo
+
+ 20 format(e20.10)
+ 21 format(i20)
+ 22 format(3e20.10)
+ 23 format(3i6, 3e20.10)
+
+ close(10)
+ close(11)
+ close(12)
+ close(13)
+ close(14)
+
+
+
+ nVerticesTotal = 0
+ do i=1,nCells
+ nVerticesTotal = nVerticesTotal + nEdgesOnCell(i)**2
+ enddo
+ write(6,*) nVerticesTotal, nCells*nVertexDegree**2
+
+ open(unit=1,file='dx/kite.dx',form='formatted',status='unknown')
+
+ a = trim('object "positions list" class array type float rank 1 shape 3 items')
+ b = trim('ascii data file kite.position.data')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,*)
+
+ a = trim('object "edge list" class array type int rank 0 items')
+ b = trim('ascii data file kite.edge.data')
+ c = trim('attribute "ref" string "positions"')
+ write(1,10) a, nVerticesTotal
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "loops list" class array type int rank 0 items')
+ b = trim('ascii data file kite.loop.data')
+ c = trim('attribute "ref" string "edges"')
+ write(1,10) a, nCells*nVertexDegree
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "face list" class array type int rank 0 items')
+ b = trim('ascii data file kite.face.data')
+ c = trim('attribute "ref" string "loops"')
+ write(1,10) a, nCells*nVertexDegree
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object 0 class array type float rank 0 items')
+ b = trim('data file kite.area.data')
+ c = trim('attribute "dep" string "faces"')
+ write(1,10) a, nCells*nVertexDegree
+ write(1,10) b
+ write(1,10) c
+ write(1,*)
+
+ a = trim('object "area" class field')
+ b = trim('component "positions" "positions list"')
+ c = trim('component "edges" "edge list"')
+ d = trim('component "loops" "loops list"')
+ e = trim('component "faces" "face list"')
+ f = trim('component "data" 0')
+ write(1,10) a
+ write(1,10) b
+ write(1,10) c
+ write(1,10) d
+ write(1,10) e
+ write(1,10) f
+
+ close(1)
+
+ open(unit=10,file='dx/kite.area.data',form='formatted',status='unknown')
+ open(unit=11,file='dx/kite.face.data',form='formatted',status='unknown')
+ open(unit=12,file='dx/kite.loop.data',form='formatted',status='unknown')
+ open(unit=13,file='dx/kite.edge.data',form='formatted',status='unknown')
+ open(unit=14,file='dx/kite.position.data',form='formatted',status='unknown')
+
+ iLoop = 0
+ iEdge = 0
+ iFace = 0
+
+ do i=1,nCells
+
+ Vert(1) = verticesOnCell(1,i)
+ Vert(2) = verticesOnCell(2,i)
+ Vert(3) = verticesOnCell(3,i)
+ Vert(4) = verticesOnCell(4,i)
+
+ if(minval(Vert).le.0) then
+ write(6,*) 'vert ',i, Vert
+ endif
+
+ Edge(1) = edgesOnCell(1,i)
+ Edge(2) = edgesOnCell(2,i)
+ Edge(3) = edgesOnCell(3,i)
+ Edge(4) = edgesOnCell(4,i)
+
+ if(minval(Edge).le.0) then
+ write(6,*) 'edge ', i, Edge
+ endif
+
+ write(10,20) kiteAreasOnVertex(3,Vert(1))
+ write(10,20) kiteAreasOnVertex(4,Vert(2))
+ write(10,20) kiteAreasOnVertex(1,Vert(3))
+ write(10,20) kiteAreasOnVertex(2,Vert(4))
+
+ if(kiteAreasOnVertex(3,Vert(1)).lt.1.0e-10) then
+ write(6,*) '3 Vert1', Vert(1), i,kiteAreasOnVertex(3,Vert(1))
+ stop
+ endif
+
+ if(kiteAreasOnVertex(4,Vert(2)).lt.1.0e-10) then
+ write(6,*) '4 Vert2'
+ stop
+ endif
+
+ if(kiteAreasOnVertex(1,Vert(3)).lt.1.0e-10) then
+ write(6,*) '1 Vert3'
+ stop
+ endif
+
+ if(kiteAreasOnVertex(2,Vert(4)).lt.1.0e-10) then
+ write(6,*) '2 Vert4'
+ stop
+ endif
+
+
+ do iVertex=1,4
+ write(11,21) iFace
+ write(12,21) iLoop
+ iFace = iFace + 1
+ iLoop = iLoop + 4
+ do j=1,4
+ write(13,21) iEdge
+ iEdge = iEdge + 1
+ enddo
+ enddo
+
+ i1 = 1; i2 = 2
+ x1 = xCell(i); y1 = yCell(i); z1 = zCell(i)
+ x2 = xEdge(Edge(i1)); y2 = yEdge(Edge(i1)); z2 = zEdge(Edge(i1))
+ x3 = xVertex(Vert(i1)); y3 = yVertex(Vert(i1)); z3 = zVertex(Vert(i1))
+ x4 = xEdge(Edge(i2)); y4 = yEdge(Edge(i2)); z4 = zEdge(Edge(i2))
+ write(14,22) x1, y1, z1
+ write(14,22) x2, y2, z2
+ write(14,22) x3, y3, z3
+ write(14,22) x4, y4, z4
+
+ i1 = 2; i2 = 3
+ x1 = xCell(i); y1 = yCell(i); z1 = zCell(i)
+ x2 = xEdge(Edge(i1)); y2 = yEdge(Edge(i1)); z2 = zEdge(Edge(i1))
+ x3 = xVertex(Vert(i1)); y3 = yVertex(Vert(i1)); z3 = zVertex(Vert(i1))
+ x4 = xEdge(Edge(i2)); y4 = yEdge(Edge(i2)); z4 = zEdge(Edge(i2))
+ write(14,22) x1, y1, z1
+ write(14,22) x2, y2, z2
+ write(14,22) x3, y3, z3
+ write(14,22) x4, y4, z4
+
+ i1 = 3; i2 = 4
+ x1 = xCell(i); y1 = yCell(i); z1 = zCell(i)
+ x2 = xEdge(Edge(i1)); y2 = yEdge(Edge(i1)); z2 = zEdge(Edge(i1))
+ x3 = xVertex(Vert(i1)); y3 = yVertex(Vert(i1)); z3 = zVertex(Vert(i1))
+ x4 = xEdge(Edge(i2)); y4 = yEdge(Edge(i2)); z4 = zEdge(Edge(i2))
+ write(14,22) x1, y1, z1
+ write(14,22) x2, y2, z2
+ write(14,22) x3, y3, z3
+ write(14,22) x4, y4, z4
+
+ i1 = 4; i2 = 1
+ x1 = xCell(i); y1 = yCell(i); z1 = zCell(i)
+ x2 = xEdge(Edge(i1)); y2 = yEdge(Edge(i1)); z2 = zEdge(Edge(i1))
+ x3 = xVertex(Vert(i1)); y3 = yVertex(Vert(i1)); z3 = zVertex(Vert(i1))
+ x4 = xEdge(Edge(i2)); y4 = yEdge(Edge(i2)); z4 = zEdge(Edge(i2))
+ write(14,22) x1, y1, z1
+ write(14,22) x2, y2, z2
+ write(14,22) x3, y3, z3
+ write(14,22) x4, y4, z4
+
+
+ enddo
+
+ close(10)
+ close(11)
+ close(12)
+ close(13)
+ close(14)
+
+
+ end subroutine write_OpenDX
</font>
</pre>