<p><b>ringler@lanl.gov</b> 2011-01-28 12:02:48 -0700 (Fri, 28 Jan 2011)</p><p><br>
This modification brings the specification of limited-area domains in line with previous design decisions.<br>
<br>
Notably, the MPAS framework assumes throughout that the positive normal direction points from cell1 to cell2, where<br>
the global index of cell1 is less than the global index of cell2. In that, we define the positive normal direction to point<br>
from the lower cell global ID to the higer cell global ID. We then use this definition to define the tangent direction as<br>
k cross n where k is the local vertical and n is the normal vector. From a vantage point at cell1 looking in the positive<br>
normal direction (i.e. toward cell2) the positive tangent direction is from right to left.<br>
<br>
basin.F takes a global domain and produces a limited area domain. The rule noted immediately above was not applied on the<br>
boundary edges. At a boundary edge, the grid cell outside of the limited-area domain is culled leaving only one real cell<br>
associated with that boundary edge. In the MPAS dynamical core, the memory location of this culled cell points to nCells+1, i.e. the last entry in the list. In keeping with the definition that the positive points from lower to higher, we chose to<br>
define the positive normal to always point out of the domain, i.e. from the interior (real cell) toward the ficticious cell<br>
at nCells+1.<br>
<br>
This leads to limited area grid where along the boundaryEdges the positive normal always points out of the domain and the<br>
tangent direction along the bounaryEdges always loops in the CCW direction. This property might of use when deriving lateral boundary velocity forcing.<br>
<br>
The basin.F code now follows this definition. It does by based on the assumption that the global mesh that is used as input<br>
follows this definition.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin/src/basin.F
===================================================================
--- branches/ocean_projects/basin/src/basin.F 2011-01-28 14:58:26 UTC (rev 708)
+++ branches/ocean_projects/basin/src/basin.F 2011-01-28 19:02:48 UTC (rev 709)
@@ -67,7 +67,8 @@
! Step 2: Set if the grid is on a sphere or not, and it's radius
character (len=16) :: on_a_sphere = 'YES '
-real*8, parameter :: sphere_radius = 6.37122e6
+!real*8, parameter :: sphere_radius = 6.37122e6
+real*8, parameter :: sphere_radius = 1.0
!character (len=16) :: on_a_sphere = 'NO '
!real*8, parameter :: sphere_radius = 0.0
@@ -96,7 +97,7 @@
real, allocatable, dimension(:) :: xCellNew, yCellNew, zCellNew, latCellNew, lonCellNew
real, allocatable, dimension(:) :: xEdgeNew, yEdgeNew, zEdgeNew, latEdgeNew, lonEdgeNew
real, allocatable, dimension(:) :: xVertexNew, yVertexNew, zVertexNew, latVertexNew, lonVertexNew
-integer, allocatable, dimension(:) :: nEdgesOnCellNew, nEdgesOnEdgeNew
+integer, allocatable, dimension(:) :: nEdgesOnCellNew, nEdgesOnEdgeNew, flipVerticesOnEdgeOrdering
integer, allocatable, dimension(:,:) :: cellsOnCellNew, edgesOnCellNew, verticesOnCellNew
integer, allocatable, dimension(:,:) :: cellsOnEdgeNew, verticesOnEdgeNew, edgesOnEdgeNew
integer, allocatable, dimension(:,:) :: cellsOnVertexNew, edgesOnVertexNew
@@ -994,7 +995,7 @@
real (kind=4), allocatable, dimension(:,:) :: ztopo
integer :: nx, ny, inx, iny, ix, iy
real :: pi, dtr, zdata, rlon, rlat, r, ymin, ymax
-real :: latmin, latmax;
+real :: latmin, latmax, lonmin, lonmax
logical :: flag, kmt_flag
pi = 4.0*atan(1.0)
dtr = pi / 180.0
@@ -1006,12 +1007,16 @@
kmt = nVertLevelsMOD
if(on_a_sphere.eq.'YES ') then
write(6,*) 'Working on a sphere'
- latmin = -70*dtr
- latmax = -25*dtr
+ latmin = -30*dtr
+ latmax = +30*dtr
+ lonmin = +10*dtr
+ lonmax = +70*dtr
write(6,*) ' lat min ', latmin
write(6,*) ' lat max ', latmax
where(latCell.lt.latmin) kmt = 0
where(latCell.gt.latmax) kmt = 0
+ where(lonCell.lt.lonmin) kmt = 0
+ where(lonCell.gt.lonmax) kmt = 0
else
ymin = minval(yCell)
write(6,*) ' minimum yCell ', ymin
@@ -1288,8 +1293,10 @@
allocate(cellsOnEdgeNew(TWONew,nEdgesNew))
allocate(boundaryEdgeNew(nVertLevelsNew,nEdgesNew))
+allocate(flipVerticesOnEdgeOrdering(nEdgesNew))
cellsOnEdgeNew(:,:) = 0
boundaryEdgeNew(:,:) = 0
+flipVerticesOnEdgeOrdering(:) = 0
do iEdge=1,nEdges
if(edgeMap(iEdge).eq.0) cycle
iEdgeNew = edgeMap(iEdge)
@@ -1304,201 +1311,207 @@
if(iCell1New.eq.0) then
cellsOnEdgeNew(1,iEdgeNew) = iCell2New
cellsOnEdgeNew(2,iEdgeNew) = iCell1New
+ flipVerticesOnEdgeOrdering(iEdgeNew) = 1
endif
enddo
allocate(verticesOnEdgeNew(TWONew,nEdgesNew))
- allocate(boundaryVertexNew(nVertLevelsNew,nVerticesNew))
- verticesOnEdgeNew(:,:) = 0
- boundaryVertexNew(:,:) = 0
- do iEdge=1,nEdges
- if(edgeMap(iEdge).eq.0) cycle
- iEdgeNew = edgeMap(iEdge)
- iVertex1 = VerticesOnEdge(1,iEdge)
- iVertex2 = VerticesOnEdge(2,iEdge)
- iVertex1New = vertexMap(iVertex1)
- iVertex2New = vertexMap(iVertex2)
- if(iVertex1New.eq.0.or.iVertex2New.eq.0) stop "verticesOnEdge"
- verticesOnEdgeNew(1,iEdgeNew) = iVertex1New
- verticesOnEdgeNew(2,iEdgeNew) = iVertex2New
- if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
- boundaryVertexNew(:,iVertex1New)=1
- boundaryVertexNew(:,iVertex2New)=1
- endif
- enddo
+allocate(boundaryVertexNew(nVertLevelsNew,nVerticesNew))
+verticesOnEdgeNew(:,:) = 0
+boundaryVertexNew(:,:) = 0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+iVertex1 = VerticesOnEdge(1,iEdge)
+iVertex2 = VerticesOnEdge(2,iEdge)
+iVertex1New = vertexMap(iVertex1)
+iVertex2New = vertexMap(iVertex2)
+if(iVertex1New.eq.0.or.iVertex2New.eq.0) stop "verticesOnEdge"
+if(flipVerticesOnEdgeOrdering(iEdgeNew).eq.0) then
+ verticesOnEdgeNew(1,iEdgeNew) = iVertex1New
+ verticesOnEdgeNew(2,iEdgeNew) = iVertex2New
+else
+ verticesOnEdgeNew(1,iEdgeNew) = iVertex2New
+ verticesOnEdgeNew(2,iEdgeNew) = iVertex1New
+endif
+if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
+ boundaryVertexNew(:,iVertex1New)=1
+ boundaryVertexNew(:,iVertex2New)=1
+endif
+enddo
- allocate(nEdgesOnEdgeNew(nEdgesNew))
- allocate(edgesOnEdgeNew(maxEdges2,nEdgesNew))
- allocate(weightsOnEdgeNew(maxEdges2,nEdgesNew))
- nEdgesOnEdgeNew(:) = 0
- edgesOnEdgeNew(:,:) = 0
- weightsOnEdgeNew(:,:) = 0.0
- do iEdge=1,nEdges
- if(edgeMap(iEdge).eq.0) cycle
- iEdgeNew = edgeMap(iEdge)
- if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
- nEdgesOnEdgeNew(iEdgeNew) = 0
- edgesOnEdgeNew(:,iEdgeNew) = 0
- weightsOnEdgeNew(:,iEdgeNew) = 0.0
- else
- nEdgesOnEdgeNew(iEdgeNew) = nEdgesOnEdge(iEdge)
- do i=1,nEdgesOnEdgeNew(iEdgeNew)
- jEdge = edgesOnEdge(i,iEdge)
- jEdgeNew = edgeMap(jEdge)
- if(jEdgeNew.eq.0) stop "jEdgeNew"
- edgesOnEdgeNew(i,iEdgeNew)=jEdgeNew
- weightsOnEdgeNew(i,iEdgeNew) = weightsOnEdge(i,iEdge)
- enddo
- endif
- enddo
+allocate(nEdgesOnEdgeNew(nEdgesNew))
+allocate(edgesOnEdgeNew(maxEdges2,nEdgesNew))
+allocate(weightsOnEdgeNew(maxEdges2,nEdgesNew))
+nEdgesOnEdgeNew(:) = 0
+edgesOnEdgeNew(:,:) = 0
+weightsOnEdgeNew(:,:) = 0.0
+do iEdge=1,nEdges
+if(edgeMap(iEdge).eq.0) cycle
+iEdgeNew = edgeMap(iEdge)
+if(boundaryEdgeNew(1,iEdgeNew).eq.1) then
+ nEdgesOnEdgeNew(iEdgeNew) = 0
+ edgesOnEdgeNew(:,iEdgeNew) = 0
+ weightsOnEdgeNew(:,iEdgeNew) = 0.0
+else
+ nEdgesOnEdgeNew(iEdgeNew) = nEdgesOnEdge(iEdge)
+ do i=1,nEdgesOnEdgeNew(iEdgeNew)
+ jEdge = edgesOnEdge(i,iEdge)
+ jEdgeNew = edgeMap(jEdge)
+ if(jEdgeNew.eq.0) stop "jEdgeNew"
+ edgesOnEdgeNew(i,iEdgeNew)=jEdgeNew
+ weightsOnEdgeNew(i,iEdgeNew) = weightsOnEdge(i,iEdge)
+ enddo
+endif
+enddo
- allocate(cellsOnCellNew(maxEdges,nCellsNew))
- allocate(nEdgesOnCellNew(nCellsNew))
- cellsOnCellNew = 0
- nEdgesOnCellNew = 0
- do iCell=1,nCells
- if(cellMap(iCell).eq.0) cycle
- iCellNew = cellMap(iCell)
- nEdgesOnCellNew(iCellNew)=nEdgesOnCell(iCell)
- do i=1,nEdgesOnCellNew(iCellNew)
- j = cellsOnCell(i,iCell)
- jNew = cellMap(j)
- cellsOnCellNew(i,iCellNew) = jNew
- enddo
- enddo
+allocate(cellsOnCellNew(maxEdges,nCellsNew))
+allocate(nEdgesOnCellNew(nCellsNew))
+cellsOnCellNew = 0
+nEdgesOnCellNew = 0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+nEdgesOnCellNew(iCellNew)=nEdgesOnCell(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = cellsOnCell(i,iCell)
+jNew = cellMap(j)
+cellsOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
- allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
- edgesOnCellNew(:,:) = 0
- do iCell=1,nCells
- if(cellMap(iCell).eq.0) cycle
- iCellNew = cellMap(iCell)
- do i=1,nEdgesOnCellNew(iCellNew)
- j = edgesOnCell(i,iCell)
- jNew = edgeMap(j)
- if(jNew.eq.0) stop "edgesOnCell"
- edgesOnCellNew(i,iCellNew) = jNew
- enddo
- enddo
+allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
+edgesOnCellNew(:,:) = 0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j = edgesOnCell(i,iCell)
+jNew = edgeMap(j)
+if(jNew.eq.0) stop "edgesOnCell"
+edgesOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
- allocate(verticesOnCellNew(maxEdgesNew,nCellsNew))
- verticesOnCellNew(:,:)=0
- do iCell=1,nCells
- if(cellMap(iCell).eq.0) cycle
- iCellNew = cellMap(iCell)
- do i=1,nEdgesOnCellNew(iCellNew)
- j=verticesOnCell(i,iCell)
- jNew = vertexMap(j)
- if(jNew.eq.0) stop "verticesOnCell"
- verticesOnCellNew(i,iCellNew) = jNew
- enddo
- enddo
+allocate(verticesOnCellNew(maxEdgesNew,nCellsNew))
+verticesOnCellNew(:,:)=0
+do iCell=1,nCells
+if(cellMap(iCell).eq.0) cycle
+iCellNew = cellMap(iCell)
+do i=1,nEdgesOnCellNew(iCellNew)
+j=verticesOnCell(i,iCell)
+jNew = vertexMap(j)
+if(jNew.eq.0) stop "verticesOnCell"
+verticesOnCellNew(i,iCellNew) = jNew
+enddo
+enddo
- allocate(cellsOnVertexNew(vertexDegreeNew,nVerticesNew))
- allocate(kiteAreasOnVertexNew(vertexDegreeNew,nVerticesNew))
- cellsOnVertexNew = 0
- kiteAreasOnVertexNew = 0
- do iVertex=1,nVertices
- if(vertexMap(iVertex).eq.0) cycle
- iVertexNew = vertexMap(iVertex)
- do i=1,vertexDegree
- j=cellsOnVertex(i,iVertex)
- jNew=cellMap(j)
- if(jNew.eq.0) then
- kiteAreasOnVertexNew(i,iVertexNew)=0
- else
- kiteAreasOnVertexNew(i,iVertexNew)=kiteAreasOnVertex(i,iVertex)
- endif
- cellsOnVertexNew(i,iVertexNew)=jNew
- enddo
- enddo
+allocate(cellsOnVertexNew(vertexDegreeNew,nVerticesNew))
+allocate(kiteAreasOnVertexNew(vertexDegreeNew,nVerticesNew))
+cellsOnVertexNew = 0
+kiteAreasOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=cellsOnVertex(i,iVertex)
+jNew=cellMap(j)
+if(jNew.eq.0) then
+ kiteAreasOnVertexNew(i,iVertexNew)=0
+else
+ kiteAreasOnVertexNew(i,iVertexNew)=kiteAreasOnVertex(i,iVertex)
+endif
+cellsOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
- areaTriangleNew = 0
- do iVertex=1,nVerticesNew
- do i=1,vertexDegree
- areaTriangleNew(iVertex) = areaTriangleNew(iVertex) + kiteAreasOnVertexNew(i,iVertex)
- enddo
- enddo
+areaTriangleNew = 0
+do iVertex=1,nVerticesNew
+do i=1,vertexDegree
+areaTriangleNew(iVertex) = areaTriangleNew(iVertex) + kiteAreasOnVertexNew(i,iVertex)
+enddo
+enddo
- allocate(edgesOnVertexNew(vertexDegreeNew, nVerticesNew))
- edgesOnVertexNew = 0
- do iVertex=1,nVertices
- if(vertexMap(iVertex).eq.0) cycle
- iVertexNew = vertexMap(iVertex)
- do i=1,vertexDegree
- j=edgesOnVertex(i,iVertex)
- jNew=edgeMap(j)
- edgesOnVertexNew(i,iVertexNew)=jNew
- enddo
- enddo
+allocate(edgesOnVertexNew(vertexDegreeNew, nVerticesNew))
+edgesOnVertexNew = 0
+do iVertex=1,nVertices
+if(vertexMap(iVertex).eq.0) cycle
+iVertexNew = vertexMap(iVertex)
+do i=1,vertexDegree
+j=edgesOnVertex(i,iVertex)
+jNew=edgeMap(j)
+edgesOnVertexNew(i,iVertexNew)=jNew
+enddo
+enddo
- ! find normals
- normalsNew = 0.0
- do iEdge=1,nEdgesNew
- cell1 = cellsOnEdgeNew(1,iEdge)
- cell2 = cellsOnEdgeNew(2,iEdge)
- if(cell1.eq.0.or.cell2.eq.0) cycle
- c1(1) = xCellNew(cell1); c1(2) = yCellNew(cell1); c1(3) = zCellNew(cell1)
- c2(1) = xCellNew(cell2); c2(2) = yCellNew(cell2); c2(3) = zCellNew(cell2)
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+! find normals
+normalsNew = 0.0
+do iEdge=1,nEdgesNew
+cell1 = cellsOnEdgeNew(1,iEdge)
+cell2 = cellsOnEdgeNew(2,iEdge)
+if(cell1.eq.0.or.cell2.eq.0) cycle
+c1(1) = xCellNew(cell1); c1(2) = yCellNew(cell1); c1(3) = zCellNew(cell1)
+c2(1) = xCellNew(cell2); c2(2) = yCellNew(cell2); c2(3) = zCellNew(cell2)
+distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- if(on_a_sphere.eq.'YES ') then
- normalsNew(1,iEdge) = c2(1) - c1(1)
- normalsNew(2,iEdge) = c2(2) - c1(2)
- normalsNew(3,iEdge) = c2(3) - c1(3)
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
- else
- if(distance.gt.0.5*Lx) then
- write(6,*) ' periodic edge ', iEdge, distance
- write(6,10) ' c1 ', c1(:)
- write(6,10) ' c2 ', c2(:)
- r = c2(1) - c1(1)
- if(r.gt.0) c2(1) = c2(1) - Lx
- if(r.lt.0) c2(1) = c2(1) + Lx
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- write(6,*) ' periodic edge fix ', iEdge, r, distance
- endif
- normalsNew(1,iEdge) = c2(1) - c1(1)
- normalsNew(2,iEdge) = c2(2) - c1(2)
- normalsNew(3,iEdge) = c2(3) - c1(3)
- distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
- normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
- endif
- enddo
- 10 format(a20,3e15.5)
+if(on_a_sphere.eq.'YES ') then
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+else
+ if(distance.gt.0.5*Lx) then
+ write(6,*) ' periodic edge ', iEdge, distance
+ write(6,10) ' c1 ', c1(:)
+ write(6,10) ' c2 ', c2(:)
+ r = c2(1) - c1(1)
+ if(r.gt.0) c2(1) = c2(1) - Lx
+ if(r.lt.0) c2(1) = c2(1) + Lx
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ write(6,*) ' periodic edge fix ', iEdge, r, distance
+ endif
+ normalsNew(1,iEdge) = c2(1) - c1(1)
+ normalsNew(2,iEdge) = c2(2) - c1(2)
+ normalsNew(3,iEdge) = c2(3) - c1(3)
+ distance = sqrt( (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + (c1(3)-c2(3))**2 )
+ normalsNew(:,iEdge) = normalsNew(:,iEdge) / distance
+endif
+enddo
+10 format(a20,3e15.5)
- end subroutine map_connectivity
+end subroutine map_connectivity
- subroutine get_dz
+subroutine get_dz
- dz( 1) = 1200.000000000000
- dz( 2) = 1595.000000000000
- dz( 3) = 2095.000000000000
- dz( 4) = 2721.000000000000
- dz( 5) = 3493.000000000000
- dz( 6) = 4431.000000000000
- dz( 7) = 5558.000000000000
- dz( 8) = 6889.000000000000
- dz( 9) = 8442.000000000000
- dz(10) = 10225.00000000000
- dz(11) = 12241.00000000000
- dz(12) = 14486.00000000000
- dz(13) = 16944.00000000000
- dz(14) = 19591.00000000000
- dz(15) = 22390.00000000000
- dz(16) = 25294.00000000000
- dz(17) = 28244.00000000000
- dz(18) = 31174.00000000000
- dz(19) = 34011.00000000000
- dz(20) = 36677.00000000000
- dz(21) = 39096.00000000000
- dz(22) = 41193.00000000000
- dz(23) = 42902.00000000000
- dz(24) = 44166.00000000000
- dz(25) = 44942.00000000000
+dz( 1) = 1200.000000000000
+dz( 2) = 1595.000000000000
+dz( 3) = 2095.000000000000
+dz( 4) = 2721.000000000000
+dz( 5) = 3493.000000000000
+dz( 6) = 4431.000000000000
+dz( 7) = 5558.000000000000
+dz( 8) = 6889.000000000000
+dz( 9) = 8442.000000000000
+dz(10) = 10225.00000000000
+dz(11) = 12241.00000000000
+dz(12) = 14486.00000000000
+dz(13) = 16944.00000000000
+dz(14) = 19591.00000000000
+dz(15) = 22390.00000000000
+dz(16) = 25294.00000000000
+dz(17) = 28244.00000000000
+dz(18) = 31174.00000000000
+dz(19) = 34011.00000000000
+dz(20) = 36677.00000000000
+dz(21) = 39096.00000000000
+dz(22) = 41193.00000000000
+dz(23) = 42902.00000000000
+dz(24) = 44166.00000000000
+dz(25) = 44942.00000000000
- dz = dz / 100.0
+dz = dz / 100.0
- end subroutine get_dz
- end program map_to_basin
+end subroutine get_dz
+end program map_to_basin
</font>
</pre>