<p><b>ringler@lanl.gov</b> 2012-03-02 12:44:19 -0700 (Fri, 02 Mar 2012)</p><p><br>
improved search algorithm of moveAhead during loop finding/culling.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/basin/src/basin.F
===================================================================
--- branches/ocean_projects/basin/src/basin.F 2012-03-02 18:30:06 UTC (rev 1579)
+++ branches/ocean_projects/basin/src/basin.F 2012-03-02 19:44:19 UTC (rev 1580)
@@ -105,7 +105,7 @@
integer :: nCellsNew, nEdgesNew, nVerticesNew
integer :: maxEdgesNew, maxEdges2New, TWONew, vertexDegreeNew, nVertLevelsNew
integer, allocatable, dimension(:) :: indexToCellIDNew, indexToEdgeIDNew, indexToVertexIDNew
-real, allocatable, dimension(:) :: xCellNew, yCellNew, zCellNew, latCellNew, lonCellNew, meshDensityNew
+real, allocatable, dimension(:) :: xCellNew, yCellNew, zCellNew, latCellNew, lonCellNew, meshDensityNew, meshSpacingNew
real, allocatable, dimension(:) :: xEdgeNew, yEdgeNew, zEdgeNew, latEdgeNew, lonEdgeNew
real, allocatable, dimension(:) :: xVertexNew, yVertexNew, zVertexNew, latVertexNew, lonVertexNew
integer, allocatable, dimension(:) :: nEdgesOnCellNew, nEdgesOnEdgeNew, flipVerticesOnEdgeOrdering
@@ -237,6 +237,7 @@
edgesOnCellNew, &
areaCellNew, &
maxLevelCellNew, &
+ meshDensityNew, &
depthCell, &
temperatureNew(1,1,:), &
kiteAreasOnVertexNew )
@@ -1313,6 +1314,7 @@
allocate(latCellNew(nCellsNew))
allocate(lonCellNew(nCellsNew))
allocate(meshDensityNew(nCellsNew))
+allocate(meshSpacingNew(nCellsNew))
allocate(xEdgeNew(nEdgesNew))
allocate(yEdgeNew(nEdgesNew))
allocate(zEdgeNew(nEdgesNew))
@@ -1348,7 +1350,7 @@
allocate(temperatureRestoreNew(nCellsNew))
allocate(salinityRestoreNew(nCellsNew))
-xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0; meshDensityNew=1.0
+xCellNew=0; yCellNew=0; zCellNew=0; latCellNew=0; lonCellNew=0; meshDensityNew=1.0; meshSpacingNew=0.0
xEdgeNew=0; yEdgeNew=0; zEdgeNew=0; latEdgeNew=0; lonEdgeNew=0
xVertexNew=0; yVertexNew=0; zVertexNew=0; latVertexNew=0; lonVertexNew=0
@@ -1527,6 +1529,7 @@
allocate(edgesOnCellNew(maxEdgesNew,nCellsNew))
edgesOnCellNew(:,:) = 0
+meshSpacingNew(:) = 0.0
do iCell=1,nCells
if(cellMap(iCell).eq.0) cycle
iCellNew = cellMap(iCell)
@@ -1535,6 +1538,7 @@
jNew = edgeMap(j)
if(jNew.eq.0) stop "edgesOnCell"
edgesOnCellNew(i,iCellNew) = jNew
+meshSpacingNew(iCellNew) = meshSpacingNew(iCellNew) + dcEdgeNew(jNew)/nEdgesOnCellNew(iCellNew)
enddo
enddo
deallocate(edgesOnCell)
Modified: branches/ocean_projects/basin/src/module_cullLoops.F
===================================================================
--- branches/ocean_projects/basin/src/module_cullLoops.F 2012-03-02 18:30:06 UTC (rev 1579)
+++ branches/ocean_projects/basin/src/module_cullLoops.F 2012-03-02 19:44:19 UTC (rev 1580)
@@ -113,11 +113,11 @@
if(KMT(iCellAhead).gt.0) then
oCell = iCellAhead
RightTurns = RightTurns + 1
- ! write(6,*) ' the cell ahead is ocean, will turn right ', lCell, oCell
+ ! 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
+ ! write(6,*) ' the cell ahead is land, will turn left ', lCell, oCell
endif
iSharedEdge = sharedEdge(oCell,lCell,nCells,maxEdges,nEdgesOnCell,edgesOnCell)
@@ -196,7 +196,7 @@
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)
+ real :: v1(3), v2(3), v3(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
@@ -224,28 +224,38 @@
land(2) = yCell(lCell)
land(3) = zCell(lCell)
- ocean = land - ocean
v1 = v1 - ocean
v2 = v2 - ocean
+ v3 = land - ocean
ocean(:) = ocean(:) / sqrt( ocean(1)**2 + ocean(2)**2 + ocean(3)**2)
+ land(:) = land(:) / sqrt( land(1)**2 + land(2)**2 + land(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)
+ v3(:) = v3(:) / sqrt( v3(1)**2 + v3(2)**2 + v3(3)**2)
- call cross_product_in_R3(ocean,v1,cross1)
- call cross_product_in_R3(ocean,v2,cross2)
+ call cross_product_in_R3(v3,v1,cross1)
+ call cross_product_in_R3(v3,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)
+ d1 = ocean(1)*cross1(1) + ocean(2)*cross1(2) + ocean(3)*cross1(3)
+ d2 = ocean(1)*cross2(1) + ocean(2)*cross2(2) + ocean(3)*cross2(3)
- if(d1.gt.0.0) then
+ kCell = 0
+
+ ! write(6,*) lCell, oCell
+ ! write(6,*) cellsOnVertex(:,iVertex1), d1
+ ! write(6,*) cellsOnVertex(:,iVertex2), d2
+ ! write(6,*) v1(:)
+ ! write(6,*) v2(:)
+
+ if(d1.gt.d2) 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
+ if(d2.gt.d1) then
do i=1,vertexDegree
kCell=cellsOnVertex(i,iVertex2)
if(kCell.ne.oCell.and.kCell.ne.lCell) iCellAhead=kCell
Modified: branches/ocean_projects/basin/src/utilities.F
===================================================================
--- branches/ocean_projects/basin/src/utilities.F 2012-03-02 18:30:06 UTC (rev 1579)
+++ branches/ocean_projects/basin/src/utilities.F 2012-03-02 19:44:19 UTC (rev 1580)
@@ -24,6 +24,7 @@
edgesOnCell, &
areaCell, &
maxLevelCell, &
+ meshSpacing, &
depthCell, &
SST, &
kiteAreasOnVertex )
@@ -48,13 +49,13 @@
integer, dimension(vertexDegree, nVertices), intent(in) :: cellsOnVertex
integer, dimension(nCells), intent(in) :: maxLevelCell
real (kind=8), dimension(nCells), intent(in) :: areaCell
- real (kind=8), dimension(nCells), intent(in) :: depthCell, SST
+ real (kind=8), dimension(nCells), intent(in) :: depthCell, SST, meshSpacing
real (kind=8), dimension(vertexDegree,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, jp1
integer :: nKitesTotal, iCell, iEdge1, iEdge2, iVertex11, iVertex12, iVertex21, iVertex22, ksave
- real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xscale, work(nCells), work1(nCells)
+ real (kind=8) :: x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xscale, work(nCells), work1(nCells), work2(nCells)
real (kind=8) :: xv, yv, zv, xc, yc, zc, dist
logical (kind=8) :: eflag
@@ -152,7 +153,7 @@
write(1,*)
a = trim('object 0 class array type float rank 0 items')
- b = trim('data file ocean.area.data')
+ b = trim('data file ocean.meshSpacing.data')
c = trim('attribute "dep" string "faces"')
write(1,10) a, nCells
write(1,10) b
@@ -174,9 +175,12 @@
close(1)
+
+ work2 = meshSpacing
work1 = depthCell
work = SST
+ open(unit= 8,file='dx/ocean.meshSpacing.data',form='formatted',status='unknown')
open(unit= 9,file='dx/ocean.depth.data',form='formatted',status='unknown')
open(unit=10,file='dx/ocean.area.data',form='formatted',status='unknown')
open(unit=11,file='dx/ocean.face.data',form='formatted',status='unknown')
@@ -187,6 +191,7 @@
iLoop = 0
iEdge = 0
do i=1,nCells
+ write(8,20) work2(i)
write(9,20) work1(i)
write(10,20) work(i)
write(11,21) i-1
</font>
</pre>