<p><b>mpetersen@lanl.gov</b> 2013-04-05 11:14:48 -0600 (Fri, 05 Apr 2013)</p><p>branch commit, tensor_operations: add computation of edgeTangentVector<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml
===================================================================
--- branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml 2013-04-04 18:01:21 UTC (rev 2729)
+++ branches/ocean_projects/tensor_operations/src/core_ocean/Registry.xml 2013-04-05 17:14:48 UTC (rev 2730)
@@ -25,6 +25,9 @@
<dim name="R3" definition="3" units="unitless"
description="The number three as a dimension."
/>
+ <dim name="SIX" definition="6" units="unitless"
+ description="The number six as a dimension."
+ />
<dim name="FIFTEEN" definition="15" units="unitless"
description="The number 15 as a dimension."
/>
@@ -717,6 +720,51 @@
<var name="rhoDisplaced" type="real" dimensions="nVertLevels nCells Time" units="kg m^{-3}"
description="potential density displaced to the mid-depth of top layer"
/>
+ <var name="strainRate2DCell" type="real" dimensions="R3 nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="2D strain rate at cell center"
+ />
+ <var name="strainRate2DVertex" type="real" dimensions="R3 nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="2D strain rate at vertices"
+ />
+ <var name="strainRate3DCell" type="real" dimensions="SIX nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="3D strain rate at cell center"
+ />
+ <var name="strainRate3DVertex" type="real" dimensions="SIX nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="3D strain rate at vertices"
+ />
+ <var name="divTensor" type="real" dimensions="nVertLevels nEdges Time" streams="o" units="m^{-1} s^{-1}"
+ description="divergence of a tensor"
+ />
+ <var name="strainRate2DCellSolution" type="real" dimensions="R3 nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="2D strain rate at cell center Solution"
+ />
+ <var name="strainRate2DVertexSolution" type="real" dimensions="R3 nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="2D strain rate at vertices"
+ />
+ <var name="strainRate3DCellSolution" type="real" dimensions="SIX nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="3D strain rate at cell center Solution"
+ />
+ <var name="strainRate3DVertexSolution" type="real" dimensions="SIX nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="3D strain rate at vertices Solution"
+ />
+ <var name="divTensorSolution" type="real" dimensions="nVertLevels nEdges Time" streams="o" units="m^{-1} s^{-1}"
+ description="divergence of a tensor"
+ />
+ <var name="strainRate2DCellDiff" type="real" dimensions="R3 nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="2D strain rate at cell center"
+ />
+ <var name="strainRate2DVertexDiff" type="real" dimensions="R3 nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="2D strain rate at vertices"
+ />
+ <var name="strainRate3DCellDiff" type="real" dimensions="SIX nVertLevels nCells Time" streams="o" units="s^{-1}"
+ description="3D strain rate at cell center error"
+ />
+ <var name="strainRate3DVertexDiff" type="real" dimensions="SIX nVertLevels nVertices Time" streams="o" units="s^{-1}"
+ description="3D strain rate at vertices error"
+ />
+ <var name="divTensorDiff" type="real" dimensions="nVertLevels nEdges Time" streams="o" units="m^{-1} s^{-1}"
+ description="divergence of a tensor"
+ />
<var name="BruntVaisalaFreqTop" type="real" dimensions="nVertLevels nCells Time" streams="o" units="s^{-2}"
description="Brunt Vaisala frequency defined at the center (horizontally) and top (vertically) of cell"
/>
@@ -887,6 +935,9 @@
<var name="edgeNormalVectors" type="real" dimensions="R3 nEdges" streams="o" units="unitless"
description="Normal vector defined at an edge."
/>
+ <var name="edgeTangentVectors" type="real" dimensions="R3 nEdges" streams="o" units="unitless"
+ description="Tangent vector defined at an edge."
+ />
<var name="localVerticalUnitVectors" type="real" dimensions="R3 nCells" streams="o" units="unitless"
description="Unit surface normal vectors defined at cell centers."
/>
Modified: branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F 2013-04-04 18:01:21 UTC (rev 2729)
+++ branches/ocean_projects/tensor_operations/src/core_ocean/mpas_ocn_mpas_core.F 2013-04-05 17:14:48 UTC (rev 2730)
@@ -5,6 +5,8 @@
use mpas_timekeeping
use mpas_dmpar
use mpas_timer
+ ! mrp added to test stress tensor:
+! use mpas_tensor
use ocn_global_diagnostics
use ocn_time_integration
@@ -121,6 +123,13 @@
call ocn_compute_max_level(domain)
+ ! mrp added just to test stress tensor ops:
+! print *, 'mpas_test_stress_tensor'
+! call mpas_test_tensor(domain)
+! print *, 'mpas_test_stress_tensor complete'
+! stop
+ ! mrp added just to test stress tensor ops end
+
if (.not.config_do_restart) call ocn_init_split_timestep(domain)
write (0,'(a,a)') ' Vertical coordinate movement is: ',trim(config_vert_coord_movement)
Modified: branches/ocean_projects/tensor_operations/src/operators/Makefile
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/Makefile 2013-04-04 18:01:21 UTC (rev 2729)
+++ branches/ocean_projects/tensor_operations/src/operators/Makefile 2013-04-05 17:14:48 UTC (rev 2730)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = mpas_rbf_interpolation.o mpas_vector_reconstruction.o mpas_spline_interpolation.o
+OBJS = mpas_rbf_interpolation.o mpas_vector_reconstruction.o mpas_spline_interpolation.o mpas_tensor.o
all: operators
@@ -10,6 +10,7 @@
mpas_vector_reconstruction.o: mpas_rbf_interpolation.o
mpas_rbf_interpolation.o:
mpas_spline_interpolation:
+mpas_tensor.o:
clean:
$(RM) *.o *.mod *.f90 libops.a
Modified: branches/ocean_projects/tensor_operations/src/operators/mpas_rbf_interpolation.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/mpas_rbf_interpolation.F 2013-04-04 18:01:21 UTC (rev 2729)
+++ branches/ocean_projects/tensor_operations/src/operators/mpas_rbf_interpolation.F 2013-04-05 17:14:48 UTC (rev 2730)
@@ -129,10 +129,10 @@
type (mesh_type), intent(inout) :: grid !< Input/Output: Grid information
integer :: nCells, nEdges
- integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
- integer :: iEdge, iCell, cell1, cell2
- real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
- real(kind=RKIND), dimension(:,:), pointer :: localVerticalUnitVectors, edgeNormalVectors
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgesOnCell
+ integer :: iEdge, iCell, cell1, cell2, vertex1, vertex2
+ real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge, xVertex, yVertex, zVertex
+ real(kind=RKIND), dimension(:,:), pointer :: localVerticalUnitVectors, edgeNormalVectors, edgeTangentVectors
real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
real(kind=RKIND), dimension(3) :: xHatPlane, yHatPlane, rHat
real(kind=RKIND) :: normalDotRHat
@@ -144,7 +144,11 @@
xEdge => grid % xEdge % array
yEdge => grid % yEdge % array
zEdge => grid % zEdge % array
+ xVertex => grid % xVertex % array
+ yVertex => grid % yVertex % array
+ zVertex => grid % zVertex % array
cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
edgesOnCell => grid % edgesOnCell % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -152,10 +156,12 @@
localVerticalUnitVectors => grid % localVerticalUnitVectors % array
edgeNormalVectors => grid % edgeNormalVectors % array
+ edgeTangentVectors => grid % edgeTangentVectors % array
cellTangentPlane => grid % cellTangentPlane % array
! init arrays
edgeNormalVectors = 0
+ edgeTangentVectors = 0
localVerticalUnitVectors = 0
! loop over all cells to be solved on this block
@@ -170,20 +176,60 @@
end if
end do
+! do iEdge = 1,nEdges
+! iCell = cellsOnEdge(1,iEdge) ! the normal vector points from the first cell toward the edge
+! if(iCell == nCells+1) then ! this is a boundary edge
+! ! the first cell bordering this edge is not real, use the second cell
+! ! The normal should always point outward at boundaries, away from the valid cell center
+! iCell = cellsOnEdge(2,iEdge)
+! end if
+! ! the normal points from the cell location to the edge location
+! edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(iCell)
+! edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(iCell)
+! edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(iCell)
+! call mpas_unit_vec_in_r3(edgeNormalVectors(:,iEdge))
+! end do
+
+ ! mrp 130329, replace above section. For 3D, I need vector to point from cell to cell,
+ ! not cell to edge. This must be modified for boundaries.
do iEdge = 1,nEdges
- iCell = cellsOnEdge(1,iEdge) ! the normal vector points from the first cell toward the edge
- if(iCell == nCells+1) then ! this is a boundary edge
- ! the first cell bordering this edge is not real, use the second cell
- ! The normal should always point outward at boundaries, away from the valid cell center
- iCell = cellsOnEdge(2,iEdge)
- end if
- ! the normal points from the cell location to the edge location
- edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(iCell)
- edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(iCell)
- edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(iCell)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ if (cell1 == nCells+1) then ! this is a boundary edge
+ ! the normal points from the cell location to the edge location
+ edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(cell2)
+ edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(cell2)
+ edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(cell2)
+
+ elseif (cell2 == nCells+1) then ! this is a boundary edge
+ ! the normal points from the cell location to the edge location
+ edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(cell1)
+ edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(cell1)
+ edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(cell1)
+
+ else ! this is not a boundary cell
+ ! the normal points from the cell 1 to cell2
+ edgeNormalVectors(1,iEdge) = xCell(cell2) - xCell(cell1)
+ edgeNormalVectors(2,iEdge) = yCell(cell2) - yCell(cell1)
+ edgeNormalVectors(3,iEdge) = zCell(cell2) - zCell(cell1)
+
+ endif
call mpas_unit_vec_in_r3(edgeNormalVectors(:,iEdge))
end do
+ ! mrp 130329, added tangent vectors. This is not required for rbf, so could be moved to a
+ ! tensor init subroutine later.
+ do iEdge = 1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ ! the tangent vector points from the vertex 1 to vertex 2
+ edgeTangentVectors(1,iEdge) = xVertex(vertex2) - xVertex(vertex1)
+ edgeTangentVectors(2,iEdge) = yVertex(vertex2) - yVertex(vertex1)
+ edgeTangentVectors(3,iEdge) = zVertex(vertex2) - zVertex(vertex1)
+ call mpas_unit_vec_in_r3(edgeTangentVectors(:,iEdge))
+ end do
+
do iCell=1,nCells
iEdge = edgesOnCell(1,iCell)
! xHat and yHat are a local basis in the plane of the horizontal cell
Added: branches/ocean_projects/tensor_operations/src/operators/mpas_tensor.F
===================================================================
--- branches/ocean_projects/tensor_operations/src/operators/mpas_tensor.F (rev 0)
+++ branches/ocean_projects/tensor_operations/src/operators/mpas_tensor.F 2013-04-05 17:14:48 UTC (rev 2730)
@@ -0,0 +1,748 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! mpas_tensor
+!
+!> \brief MPAS tensor operations
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> the strain rate tensor, the divergence of a tensor,
+!> and a testing routine to verify these work properly.
+!
+!-----------------------------------------------------------------------
+
+module mpas_tensor
+
+ use mpas_grid_types
+ use mpas_constants
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: mpas_strain_rate, &
+ mpas_div_tensor, &
+ mpas_test_tensor
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine mpas_strain_rate
+!
+!> \brief Computes strain rate at cells and vertices
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the strain rate at the cell centers and
+!> vertices using the weak derivative.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_strain_rate(grid, u,v, strainRate2DCell, strainRate2DVertex)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u, &!< Input: Horizontal velocity normal to edge
+ v !< Input: Horizontal velocity tangent to edge
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(out) :: &
+ strainRate2DCell !< Output: strain rate tensor at cell center
+ real (kind=RKIND), dimension(:,:,:), intent(out) :: &
+ strainRate2DVertex !< Output: strain rate tensor at vertex
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, nCells,nVertices, cell1, cell2, vertex1, vertex2, k
+
+ integer, dimension(:), pointer :: maxLevelEdgeBot
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge
+
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, &
+ cos_pwr1, sin_pwr1, cos_pwr2, sin_pwr2, cossin, uDvEdge, uDcEdge, eps11, eps22, eps12, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, areaTriangle, angleEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ nEdges = grid % nEdges
+ nCells = grid % nCells
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ angleEdge => grid % angleEdge % array
+
+ strainRate2DCell = 0.0
+ strainRate2DVertex = 0.0
+
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaTri1 = 1.0 / areaTriangle(vertex1)
+ invAreaTri2 = 1.0 / areaTriangle(vertex2)
+
+ !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
+ invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0)
+ invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0)
+
+ cos_pwr1 = cos(angleEdge(iEdge))
+ cos_pwr2 = cos_pwr1**2
+ sin_pwr1 = sin(angleEdge(iEdge))
+ sin_pwr2 = sin_pwr1**2
+ cossin = cos_pwr1*sin_pwr1
+
+ ! change later:
+ !do k=1,maxLevelEdgeBot(iEdge)
+ do k=1,nVertLevels
+
+ ! strain rate tensor at cell center
+ eps11 = u(k,iEdge)*cos_pwr2 - v(k,iEdge)*cossin
+ eps22 = u(k,iEdge)*sin_pwr2 + v(k,iEdge)*cossin
+ eps12 = u(k,iEdge)*cossin + v(k,iEdge)*0.5*(cos_pwr2-sin_pwr2)
+
+ ! index 1, 2, and 3 is for 11, 22, and 12
+ strainRate2DCell(1,k,cell1) = strainRate2DCell(1,k,cell1) + eps11*invAreaCell1 * dvEdge(iEdge)
+ strainRate2DCell(1,k,cell2) = strainRate2DCell(1,k,cell2) - eps11*invAreaCell2 * dvEdge(iEdge)
+ strainRate2DCell(2,k,cell1) = strainRate2DCell(2,k,cell1) + eps22*invAreaCell1 * dvEdge(iEdge)
+ strainRate2DCell(2,k,cell2) = strainRate2DCell(2,k,cell2) - eps22*invAreaCell2 * dvEdge(iEdge)
+ strainRate2DCell(3,k,cell1) = strainRate2DCell(3,k,cell1) + eps12*invAreaCell1 * dvEdge(iEdge)
+ strainRate2DCell(3,k,cell2) = strainRate2DCell(3,k,cell2) - eps12*invAreaCell2 * dvEdge(iEdge)
+
+ ! strain rate tensor at vertex
+ eps11 = -v(k,iEdge)*sin_pwr2 + u(k,iEdge)*cossin
+ eps22 = -v(k,iEdge)*cos_pwr2 - u(k,iEdge)*cossin
+ eps12 = v(k,iEdge)*cossin - u(k,iEdge)*0.5*(cos_pwr2-sin_pwr2)
+
+ strainRate2DVertex(1,k,vertex1) = strainRate2DVertex(1,k,vertex1) - eps11*invAreaTri1* dcEdge(iEdge)
+ strainRate2DVertex(1,k,vertex2) = strainRate2DVertex(1,k,vertex2) + eps11*invAreaTri2* dcEdge(iEdge)
+ strainRate2DVertex(2,k,vertex1) = strainRate2DVertex(2,k,vertex1) - eps22*invAreaTri1* dcEdge(iEdge)
+ strainRate2DVertex(2,k,vertex2) = strainRate2DVertex(2,k,vertex2) + eps22*invAreaTri2* dcEdge(iEdge)
+ strainRate2DVertex(3,k,vertex1) = strainRate2DVertex(3,k,vertex1) - eps12*invAreaTri1* dcEdge(iEdge)
+ strainRate2DVertex(3,k,vertex2) = strainRate2DVertex(3,k,vertex2) + eps12*invAreaTri2* dcEdge(iEdge)
+
+ end do
+ end do
+
+! print '(a,100f8.3)', 'xCell ',grid % xCell % array(6:7)
+! print '(a,100f8.3)', 'yCell ',grid % yCell % array(6:7)
+! print '(a,100f8.3)', 'xEdge ',grid % xEdge % array(1:nEdges)
+! print '(a,100f8.3)', 'yEdge ',grid % yEdge % array(1:nEdges)
+! print '(a,100f8.3)', 'angleEdge ',angleEdge(1:nEdges)
+! print '(a,100f8.3)', 'u ',u(:,1:nEdges)
+! print '(a,100f8.3)', 'v ',v(:,1:nEdges)
+! print '(a,100f8.3)', 'eps Nor 1 ',strainRate2DCell(1,:,6:7)
+! print '(a,100f8.3)', 'eps Nor 2 ',strainRate2DCell(2,:,6:7)
+! print '(a,100f8.3)', 'eps Nor 3 ',strainRate2DCell(3,:,6:7)
+! print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(10:11)
+! print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(10:11)
+! print '(a,100f8.3)', 'eps Shr 1 ',strainRate2DVertex(1,:,10:11)
+! print '(a,100f8.3)', 'eps Shr 2 ',strainRate2DVertex(2,:,10:11)
+! print '(a,100f8.3)', 'eps Shr 3 ',strainRate2DVertex(3,:,10:11)
+
+if (1==2) then
+ print '(a,100f8.3)', 'xCell ',grid % xCell % array(1:nCells)
+ print '(a,100f8.3)', 'yCell ',grid % yCell % array(1:nCells)
+ print '(a,100f8.3)', 'xEdge ',grid % xEdge % array(1:nEdges)
+ print '(a,100f8.3)', 'yEdge ',grid % yEdge % array(1:nEdges)
+ print '(a,100f8.3)', 'angleEdge ',angleEdge(1:nEdges)
+ print '(a,100f8.3)', 'u ',u(:,1:nEdges)
+ print '(a,100f8.3)', 'v ',v(:,1:nEdges)
+ print '(a,100f8.3)', 'eps Nor 1 ',strainRate2DCell(1,:,1:nCells)
+ print '(a,100f8.3)', 'eps Nor 2 ',strainRate2DCell(2,:,1:nCells)
+ print '(a,100f8.3)', 'eps Nor 3 ',strainRate2DCell(3,:,1:nCells)
+ print '(a,100f8.3)', 'xVertex ',grid % xVertex % array(1:nVertices)
+ print '(a,100f8.3)', 'yVertex ',grid % yVertex % array(1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 1 ',strainRate2DVertex(1,:,1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 2 ',strainRate2DVertex(2,:,1:nVertices)
+ print '(a,100f8.3)', 'eps Shr 3 ',strainRate2DVertex(3,:,1:nVertices)
+endif
+
+
+ end subroutine mpas_strain_rate!}}}
+
+!***********************************************************************
+!
+! routine mpas_div_tensor
+!
+!> \brief Computes divergence of the stress tensor
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the divergence of the stress tensor
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_div_tensor(grid, tensorCell, tensorVertex, divTensor)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tensorCell !< Input: tensor at cell center
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tensorVertex !< Input: tensor at vertex
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ divTensor !< Output: normal component of divergence of the stress tensor at edge
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2, k
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, &
+ cos_pwr1, sin_pwr1, cos_pwr2, sin_pwr2, cossin, uDvEdge, uDcEdge, eps11, eps22, eps12, &
+ invdcEdge, invdvEdge
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, areaTriangle, angleEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ nEdgesSolve = grid % nEdgesSolve
+
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ angleEdge => grid % angleEdge % array
+ edgeMask => grid % edgeMask % array
+
+ do iEdge=1,nEdgesSolve
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ cos_pwr1 = cos(angleEdge(iEdge))
+ cos_pwr2 = cos_pwr1**2
+ sin_pwr1 = sin(angleEdge(iEdge))
+ sin_pwr2 = sin_pwr1**2
+ cossin = cos_pwr1*sin_pwr1
+
+ invdcEdge = 1.0 / dcEdge(iEdge)
+ invdvEdge = 1.0 / dvEdge(iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ divTensor(k,iEdge) = edgeMask(k,iEdge) * invdcEdge &
+ * ( cos_pwr2*tensorCell(1,k,cell2) - 2.0*cossin*tensorCell(3,k,cell2) + sin_pwr2*tensorCell(2,k,cell2) &
+ -(cos_pwr2*tensorCell(1,k,cell1) - 2.0*cossin*tensorCell(3,k,cell1) + sin_pwr2*tensorCell(2,k,cell1)) ) &
+ + edgeMask(k,iEdge) * invdvEdge &
+ * ( (cos_pwr2-sin_pwr2)*tensorVertex(3,k,vertex2) + cossin*(tensorVertex(1,k,vertex2) - tensorVertex(2,k,vertex2)) &
+ -((cos_pwr2-sin_pwr2)*tensorVertex(3,k,vertex1) + cossin*(tensorVertex(1,k,vertex1) - tensorVertex(2,k,vertex1))))
+ end do
+
+ end do
+
+ end subroutine mpas_div_tensor!}}}
+
+!***********************************************************************
+!
+! routine mpas_test_tensor
+!
+!> \brief Tests strain rate and tensor divergence operators
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine tests strain rate and tensor divergence operators.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_test_tensor(domain) !{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (domain_type), intent(inout) :: domain
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ type (block_type), pointer :: block
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, iCell, iEdge, iVertex, k, i, p
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnEdge
+
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, &
+ cos_pwr1, sin_pwr1, cos_pwr2, sin_pwr2, cossin, uDvEdge, uDcEdge, eps11, eps22, eps12, &
+ invdcEdge, invdvEdge, ux, uy, eoe, cn, cs, r, theta, rot, f, g1, g2, fcos, pi2l, ld, &
+ divSigma_x, divSigma_y
+ real (kind=RKIND), dimension(:), pointer :: angleEdge, &
+ xCell, yCell, xEdge, yEdge, xVertex, yVertex
+ real (kind=RKIND), dimension(:,:), pointer :: u, v, divTensor, weightsOnEdge
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRate2DCell, strainRate2DVertex
+! real (kind=RKIND), dimension(:,:,:), allocatable :: strainRate2DCellSolution, strainRate2DVertexSolution, divTensorSolution
+ real (kind=RKIND), dimension(:,:), pointer :: divTensorSolution
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRate2DCellSolution, strainRate2DVertexSolution
+ character(len=16) :: test_function_type
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+! maxLevelCell => block % mesh % maxLevelCell % array
+! maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+! maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+! maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+! maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+! cellsOnEdge => block % mesh % cellsOnEdge % array
+! cellsOnVertex => block % mesh % cellsOnVertex % array
+! verticesOnEdge => block % mesh % verticesOnEdge % array
+! boundaryEdge => block % mesh % boundaryEdge % array
+! boundaryCell => block % mesh % boundaryCell % array
+! boundaryVertex => block % mesh % boundaryVertex % array
+! edgeMask => block % mesh % edgeMask % array
+! cellMask => block % mesh % cellMask % array
+! vertexMask => block % mesh % vertexMask % array
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+
+ xCell => block % mesh % xCell % array
+ yCell => block % mesh % yCell % array
+ xEdge => block % mesh % xEdge % array
+ yEdge => block % mesh % yEdge % array
+ xVertex => block % mesh % xVertex % array
+ yVertex => block % mesh % yVertex % array
+ angleEdge => block % mesh % angleEdge % array
+
+ nEdgesOnEdge => block % mesh % nEdgesOnEdge % array
+ edgesOnEdge => block % mesh % edgesOnEdge % array
+ weightsOnEdge => block % mesh % weightsOnEdge % array
+
+ u => block % state % time_levs(1) % state % u % array
+ v => block % state % time_levs(1) % state % v % array
+ strainRate2DCell => block % state % time_levs(1) % state % strainRate2DCell % array
+ strainRate2DVertex => block % state % time_levs(1) % state % strainRate2DVertex % array
+ divTensor => block % state % time_levs(1) % state % divTensor % array
+ strainRate2DCellSolution => block % state % time_levs(1) % state % strainRate2DCellSolution % array
+ strainRate2DVertexSolution => block % state % time_levs(1) % state % strainRate2DVertexSolution % array
+ divTensorSolution => block % state % time_levs(1) % state % divTensorSolution % array
+
+! allocate(strainRate2DCellSolution(3,nVertLevels,nCells),strainRate2DVertexSolution(3,nVertLevels,nVertices),divTensorSolution(3,nVertLevels,nEdges))
+
+ cn = 1.555
+ cs = 2.555
+ rot = 3.0 !pii/4.0 !0.0*pii/2
+ p = 3 ! power for polynomial function
+! test_function_type = 'linear_x'
+! test_function_type = 'linear_y'
+! test_function_type = 'linear_arb_rot'
+! test_function_type = 'sin_arb_rot'
+! test_function_type = 'power_x'
+! test_function_type = 'power_y'
+ test_function_type = 'power_arb_rot'
+
+ ld = maxval(xEdge) ! note: this will not work in parallel
+ pi2l = 2*pii/ld
+ g1 = cn*cos(rot) - cs*sin(rot)
+ g2 = cn*sin(rot) + cs*cos(rot)
+
+ if (test_function_type.eq.'linear_x') then
+ ! linear function in x
+ do iEdge = 1,nEdges
+ ux = cn*xEdge(iEdge)
+ uy = cs*xEdge(iEdge)
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+ enddo
+
+ strainRate2DCellSolution(1,:,:) = cn
+ strainRate2DCellSolution(2,:,:) = 0.0
+ strainRate2DCellSolution(3,:,:) = 0.5*cs
+
+ strainRate2DVertexSolution(1,:,:) = cn
+ strainRate2DVertexSolution(2,:,:) = 0.0
+ strainRate2DVertexSolution(3,:,:) = 0.5*cs
+
+ divTensorSolution = 0.0
+
+ elseif (test_function_type.eq.'linear_y') then
+
+ ! linear function in y
+ do iEdge = 1,nEdges
+ ux = -cs*yEdge(iEdge)
+ uy = cn*yEdge(iEdge)
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+ enddo
+
+ strainRate2DCellSolution(1,:,:) = 0.0
+ strainRate2DCellSolution(2,:,:) = cn
+ strainRate2DCellSolution(3,:,:) = -0.5*cs
+
+ strainRate2DVertexSolution(1,:,:) = 0.0
+ strainRate2DVertexSolution(2,:,:) = cn
+ strainRate2DVertexSolution(3,:,:) = -0.5*cs
+
+ divTensorSolution = 0.0
+
+ elseif (test_function_type.eq.'linear_arb_rot') then
+
+ ! linear function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+
+ f = r*cos(theta-rot)
+ ux = f*g1
+ uy = f*g2
+
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ enddo
+
+ strainRate2DCellSolution(1,:,:) = cos(rot)*g1
+ strainRate2DCellSolution(2,:,:) = sin(rot)*g2
+ strainRate2DCellSolution(3,:,:) = 0.5*(cos(rot)*g2 + sin(rot)*g1)
+
+ strainRate2DVertexSolution(1,:,:) = cos(rot)*g1
+ strainRate2DVertexSolution(2,:,:) = sin(rot)*g2
+ strainRate2DVertexSolution(3,:,:) = 0.5*(cos(rot)*g2 + sin(rot)*g1)
+
+ divTensorSolution = 0.0
+
+ elseif (test_function_type.eq.'power_x') then
+
+ ! power function in x: x^p
+ do iEdge = 1,nEdges
+ ux = cn*xEdge(iEdge)**p
+ uy = cs*xEdge(iEdge)**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = cn *p*(p-1)*xEdge(iEdge)**(p-2)
+ divSigma_y = 0.5*cs*p*(p-1)*xEdge(iEdge)**(p-2)
+
+ divTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ strainRate2DCellSolution(1,:,iCell) = cn *p*xCell(iCell)**(p-1)
+ strainRate2DCellSolution(2,:,iCell) = 0.0
+ strainRate2DCellSolution(3,:,iCell) = 0.5*cs*p*xCell(iCell)**(p-1)
+ end do
+
+ do iVertex = 1,nVertices
+ strainRate2DVertexSolution(1,:,iVertex) = cn *p*xVertex(iVertex)**(p-1)
+ strainRate2DVertexSolution(2,:,iVertex) = 0.0
+ strainRate2DVertexSolution(3,:,iVertex) = 0.5*cs*p*xVertex(iVertex)**(p-1)
+ end do
+
+ elseif (test_function_type.eq.'power_y') then
+
+ ! power function in y: y^n
+ do iEdge = 1,nEdges
+ ux = -cs*yEdge(iEdge)**p
+ uy = cn*yEdge(iEdge)**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = -0.5*cs*p*(p-1)*yEdge(iEdge)**(p-2)
+ divSigma_y = cn *p*(p-1)*yEdge(iEdge)**(p-2)
+
+ divTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ strainRate2DCellSolution(1,:,iCell) = 0.0
+ strainRate2DCellSolution(2,:,iCell) = cn *p*yCell(iCell)**(p-1)
+ strainRate2DCellSolution(3,:,iCell) = -0.5*cs*p*yCell(iCell)**(p-1)
+ end do
+
+ do iVertex = 1,nVertices
+ strainRate2DVertexSolution(1,:,iVertex) = 0.0
+ strainRate2DVertexSolution(2,:,iVertex) = cn *p*yVertex(iVertex)**(p-1)
+ strainRate2DVertexSolution(3,:,iVertex) = -0.5*cs*p*yVertex(iVertex)**(p-1)
+ end do
+
+ elseif (test_function_type.eq.'power_arb_rot') then
+
+ ! power function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ f = r*cos(theta-rot)
+
+ ux = g1*f**p
+ uy = g2*f**p
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = p*(p-1)*f**(p-2) *(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+ divSigma_y = p*(p-1)*f**(p-2) *(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+
+ divTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ r = sqrt(xCell(iCell)**2 + yCell(iCell)**2)
+ theta = atan(yCell(iCell)/xCell(iCell))
+ f = r*cos(theta-rot)
+
+ strainRate2DCellSolution(1,:,iCell) = p *f**(p-1) *cos(rot)*g1
+ strainRate2DCellSolution(2,:,iCell) = p *f**(p-1) *sin(rot)*g2
+ strainRate2DCellSolution(3,:,iCell) = p *f**(p-1) *(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+ do iVertex = 1,nVertices
+ r = sqrt(xVertex(iVertex)**2 + yVertex(iVertex)**2)
+ theta = atan(yVertex(iVertex)/xVertex(iVertex))
+ f = r*cos(theta-rot)
+
+ strainRate2DVertexSolution(1,:,iVertex) = p *f**(p-1) *cos(rot)*g1
+ strainRate2DVertexSolution(2,:,iVertex) = p *f**(p-1) *sin(rot)*g2
+ strainRate2DVertexSolution(3,:,iVertex) = p *f**(p-1) *(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+
+ elseif (test_function_type.eq.'sin_arb_rot') then
+
+ ! sin function, arbitrary rotation
+
+ do iEdge = 1,nEdges
+ ! arbitrary rotation
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ f = sin(pi2l*r*cos(theta-rot))
+
+ ux = f*g1
+ uy = f*g2
+ u(:,iEdge) = ux*cos(angleEdge(iEdge)) + uy*sin(angleEdge(iEdge))
+! v(:,iEdge) = ux*sin(angleEdge(iEdge)) - uy*cos(angleEdge(iEdge))
+
+ ! solution for div(sigma)
+ divSigma_x = -pi2l**2*f*(cos(rot)**2*g1 + 0.5*(sin(rot)**2*g1 + sin(rot)*cos(rot)*g2) )
+ divSigma_y = -pi2l**2*f*(sin(rot)**2*g2 + 0.5*(cos(rot)**2*g2 + sin(rot)*cos(rot)*g1) )
+
+ divTensorSolution(:,iEdge) = divSigma_x*cos(angleEdge(iEdge)) + divSigma_y*sin(angleEdge(iEdge))
+ enddo
+
+ do iCell = 1,nCells
+ r = sqrt(xCell(iCell)**2 + yCell(iCell)**2)
+ theta = atan(yCell(iCell)/xCell(iCell))
+ fcos = cos(pi2l*r*cos(theta-rot))
+
+ strainRate2DCellSolution(1,:,iCell) = pi2l*fcos*cos(rot)*g1
+ strainRate2DCellSolution(2,:,iCell) = pi2l*fcos*sin(rot)*g2
+ strainRate2DCellSolution(3,:,iCell) = pi2l*fcos*(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+ do iVertex = 1,nVertices
+ r = sqrt(xVertex(iVertex)**2 + yVertex(iVertex)**2)
+ theta = atan(yVertex(iVertex)/xVertex(iVertex))
+ fcos = cos(pi2l*r*cos(theta-rot))
+
+ strainRate2DVertexSolution(1,:,iVertex) = pi2l*fcos*cos(rot)*g1
+ strainRate2DVertexSolution(2,:,iVertex) = pi2l*fcos*sin(rot)*g2
+ strainRate2DVertexSolution(3,:,iVertex) = pi2l*fcos*(cos(rot)*g2+sin(rot)*g1)/2.0
+ end do
+
+ else
+ print *, 'bad choice of test_function_type:',test_function_type
+ stop
+ endif
+
+ do iEdge = 1,nEdges
+ ! Compute v (tangential) velocities
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k=1,nVertLevels
+ !do k = 1,maxLevelEdgeTop(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
+ end do
+ end do
+
+ call mpas_strain_rate(block % mesh, u,v, strainRate2DCell, strainRate2DVertex)
+
+ call mpas_div_tensor(block % mesh, strainRate2DCell, strainRate2DVertex, divTensor)
+
+if (1==2) then
+ print *, 'nedges, size(u)',nEdges, size(u)
+ print '(a,3f20.10)', 'min max mean u ',minval(u(:,1:nEdges)),maxval(u(:,1:nEdges)),sum(u(:,1:nEdges))/nEdges/nVertLevels
+ print '(a,3f20.10)', 'min max mean v ',minval(v(:,1:nEdges)),maxval(v(:,1:nEdges)),sum(v(:,1:nEdges))/nEdges/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DCell1',minval(strainRate2DCell(1,:,1:nCells)),maxval(strainRate2DCell(1,:,1:nCells)),sum(strainRate2DCell(1,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DCell2',minval(strainRate2DCell(2,:,1:nCells)),maxval(strainRate2DCell(2,:,1:nCells)),sum(strainRate2DCell(2,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DCell3',minval(strainRate2DCell(3,:,1:nCells)),maxval(strainRate2DCell(3,:,1:nCells)),sum(strainRate2DCell(3,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DVertex1 ',minval(strainRate2DVertex(1,:,1:nVertices)),maxval(strainRate2DVertex(1,:,1:nVertices)),sum(strainRate2DVertex(1,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DVertex2 ',minval(strainRate2DVertex(2,:,1:nVertices)),maxval(strainRate2DVertex(2,:,1:nVertices)),sum(strainRate2DVertex(2,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRate2DVertex3 ',minval(strainRate2DVertex(3,:,1:nVertices)),maxval(strainRate2DVertex(3,:,1:nVertices)),sum(strainRate2DVertex(3,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean divTensor ',minval(divTensor(:,1:nEdges)),maxval(divTensor(:,1:nEdges)),sum(divTensor(:,1:nEdges))/nEdges/nVertLevels
+endif
+
+ block % state % time_levs(1) % state % strainRate2DCellDiff % array = abs(strainRate2DCell- strainRate2DCellSolution)
+ block % state % time_levs(1) % state % strainRate2DVertexDiff % array = abs(strainRate2DVertex- strainRate2DVertexSolution)
+ block % state % time_levs(1) % state % divTensorDiff % array = abs(divTensor- divTensorSolution)
+
+ print '(a,3es20.10)', 'max abs diff strainRate2DCell1',maxval(abs(strainRate2DCell(1,:,1:nCells) - strainRate2DCellSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff strainRate2DCell2',maxval(abs(strainRate2DCell(2,:,1:nCells) - strainRate2DCellSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff strainRate2DCell3',maxval(abs(strainRate2DCell(3,:,1:nCells) - strainRate2DCellSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff strainRate2DVertex1 ',maxval(abs(strainRate2DVertex(1,:,1:nVertices) - strainRate2DVertexSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff strainRate2DVertex2 ',maxval(abs(strainRate2DVertex(2,:,1:nVertices) - strainRate2DVertexSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff strainRate2DVertex3 ',maxval(abs(strainRate2DVertex(3,:,1:nVertices) - strainRate2DVertexSolution(1,:,1:nCells)))
+ print '(a,3es20.10)', 'max abs diff divTensor ',maxval(abs(divTensor(:,1:nEdges) - divTensorSolution(:,1:nEdges)))
+
+
+! deallocate(strainRate2DCellSolution,strainRate2DVertexSolution,divTensorSolution)
+
+ block => block % next
+ end do
+
+ end subroutine mpas_test_tensor!}}}
+
+!***********************************************************************
+!
+! routine mpas_tensor_init
+!
+!> \brief Initializes flags used within tendency routines.
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes flags related to quantities computed within
+!> other tendency routines.
+!
+!-----------------------------------------------------------------------
+! subroutine mpas_tensor_init(err)!{{{
+! integer, intent(out) :: err !< Output: Error flag
+
+! end subroutine mpas_tensor_init!}}}
+
+!***********************************************************************
+
+end module mpas_tensor
+
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
+
</font>
</pre>