<p><b>mpetersen@lanl.gov</b> 2012-09-27 12:32:06 -0600 (Thu, 27 Sep 2012)</p><p>Added mpas_stress_tensor to operator directory on this branch. Curently tested with linear velocity functions on a quad grid.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/strain_rate/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/strain_rate/src/core_ocean/Registry 2012-09-26 20:56:32 UTC (rev 2169)
+++ branches/ocean_projects/strain_rate/src/core_ocean/Registry 2012-09-27 18:32:06 UTC (rev 2170)
@@ -276,6 +276,9 @@
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
+var persistent real strainRateNormal ( R3 nVertLevels nCells Time ) 2 o strainRateNormal state - -
+var persistent real strainRateShear ( R3 nVertLevels nVertices Time ) 2 o strainRateShear state - -
+var persistent real divStressTensor ( nVertLevels nEdges Time ) 2 o divStressTensor state - -
% Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
Modified: branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F 2012-09-26 20:56:32 UTC (rev 2169)
+++ branches/ocean_projects/strain_rate/src/core_ocean/mpas_ocn_mpas_core.F 2012-09-27 18:32:06 UTC (rev 2170)
@@ -5,6 +5,8 @@
use mpas_timekeeping
use mpas_dmpar
use mpas_timer
+ ! mrp added to test stress tensor:
+ use mpas_stress_tensor
use ocn_global_diagnostics
use ocn_test_cases
@@ -110,6 +112,13 @@
call ocn_init_z_level(domain)
+ ! mrp added just to test stress tensor ops:
+ print *, 'mpas_test_stress_tensor'
+ call mpas_test_stress_tensor(domain)
+ print *, 'mpas_test_stress_tensor complete'
+ stop
+ ! mrp added just to test stress tensor ops end
+
if (config_enforce_zstar_at_restart) then
call ocn_init_h_zstar(domain)
endif
@@ -569,9 +578,9 @@
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
referenceBottomDepth(1) = hZLevel(1)
- do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
- end do
+ !do k = 2,nVertLevels
+ ! referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
+ !end do
! TopOfCell needed where zero depth for the very top may be referenced.
referenceBottomDepthTopOfCell(1) = 0.0
Modified: branches/ocean_projects/strain_rate/src/operators/Makefile
===================================================================
--- branches/ocean_projects/strain_rate/src/operators/Makefile 2012-09-26 20:56:32 UTC (rev 2169)
+++ branches/ocean_projects/strain_rate/src/operators/Makefile 2012-09-27 18:32:06 UTC (rev 2170)
@@ -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_stress_tensor.o
all: operators
@@ -10,6 +10,7 @@
mpas_vector_reconstruction.o: mpas_rbf_interpolation.o
mpas_rbf_interpolation.o:
mpas_spline_interpolation:
+mpas_stress_tensor.o:
clean:
$(RM) *.o *.mod *.f90 libops.a
Added: branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F
===================================================================
--- branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F (rev 0)
+++ branches/ocean_projects/strain_rate/src/operators/mpas_stress_tensor.F 2012-09-27 18:32:06 UTC (rev 2170)
@@ -0,0 +1,532 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! mpas_stress_tensor
+!
+!> \brief MPAS strain rate and stress tensor computation
+!> \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_stress_tensor
+
+ use mpas_grid_types
+ use mpas_constants
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: mpas_strain_rate, &
+ mpas_div_stress_tensor, &
+ mpas_test_stress_tensor
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine mpas_strain_rate
+!
+!> \brief Computes normal and shear strain rate
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the normal strain rate at the cell centers
+!> and the shear strain rate at vertices using the weak derivative.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_strain_rate(grid, u,v, strainRateNormal, strainRateShear)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! 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) :: &
+ strainRateNormal !< Output: Normal component of strain rate tensor
+ real (kind=RKIND), dimension(:,:,:), intent(out) :: &
+ strainRateShear !< Output: Shear component of strain rate tensor
+
+ !-----------------------------------------------------------------
+ !
+ ! 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, &
+ cos1, sin1, cos2, sin2, 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
+
+ strainRateNormal = 0.0
+ strainRateShear = 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)
+
+ cos1 = cos(angleEdge(iEdge))
+ cos2 = cos1**2
+ sin1 = sin(angleEdge(iEdge))
+ sin2 = sin1**2
+ cossin = cos1*sin1
+
+ ! change later:
+ !do k=1,maxLevelEdgeBot(iEdge)
+ do k=1,nVertLevels
+ !eps11 = uDvEdge*cos2
+ !eps22 = uDvEdge*sin2
+ !eps12 = uDvEdge*cossin
+
+ ! using u and v both:
+ eps11 = u(k,iEdge)*cos2 - v(k,iEdge)*cossin
+ eps22 = u(k,iEdge)*sin2 + v(k,iEdge)*cossin
+ eps12 = u(k,iEdge)*cossin + v(k,iEdge)*0.5*(cos2-sin2)
+
+ !eps11 = u(k,iEdge)*cos2
+ !eps22 = u(k,iEdge)*sin2
+ !eps12 = u(k,iEdge)*cossin
+
+ strainRateNormal(1,k,cell1) = strainRateNormal(1,k,cell1) + eps11*invAreaCell1 * dvEdge(iEdge)
+ strainRateNormal(1,k,cell2) = strainRateNormal(1,k,cell2) - eps11*invAreaCell2 * dvEdge(iEdge)
+ strainRateNormal(2,k,cell1) = strainRateNormal(2,k,cell1) + eps22*invAreaCell1 * dvEdge(iEdge)
+ strainRateNormal(2,k,cell2) = strainRateNormal(2,k,cell2) - eps22*invAreaCell2 * dvEdge(iEdge)
+ strainRateNormal(3,k,cell1) = strainRateNormal(3,k,cell1) + eps12*invAreaCell1 * dvEdge(iEdge)
+ strainRateNormal(3,k,cell2) = strainRateNormal(3,k,cell2) - eps12*invAreaCell2 * dvEdge(iEdge)
+
+ ! using u and v both:
+ eps11 = -v(k,iEdge)*sin2 + u(k,iEdge)*cossin
+ eps22 = -v(k,iEdge)*cos2 - u(k,iEdge)*cossin
+ eps12 = v(k,iEdge)*cossin - u(k,iEdge)*0.5*(cos2-sin2)
+
+ strainRateShear(1,k,vertex1) = strainRateShear(1,k,vertex1) - eps11*invAreaTri1* dcEdge(iEdge)
+ strainRateShear(1,k,vertex2) = strainRateShear(1,k,vertex2) + eps11*invAreaTri2* dcEdge(iEdge)
+ strainRateShear(2,k,vertex1) = strainRateShear(2,k,vertex1) - eps22*invAreaTri1* dcEdge(iEdge)
+ strainRateShear(2,k,vertex2) = strainRateShear(2,k,vertex2) + eps22*invAreaTri2* dcEdge(iEdge)
+ strainRateShear(3,k,vertex1) = strainRateShear(3,k,vertex1) - eps12*invAreaTri1* dcEdge(iEdge)
+ strainRateShear(3,k,vertex2) = strainRateShear(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 ',strainRateNormal(1,:,6:7)
+ print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,6:7)
+ print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(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 ',strainRateShear(1,:,10:11)
+ print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,10:11)
+ print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,10:11)
+
+ ! 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 ',strainRateNormal(1,:,1:nCells)
+ ! print '(a,100f8.3)', 'eps Nor 2 ',strainRateNormal(2,:,1:nCells)
+ ! print '(a,100f8.3)', 'eps Nor 3 ',strainRateNormal(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 ',strainRateShear(1,:,1:nVertices)
+ ! print '(a,100f8.3)', 'eps Shr 2 ',strainRateShear(2,:,1:nVertices)
+ ! print '(a,100f8.3)', 'eps Shr 3 ',strainRateShear(3,:,1:nVertices)
+
+
+
+ end subroutine mpas_strain_rate!}}}
+
+!***********************************************************************
+!
+! routine mpas_div_stress_tensor
+!
+!> \brief Computes normal and shear strain rate
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the normal strain rate at the cell centers
+!> and the shear strain rate at vertices using the weak derivative.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_div_stress_tensor(grid, sigmaN, sigmaS, divStressTensor)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ sigmaN !< Input: Normal component of stress tensor
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ sigmaS !< Input: Shear component of stress tensor
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ divStressTensor !< Output: divergence of the stress tensor
+
+ !-----------------------------------------------------------------
+ !
+ ! 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, &
+ cos1, sin1, cos2, sin2, 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)
+
+ cos1 = cos(angleEdge(iEdge))
+ cos2 = cos1**2
+ sin1 = sin(angleEdge(iEdge))
+ sin2 = sin1**2
+ cossin = cos1*sin1
+
+ invdcEdge = 1.0 / dcEdge(iEdge)
+ invdvEdge = 1.0 / dvEdge(iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ divStressTensor(k,iEdge) = edgeMask(k,iEdge) * invdcEdge &
+ * ( cos2*sigmaN(1,k,cell2) - 2.0*cossin*sigmaN(3,k,cell2) + sin2*sigmaN(2,k,cell2) &
+ -(cos2*sigmaN(1,k,cell1) - 2.0*cossin*sigmaN(3,k,cell1) + sin2*sigmaN(2,k,cell1)) ) &
+ + edgeMask(k,iEdge) * invdvEdge &
+ * ( (cos2-sin2)*sigmaS(3,k,vertex2) + cossin*(sigmaS(1,k,vertex2) - sigmaS(2,k,vertex2)) &
+ -((cos2-sin2)*sigmaS(3,k,vertex1) + cossin*(sigmaS(1,k,vertex1) - sigmaS(2,k,vertex1))))
+ end do
+
+ end do
+
+ end subroutine mpas_div_stress_tensor!}}}
+
+!***********************************************************************
+!
+! routine mpas_test_stress_tensor
+!
+!> \brief Tests strain rate and stress tensor divergence operators
+!> \author Mark Petersen
+!> \date 25 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine tests strain rate and stress tensor divergence operators.
+!
+!-----------------------------------------------------------------------
+
+ subroutine mpas_test_stress_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
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnEdge
+
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, &
+ cos1, sin1, cos2, sin2, cossin, uDvEdge, uDcEdge, eps11, eps22, eps12, &
+ invdcEdge, invdvEdge, ux, uy, a, eoe, cn, cs, r, theta, rot
+ real (kind=RKIND), dimension(:), pointer :: angleEdge, &
+ xCell, yCell, xEdge, yEdge, xVertex, yVertex
+ real (kind=RKIND), dimension(:,:), pointer :: u, v, divStressTensor, weightsOnEdge
+ real (kind=RKIND), dimension(:,:,:), pointer :: strainRateNormal, strainRateShear
+
+
+ !-----------------------------------------------------------------
+ !
+ ! 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
+ strainRateNormal => block % state % time_levs(1) % state % strainRateNormal % array
+ strainRateShear => block % state % time_levs(1) % state % strainRateShear % array
+ divStressTensor => block % state % time_levs(1) % state % divStressTensor % array
+
+ cn = 1.5
+ cs = 0.7
+
+ do iEdge = 1,nEdges
+ a = angleEdge(iEdge)
+
+ !ux = cn*xEdge(iEdge)
+ !uy = cs*xEdge(iEdge)
+ !ux = -cs*yEdge(iEdge)
+ !uy = cn*yEdge(iEdge)
+ r = sqrt(xEdge(iEdge)**2 + yEdge(iEdge)**2)
+ theta = atan(yEdge(iEdge)/xEdge(iEdge))
+ rot = 4*pii/5
+
+ ux = r*cos(theta-rot)*(cn*cos(rot) - cs*sin(rot))
+ uy = r*cos(theta-rot)*(cn*sin(rot) + cs*cos(rot))
+
+ do k=1,nVertLevels
+ !do k = 1,maxLevelEdgeTop(iEdge)
+ u(k,iEdge) = ux*cos(a) + uy*sin(a)
+ end do
+ end do
+
+
+ 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, strainRateNormal, strainRateShear)
+
+ call mpas_div_stress_tensor(block % mesh, strainRateNormal, strainRateShear, divStressTensor)
+
+ 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 strainRateNormal1',minval(strainRateNormal(1,:,1:nCells)),maxval(strainRateNormal(1,:,1:nCells)),sum(strainRateNormal(1,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRateNormal2',minval(strainRateNormal(2,:,1:nCells)),maxval(strainRateNormal(2,:,1:nCells)),sum(strainRateNormal(2,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRateNormal3',minval(strainRateNormal(3,:,1:nCells)),maxval(strainRateNormal(3,:,1:nCells)),sum(strainRateNormal(3,:,1:nCells))/nCells/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRateShear1 ',minval(strainRateShear(1,:,1:nVertices)),maxval(strainRateShear(1,:,1:nVertices)),sum(strainRateShear(1,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRateShear2 ',minval(strainRateShear(2,:,1:nVertices)),maxval(strainRateShear(2,:,1:nVertices)),sum(strainRateShear(2,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean strainRateShear3 ',minval(strainRateShear(3,:,1:nVertices)),maxval(strainRateShear(3,:,1:nVertices)),sum(strainRateShear(3,:,1:nVertices))/nVertices/nVertLevels
+ print '(a,3f20.10)', 'min max mean divStressTensor ',minval(divStressTensor(:,1:nEdges)),maxval(divStressTensor(:,1:nEdges)),sum(divStressTensor(:,1:nEdges))/nEdges/nVertLevels
+
+ block => block % next
+ end do
+
+ end subroutine mpas_test_stress_tensor!}}}
+
+!***********************************************************************
+!
+! routine mpas_stress_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_stress_tensor_init(err)!{{{
+! integer, intent(out) :: err !< Output: Error flag
+
+! end subroutine mpas_stress_tensor_init!}}}
+
+!***********************************************************************
+
+end module mpas_stress_tensor
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
</font>
</pre>