<p><b>ringler@lanl.gov</b> 2011-10-14 11:39:11 -0600 (Fri, 14 Oct 2011)</p><p><br>
module_cullLoops is working in beta mode.<br>
it successfully removed all inland seas in the 10242 mesh<br>
the default in basin.F is true, i.e. inland seas are removed<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin/src/basin.F
===================================================================
--- branches/ocean_projects/basin/src/basin.F        2011-10-14 15:44:37 UTC (rev 1076)
+++ branches/ocean_projects/basin/src/basin.F        2011-10-14 17:39:11 UTC (rev 1077)
@@ -5,7 +5,7 @@
 use read_TS
 use write_netcdf
 use utilities
-!use cullLoops
+use cullLoops
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Program: basin.F
@@ -83,6 +83,7 @@
 !real*8, parameter :: sphere_radius = 0.0
 
 logical, parameter :: real_bathymetry=.true.
+logical, parameter :: eliminate_inland_seas=.true.
 logical, parameter :: l_woce = .true.
 
 
@@ -1193,10 +1194,12 @@
     if(flag) kmt(iCell)=0
 enddo
 
-!call eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &amp;
-                    !nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &amp;
-                    !xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &amp;
-                    !KMT)
+if(eliminate_inland_seas) then
+call eliminateLoops(nCells,nEdges,nVertices,maxEdges,vertexDegree, &amp;
+                    nEdgesOnCell, cellsOnCell, verticesOnEdge, cellsOnVertex, edgesOnCell, lonCell, latCell, &amp;
+                    xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex, &amp;
+                    KMT)
+endif
 
 if(real_bathymetry) then
     where(kmt.eq.1) kmt=3

Modified: branches/ocean_projects/basin/src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/basin/src/module_cullLoops.F        2011-10-14 15:44:37 UTC (rev 1076)
+++ branches/ocean_projects/basin/src/module_cullLoops.F        2011-10-14 17:39:11 UTC (rev 1077)
@@ -25,14 +25,17 @@
        integer, intent(inout) :: KMT(ncells)
 
        ! local workspace
-       integer :: iCell, jCell, oCell, lCell, iEdge, i, kCell, iSharedEdge, iStartEdge, iSave
+       integer :: iCell, jCell, oCell, lCell, iEdge, i, kCell, iSharedEdge, iStartEdge, iSave, iSweep
        integer :: iEdgeCounter, nEdgesInLoop(nCells), iCellAhead, LeftTurns, RightTurns
-       logical :: connected, atBoundary
-       real :: lat, rlat, rCenter(3), s(3), t(3), q(3), rCross
+       logical :: connected, atBoundary, moveSouth, moveEast, atGrenwich
+       real :: lat, rlat, rlon, rCenter(3), s(3), t(3), q(3), rCross, mylon, mylat, pi
 
+       pi = 4.0*atan(1.0)
+
        ! 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
+       moveSouth = .true.
        do iCell=1,nCells
 
          write(6,*) 'working on : ',iCell, KMT(iCell)
@@ -66,7 +69,9 @@
 
             ! 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
+            if(moveSouth) then
             rlat = 10000.0
+            mylat = latCell(jCell)
             do i=1,nEdgesOnCell(jCell)
                kCell = cellsOnCell(i,jCell)
                if(latCell(kCell).lt.rlat) then
@@ -75,8 +80,9 @@
                endif
             enddo
             jCell = iSave
+            endif
 
-            if(.not.atBoundary) write(6,*) '      pushing on to the south : ', jCell
+            if(moveSouth.and..not.atBoundary) write(6,*) '      pushing on to the south : ', jCell
               
          enddo ! .not.atBoundary
 
@@ -151,10 +157,9 @@
 
          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
+         if(LeftTurns-RightTurns.gt.0.and.rCross.gt.0.0) then
              iCellMask(iCell) = 1
              write(50,11) iCell, lonCell(iCell), latCell(iCell)
              11 format(i5,2f12.4)
@@ -163,6 +168,18 @@
        enddo ! iCell
 
        ! cull all inland seas
+       do iSweep=1,nCells/10
+         write(6,*) iSweep
+         do iCell=1,nCells
+          if(iCellMask(iCell).eq.1) then
+            do i=1,nEdgesOnCell(iCell)
+              kCell=cellsOnCell(i,iCell)
+              if(KMT(kCell).gt.0) iCellMask(kCell)=1
+            enddo
+           endif
+          enddo
+        enddo
+
        write(6,*) ' total cells culled ', sum(iCellMask)
        where(iCellMask(:).eq.1) KMT(:)=0
 

</font>
</pre>