<p><b>gaw06e@fsu.edu</b> 2011-06-04 17:37:20 -0600 (Sat, 04 Jun 2011)</p><p>- added alter_grid_for_triangle_borders (and support functions) which is called in mpas_init_block so that grid is altered just once per block. alter_grid_for_triangle_borders will modify areaCell and dvEdge as appropriately in support of 'triangle borders'<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F
===================================================================
--- branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F 2011-06-02 22:51:55 UTC (rev 873)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F 2011-06-04 23:37:20 UTC (rev 874)
@@ -53,7 +53,8 @@
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
-
+ call alter_grid_for_triangle_borders(block % mesh)
+
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call rbfInterp_initialize(mesh)
@@ -214,4 +215,188 @@
end subroutine mpas_core_finalize
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Alters grid for use with triangle border shallow water core. Specifically,
+ ! boundaryCell = 1, areaCell is appropriately reduced. In addition, where
+ ! boundary Edge = 1, dvEdge is appropriately reduced.
+ !
+ ! Input: grid - contains grid metadata
+ !
+ ! Output: grid, upon returning will have altered values
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine alter_grid_for_triangle_borders(grid)
+
+ use grid_types
+ use constants
+
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+
+ integer :: nCells, nEdges, vertexDegree
+ integer :: iEdge, vOEa, vOEb, bVa, bVb
+ integer :: iCell,iVerticesOnCell, iCellsOnVertex
+ integer :: currentVertex
+ real (kind=RKIND) :: dvEdgeTemp
+ real (kind=RKIND), dimension(3) :: pA, pB
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:), pointer :: xEdge, yEdge, zEdge
+ real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
+ real (kind=RKIND), dimension(:,:), pointer :: kiteAreasOnVertex
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: verticesOnCell, verticesOnEdge, cellsOnVertex
+ integer, dimension(:,:), pointer :: boundaryCell, boundaryEdge, boundaryVertex
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ vertexDegree = grid % vertexDegree
+
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ xEdge => grid % xEdge % array
+ yEdge => grid % yEdge % array
+ zEdge => grid % zEdge % array
+ xVertex => grid % xVertex % array
+ yVertex => grid % yVertex % array
+ zVertex => grid % zVertex % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ verticesOnCell => grid % verticesOnCell % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ boundaryCell => grid % boundaryCell % array
+ boundaryEdge => grid % boundaryEdge % array
+ boundaryVertex => grid % boundaryVertex % array
+
+ ! alter dvEdge values to represent the portion of the edge which lies
+ ! in the computational domain
+ dvEdgeDo: do iEdge=1,nEdges
+
+ appropriateEdgeIf: if (boundaryEdge(1,iEdge).eq.1) then
+
+ vOEa = verticesOnEdge(1,iEdge)
+ vOEb = verticesOnEdge(2,iEdge)
+ bVa = boundaryVertex(1,vOEa)
+ bVb = boundaryVertex(1,vOEb)
+
+ ! get and normalize coordinate for iEdge
+ pA(1) = xEdge(iEdge)
+ pA(2) = yEdge(iEdge)
+ pA(3) = zEdge(iEdge)
+ pA = pA / sqrt( dot_product(pA,pA) )
+
+ if (vOEa.eq.2) then
+
+ ! get and normalizee coordinate for vOEb
+ pB(1) = xVertex(vOEb)
+ pB(2) = yVertex(vOEb)
+ pB(3) = zVertex(vOEb)
+ pB = pB / sqrt( dot_product(pB,pB) )
+
+ ! find arc distance between two coordinates
+ dvEdgeTemp = arcDistanceR3( pA, pB )
+
+ ! set dvEdge to this value (multiplied by rEarth)
+ dvEdge(iEdge) = a * dvEdgeTemp
+
+ elseif (vOEb.eq.2) then
+
+ ! get and normalize coordinate for vOEa
+ pB(1) = xVertex(vOEa)
+ pB(2) = yVertex(vOEa)
+ pB(3) = zVertex(vOEa)
+ pB = pB / sqrt( dot_product(pB,pB) )
+
+ ! find arc distance between two coordinates
+ dvEdgeTemp = arcDistanceR3( pA, pB )
+
+ ! set dvEdge to this value ( multiplied by rEarth)
+ dvEdge(iEdge) = a * dvEdgeTemp
+
+ else
+ write(0,*) 'Error: alter_grid_for_triangle_borders problem ', &
+ ' with verticesOnEdge at edge ',iEdge,' with vOE values ', &
+ vOEA,' ',vOEB
+ end if
+
+ end if appropriateEdgeIf
+
+ end do dvEdgeDo
+
+ ! alter areaCell values (where boundaryCell == 1) to represent
+ ! the portion of the voronoi cell inside the computational domain
+ !
+ ! note: assumes that the ordering of kiteAreasOnVertex is
+ ! the same as cellsOnVertex
+
+ ! initialize areaCell
+ areaCell = 0
+
+ ! loop over all cells
+ areaCellLoop: do iCell=1,nCells
+
+ ! is this current cell a boundary cell?
+ appropriateCellIf: if (boundaryCell(1,iCell).eq.1) then
+
+ ! if it is, loop over the verticesOnCell for the current cell
+ verticesLoop: do iVerticesOnCell=1,nEdgesOnCell(iCell)
+
+ currentVertex = verticesOnCell(iVerticesOnCell,iCell)
+
+ ! for the current vertex, is is on the computational side?
+ boundaryVertexIf: if (boundaryVertex(1,currentVertex).eq.0) then
+
+ ! if it is, loop over it's consituent cells
+ triangleLoop: do iCellsOnVertex=1,vertexDegree
+
+ ! if the cellsOnVertex constituent matches the current cell,
+ ! add the corresponding kite value to the areaCell sum for
+ ! the current cell
+ cellMatchIf: if (cellsOnVertex(iCellsOnVertex,currentVertex).eq.iCell) then
+ areaCell(iCell) = areaCell(iCell) + kiteAreasOnVertex(iCellsOnVertex,currentVertex)
+ end if cellMatchIf
+
+ end do triangleLoop
+
+ end if boundaryVertexIf
+
+ end do verticesLoop
+
+ end if appropriateCellIf
+
+ end do areaCellLoop
+
+ end subroutine alter_grid_for_triangle_borders
+
+ ! returns arc distance between two points on a unit sphere
+ ! don't forget to multiply by rEarth
+ real(kind=RKIND) function arcDistanceR3( x, y )
+
+ implicit none
+
+ real(kind=RKIND), dimension(3), intent(in) :: x
+ real(kind=RKIND), dimension(3), intent(in) :: y
+
+ arcDistanceR3 = arc_cosine( dot_product( x, y ) )
+
+ end function arcDistanceR3
+
+ ! truncated argument arc cosine
+ ! adapted from John Burkardt, 2000
+ real(kind=RKIND) function arc_cosine( x )
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: x
+ real (kind=RKIND) c2
+
+ c2 = x
+ c2 = max ( c2, -1.0D+00 )
+ c2 = min ( c2, +1.0D+00 )
+
+ arc_cosine = acos ( c2 )
+
+ end function arc_cosine
+
end module mpas_core
</font>
</pre>