<p><b>ringler@lanl.gov</b> 2011-10-14 09:12:55 -0600 (Fri, 14 Oct 2011)</p><p><br>
adding a module to eliminate inland seas. still in testing mode.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/basin/src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/basin/src/module_cullLoops.F (rev 0)
+++ branches/ocean_projects/basin/src/module_cullLoops.F 2011-10-14 15:12:55 UTC (rev 1074)
@@ -0,0 +1,276 @@
+module cullLoops
+
+ public :: eliminateLoops
+
+ contains
+
+ subroutine eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+ nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &
+ xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &
+ KMT)
+
+ implicit none
+
+ ! intent (in)
+ integer :: nCells, nEdges, nVertices, maxEdges, vertexDegree
+ integer :: nEdgesOnCell(nCells), cellsOnCell(maxEdges,nCells), verticesOnEdge(2,nEdges)
+ integer :: cellsOnVertex(vertexDegree,nVertices), edgesOnCell(maxEdges,nCells)
+ real :: lonCell(nCells), latCell(nCells)
+ real :: xCell(nCells), yCell(nCells), zCell(nCells)
+ real :: xEdge(nEdges), yEdge(nEdges), zEdge(nEdges)
+ real :: xVertex(nVertices), yVertex(nVertices), zVertex(nVertices)
+ integer :: edgeList(nEdges), iCellMask(nCells)
+
+ ! intent(inout)
+ integer, intent(inout) :: KMT(ncells)
+
+ ! local workspace
+ integer :: iCell, jCell, oCell, lCell, iEdge, i, kCell, iSharedEdge, iStartEdge, iSave
+ integer :: iEdgeCounter, nEdgesInLoop(nCells), iCellAhead, LeftTurns, RightTurns
+ logical :: connected, atBoundary
+ real :: lat, rlat, rCenter(3), s(3), t(3), q(3), rCross
+
+ ! we loop over all cells and count the number of edges in the loop containing iCell
+ ! there is no coupling between iCell, so this can be inside an openMP directive
+ iCellMask(:) = 0
+ do iCell=1,nCells
+
+ write(6,*) 'working on : ',iCell, KMT(iCell)
+
+ ! skip over land cells
+ if(KMT(iCell).eq.0) then
+ write(6,*) ' skipping : ', iCell
+ cycle
+ endif
+
+ ! the working cell will be jCell, so set jCell=iCell to start
+ jCell = iCell
+ write(6,*) 'setting jCell: ', jCell
+
+ atBoundary=.false. ! are we at a boundary?
+ lCell = -1 ! when at a boundary, what is the index of the land cell or our right?
+ oCell = -1 ! when at a boundary, what is the index of the ocean cell to our left?
+
+ do while (.not.atBoundary)
+
+ ! check to see if any edges of jCell are along the boundary
+ do i=1,nEdgesOnCell(jCell)
+ kCell = cellsOnCell(i,jCell)
+ if(KMT(kCell).eq.0) then
+ lCell = kCell ! this is a land cell
+ oCell = jCell ! this is an ocean cell
+ atBoundary = .true.
+ write(6,*) ' found boundary : ',lCell,oCell
+ endif
+ enddo
+
+ ! choose the next cell to be the one with the most southern latitude
+ ! this jCell will only be used if atBoundary=.false., thus jCell must be an ocean cell
+ rlat = 10000.0
+ do i=1,nEdgesOnCell(jCell)
+ kCell = cellsOnCell(i,jCell)
+ if(latCell(kCell).lt.rlat) then
+ rlat = latCell(kCell)
+ iSave = kCell
+ endif
+ enddo
+ jCell = iSave
+
+ if(.not.atBoundary) write(6,*) ' pushing on to the south : ', jCell
+
+ enddo ! .not.atBoundary
+
+ ! OK, we hit a boundary ..... trace out the full loop in the CCW direction
+ write(6,*) ' OK we hit a boundary, let us trace out this loop '
+ write(6,*) ' ocean cell ', oCell, KMT(oCell)
+ write(6,*) ' land cell ', lCell, KMT(lCell)
+
+ ! start the counter at 1
+ iEdgeCounter = 1
+ edgeList(:) = 0
+
+ ! find the shared edge where we are starting and save the starting edge index
+ iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
+ iStartEdge = iSharedEdge
+ edgeList(iEdgeCounter) = iSharedEdge
+
+ connected = .false.
+ LeftTurns = 0; RightTurns = 0
+ do while (.not.connected)
+
+ call moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
+ oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+ verticesOnEdge,cellsOnVertex,iCellAhead)
+
+ ! if the cell ahead is ocean, then the boundary is shared between lCell and iCellAhead
+ ! if the cell ahead is land, then the boundary is shared between oCell and iCellAhead
+ if(KMT(iCellAhead).gt.0) then
+ oCell = iCellAhead
+ RightTurns = RightTurns + 1
+ write(6,*) ' the cell ahead is ocean, will turn right ', lCell, oCell
+ else
+ lCell = iCellAhead
+ LeftTurns = LeftTurns + 1
+ write(6,*) ' the cell ahead is land, will turn left ', lCell, oCell
+ endif
+ iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
+
+ ! check to see if we are where we started
+ if(iSharedEdge.eq.iStartEdge) then
+ connected=.true.
+ write(6,*) ' we are back to the where we started '
+ else
+ iEdgeCounter=iEdgeCounter+1
+ edgeList(iEdgeCounter) = iSharedEdge
+ endif
+
+ enddo ! .not.connected
+
+ ! OK, we now have a loop .... but did we circle an inland see or a land mass?
+ rCenter(:) = 0.0
+ do iEdge=1,iEdgeCounter
+ rCenter(1) = rCenter(1) + xEdge(edgeList(iEdge))/iEdgeCounter
+ rCenter(2) = rCenter(2) + yEdge(edgeList(iEdge))/iEdgeCounter
+ rCenter(3) = rCenter(3) + zEdge(edgeList(iEdge))/iEdgeCounter
+ enddo
+ rCenter(:) = rCenter(:) / sqrt ( rCenter(1)**2 + rCenter(2)**2 + rCenter(3)**2 )
+
+ rCross = 0.0
+ do iEdge=1,iEdgeCounter-1
+ t(1) = xEdge(edgeList(iEdge+1)) - xEdge(edgeList(iEdge))
+ t(2) = yEdge(edgeList(iEdge+1)) - yEdge(edgeList(iEdge))
+ t(3) = zEdge(edgeList(iEdge+1)) - zEdge(edgeList(iEdge))
+ s(1) = rCenter(1) - xEdge(edgeList(iEdge))
+ s(2) = rCenter(2) - yEdge(edgeList(iEdge))
+ s(3) = rCenter(3) - zEdge(edgeList(iEdge))
+ t(:) = t(:) / sqrt( t(1)**2 + t(2)**2 + t(3)**2 )
+ s(:) = s(:) / sqrt( s(1)**2 + s(2)**2 + s(3)**2 )
+ call cross_product_in_R3(t,s,q)
+ rCross = rCross + q(1)*rCenter(1) + q(2)*rCenter(2) + q(3)*rCenter(3)
+ enddo
+
+ write(6,*)
+ write(6,*) ' edges and cross ', iEdgeCounter, rCross, LeftTurns, RightTurns
+ if(rCross.gt.0.0) write(6,*) ' inland sea was found ', LeftTurns-RightTurns
+ write(6,*)
+
+ if(rCross.gt.0.0) then
+ iCellMask(iCell) = 1
+ write(50,11) iCell, lonCell(iCell), latCell(iCell)
+ 11 format(i5,2f12.4)
+ endif
+
+ enddo ! iCell
+
+ ! cull all inland seas
+ write(6,*) ' total cells culled ', sum(iCellMask)
+ where(iCellMask(:).eq.1) KMT(:)=0
+
+ end subroutine eliminateLoops
+
+
+ subroutine moveAhead(xCell,yCell,zCell,xVertex,yVertex,zVertex, &
+ oCell,lCell,iSharedEdge,nCells,nEdges,nVertices,maxEdges,vertexDegree, &
+ verticesOnEdge,cellsOnVertex,iCellAhead)
+ implicit none
+ integer, intent(in) :: oCell, lCell,iSharedEdge, nCells, nEdges, nVertices, maxEdges, vertexDegree
+ integer, intent(in) :: verticesOnEdge(2,nEdges), cellsOnVertex(vertexDegree,nVertices)
+ real, intent(in), dimension(nCells) :: xCell, yCell, zCell
+ real, intent(in), dimension(nVertices) :: xVertex, yVertex, zVertex
+ integer, intent(out) :: iCellAhead
+ integer :: iVertex1,iVertex2, i, kCell
+ real :: v1(3), v2(3), ocean(3), land(3), d1, d2, cross1(3), cross2(3)
+
+ ! solution assumes a CCW ordering of cellsOnVertex
+ ! the cell ahead will be connected to the vertex that lists cellsOnVertex with lCell following oCell
+
+ ! the vertex moving in the CCW direction has to be one of the two vertices connected to iSharedEdge
+ iVertex1 = verticesOnEdge(1,iSharedEdge)
+ iVertex2 = verticesOnEdge(2,iSharedEdge)
+ !write(6,*) ' iVertex1, iVertex2 ', iVertex1, iVertex2
+ !write(6,*) cellsOnVertex(:,iVertex1)
+ !write(6,*) cellsOnVertex(:,iVertex2)
+
+ v1(1)=xVertex(iVertex1)
+ v1(2)=yVertex(iVertex1)
+ v1(3)=zVertex(iVertex1)
+
+ v2(1)=xVertex(iVertex2)
+ v2(2)=yVertex(iVertex2)
+ v2(3)=zVertex(iVertex2)
+
+ ocean(1) = xCell(oCell)
+ ocean(2) = yCell(oCell)
+ ocean(3) = zCell(oCell)
+
+ land(1) = xCell(lCell)
+ land(2) = yCell(lCell)
+ land(3) = zCell(lCell)
+
+ ocean = land - ocean
+ v1 = v1 - ocean
+ v2 = v2 - ocean
+
+ ocean(:) = ocean(:) / sqrt( ocean(1)**2 + ocean(2)**2 + ocean(3)**2)
+ v1(:) = v1(:) / sqrt( v1(1)**2 + v1(2)**2 + v1(3)**2)
+ v2(:) = v2(:) / sqrt( v2(1)**2 + v2(2)**2 + v2(3)**2)
+
+ call cross_product_in_R3(ocean,v1,cross1)
+ call cross_product_in_R3(ocean,v2,cross2)
+
+ d1 = (ocean(1)+land(1))*cross1(1) + (ocean(2)+land(2))*cross1(2) + (ocean(3)+land(3))*cross1(3)
+ d2 = (ocean(1)+land(1))*cross2(1) + (ocean(2)+land(2))*cross2(2) + (ocean(3)+land(3))*cross2(3)
+
+ if(d1.gt.0.0) then
+ do i=1,vertexDegree
+ kCell=cellsOnVertex(i,iVertex1)
+ if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
+ enddo
+ endif
+
+ if(d2.gt.0.0) then
+ do i=1,vertexDegree
+ kCell=cellsOnVertex(i,iVertex2)
+ if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
+ enddo
+ endif
+
+ end subroutine moveAhead
+
+
+
+ function sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
+ implicit none
+ integer, intent(in) :: oCell, lCell, nCells, maxEdges, nEdgesOnCell(nCells), edgesOnCell(maxEdges,nCells)
+ integer :: sharedEdge
+ integer :: i,j,iEdge,jEdge
+
+ sharedEdge=-1
+ do i=1,nEdgesOnCell(oCell)
+ iEdge = edgesOnCell(i,oCell)
+ do j=1,nEdgesOnCell(lCell)
+ jEdge = edgesOnCell(j,lCell)
+ if(iEdge.eq.jEdge) then
+ sharedEdge = jEdge
+ exit
+ endif
+ enddo
+ enddo
+
+ if(SharedEdge.eq.-1) then
+ write(6,*) ' problem with SharedEdge ',oCell,lCell
+ stop
+ endif
+
+ end function sharedEdge
+
+ subroutine cross_product_in_R3(p_1,p_2,p_out)
+ real , intent(in) :: p_1 (3), p_2 (3)
+ real , intent(out) :: p_out (3)
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+ end subroutine cross_product_in_R3
+
+
+end module cullLoops
</font>
</pre>