<p><b>ringler@lanl.gov</b> 2010-04-06 15:53:34 -0600 (Tue, 06 Apr 2010)</p><p><br>
code changed: pop_to_mpas.F<br>
<br>
reason: added support to visualize the sub-areas within a single cell<br>
</p><hr noshade><pre><font color="gray">Modified: 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 2010-04-06 16:29:36 UTC (rev 181)
+++ branches/ocean_projects/grid_gen_ocean/pop_to_mpas/pop_to_mpas.F 2010-04-06 21:53:34 UTC (rev 182)
@@ -277,50 +277,50 @@
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
+ ! 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,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
+ ! ! 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))
@@ -838,7 +838,27 @@
AreaRight = triangle_area(latlon1, latlon2, latlon3, a_unit)
kiteAreasOnVertex(2,iVert4) = AreaLeft + AreaRight
+ if(kiteAreasOnVertex(3,iVert1).lt.1.0e-10) then
+ write(6,*) '3 iVert1'
+ stop
+ endif
+ if(kiteAreasOnVertex(4,iVert2).lt.1.0e-10) then
+ write(6,*) '4 iVert2'
+ stop
+ endif
+
+ if(kiteAreasOnVertex(1,iVert3).lt.1.0e-10) then
+ write(6,*) '1 iVert3'
+ stop
+ endif
+
+ if(kiteAreasOnVertex(2,iVert4).lt.1.0e-10) then
+ write(6,*) '2 iVert4'
+ stop
+ endif
+
+
areaCell(iCell) = kiteAreasOnVertex(3,iVert1) + &
kiteAreasOnVertex(4,iVert2) + &
kiteAreasOnVertex(1,iVert3) + &
@@ -879,25 +899,26 @@
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
+ ! 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
+ ! 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,:))
@@ -985,20 +1006,31 @@
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, &
- areaCell )
+ edgesOnCell, &
+ areaCell, &
+ kiteAreasOnVertex )
xCell = xCell/xscale
yCell = yCell/xscale
@@ -1006,8 +1038,34 @@
xVertex = xVertex/xscale
yVertex = yVertex/xscale
zVertex = zVertex/xscale
+ xEdge = xEdge/xscale
+ yEdge = yEdge/xscale
+ zEdge = zEdge/xscale
+ 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,:))
+
+
+
!
! multiply by radius to convert from unit sphere to real earth
!
@@ -1136,30 +1194,45 @@
subroutine write_OpenDX( nCells, &
nVertices, &
+ nEdges, &
nVertexDegree, &
+ maxEdges, &
xCell, &
yCell, &
zCell, &
xVertex, &
yVertex, &
zVertex, &
+ xEdge, &
+ yEdge, &
+ zEdge, &
nEdgesOnCell, &
verticesOnCell, &
- areaCell )
+ edgesOnCell, &
+ areaCell, &
+ kiteAreasOnVertex )
- integer, intent(in) :: nCells, nVertices, nVertexDegree
+ 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(nVertexDegree,nCells), intent(in) :: verticesOnCell
+ 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
+ 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)
@@ -1200,7 +1273,7 @@
do i=1,nCells
nVerticesTotal = nVerticesTotal + nEdgesOnCell(i)
enddo
- write(6,*) nVerticesTotal, nCells*nVertexDegree
+ write(6,*) nVerticesTotal, nCells*4
open(unit=1,file='dx/pop.dx',form='formatted',status='unknown')
@@ -1275,6 +1348,7 @@
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
@@ -1292,4 +1366,183 @@
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>