<p><b>dwj07@fsu.edu</b> 2011-09-23 15:26:14 -0600 (Fri, 23 Sep 2011)</p><p><br>
        -- BRANCH COMMIT --<br>
<br>
        Splitting tendency subroutines into their own module to prepare for timestepping modules.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/Makefile        2011-09-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/Makefile        2011-09-23 21:26:14 UTC (rev 1019)
@@ -35,6 +35,7 @@
            module_OcnVmixCoefsTanh.o \
            module_OcnRestoring.o \
            module_OcnEquationOfState.o \
+           module_OcnTendency.o \
        module_time_integration.o \
        module_global_diagnostics.o
 
@@ -48,39 +49,42 @@
 
 module_advection.o:
 
-module_time_integration.o:  module_OcnThickHadv.o \
-                                                        module_OcnThickVadv.o \
-                                                        module_OcnVelCoriolis.o \
-                                                        module_OcnVelVadv.o \
-                                                        module_OcnVelHmix.o \
-                                                        module_OcnVelHmixDel2.o \
-                                                        module_OcnVelHmixDel4.o \
-                                                        module_OcnVelForcing.o \
-                                                        module_OcnVelForcingWindStress.o \
-                                                        module_OcnVelForcingBottomDrag.o \
-                                                        module_OcnVelPressureGrad.o \
-                                                        module_OcnTracerVadv.o \
-                                                        module_OcnTracerVadvSpline.o \
-                                                        module_OcnTracerVadvSpline2.o \
-                                                        module_OcnTracerVadvSpline3.o \
-                                                        module_OcnTracerVadvStencil.o \
-                                                        module_OcnTracerVadvStencil2.o \
-                                                        module_OcnTracerVadvStencil3.o \
-                                                        module_OcnTracerVadvStencil4.o \
-                                                        module_OcnTracerHadv.o \
-                                                        module_OcnTracerHadv2.o \
-                                                        module_OcnTracerHadv3.o \
-                                                        module_OcnTracerHadv4.o \
-                                                        module_OcnTracerHmix.o \
-                                                        module_OcnTracerHmixDel2.o \
-                                                        module_OcnTracerHmixDel4.o \
-                                                        module_OcnVmix.o \
-                                                        module_OcnVmixCoefsConst.o \
-                                                        module_OcnVmixCoefsRich.o \
-                                                        module_OcnVmixCoefsTanh.o \
-                                                        module_OcnRestoring.o \
-                                                        module_OcnEquationOfState.o
+module_time_integration.o: module_OcnThickHadv.o \
+                               module_OcnThickVadv.o \
+                                                   module_OcnVelCoriolis.o \
+                                                   module_OcnVelVadv.o \
+                                                   module_OcnVelHmix.o \
+                                                   module_OcnVelHmixDel2.o \
+                                                   module_OcnVelHmixDel4.o \
+                                                   module_OcnVelForcing.o \
+                                                   module_OcnVelForcingWindStress.o \
+                                                   module_OcnVelForcingBottomDrag.o \
+                                                   module_OcnVelPressureGrad.o \
+                                                   module_OcnTracerVadv.o \
+                                                   module_OcnTracerVadvSpline.o \
+                                                   module_OcnTracerVadvSpline2.o \
+                                                   module_OcnTracerVadvSpline3.o \
+                                                   module_OcnTracerVadvStencil.o \
+                                                   module_OcnTracerVadvStencil2.o \
+                                                   module_OcnTracerVadvStencil3.o \
+                                                   module_OcnTracerVadvStencil4.o \
+                                                   module_OcnTracerHadv.o \
+                                                   module_OcnTracerHadv2.o \
+                                                   module_OcnTracerHadv3.o \
+                                                   module_OcnTracerHadv4.o \
+                                                   module_OcnTracerHmix.o \
+                                                   module_OcnTracerHmixDel2.o \
+                                                   module_OcnTracerHmixDel4.o \
+                                                   module_OcnVmix.o \
+                                                   module_OcnVmixCoefsConst.o \
+                                                   module_OcnVmixCoefsRich.o \
+                                                   module_OcnVmixCoefsTanh.o \
+                                                   module_OcnRestoring.o \
+                                                   module_OcnEquationOfState.o \
+                                                   module_OcnTendency.o 
 
+module_OcnTendency.o:
+
 module_global_diagnostics.o: 
 
 module_OcnThickHadv.o:
@@ -182,6 +186,7 @@
                                         module_OcnVmixCoefsRich.o \
                                         module_OcnVmixCoefsTanh.o \
                                         module_OcnEquationOfState.o \
+                                        module_OcnTendency.o \
                                         module_time_integration.o
 
 clean:

Added: branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F                                (rev 0)
+++ branches/ocean_projects/performance/src/core_ocean/module_OcnTendency.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -0,0 +1,1380 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  OcnTendency
+!
+!&gt; \brief MPAS ocean tendency driver
+!&gt; \author Doug Jacobsen
+!&gt; \date   23 September 2011
+!&gt; \version SVN:$Id:$
+!&gt; \details
+!&gt;  This module contains the routines for computing
+!&gt;  various tendencies for the ocean.
+!
+!-----------------------------------------------------------------------
+
+module OcnTendency
+
+   use grid_types
+   use configure
+   use constants
+   use timer
+
+   use OcnThickHadv
+   use OcnThickVadv
+
+   use OcnVelCoriolis
+   use OcnVelPressureGrad
+   use OcnVelVadv
+   use OcnVelHmix
+   use OcnVelForcing
+
+   use OcnTracerHadv
+   use OcnTracerVadv
+   use OcnTracerHmix
+   use OcnRestoring
+
+   use OcnEquationOfState
+   use OcnVmix
+
+
+   implicit none
+   private
+   save
+
+   !--------------------------------------------------------------------
+   !
+   ! Public parameters
+   !
+   !--------------------------------------------------------------------
+
+   !--------------------------------------------------------------------
+   !
+   ! Public member functions
+   !
+   !--------------------------------------------------------------------
+
+   public :: OcnTendH, &amp;
+             OcnTendU, &amp;
+             OcnTendScalar, &amp;
+             OcnDiagnosticSolve, &amp;
+             OcnWtop, &amp;
+             OcnFUPerp
+
+   !--------------------------------------------------------------------
+   !
+   ! Private module variables
+   !
+   !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+!  routine OcnTendH
+!
+!&gt; \brief   Computes thickness tendency
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes the thickness tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnTendH(tend, s, d, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Compute height and normal wind tendencies, as well as diagnostic variables
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !
+   ! Output: tend - computed tendencies for prognostic variables
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (tend_type), intent(inout) :: tend
+      type (state_type), intent(in) :: s
+      type (diagnostics_type), intent(in) :: d
+      type (mesh_type), intent(in) :: grid
+
+      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
+        vertex1, vertex2, eoe, i, j, err
+
+! mrp 110512 I just split compute_tend into compute_tend_u and OcnTendH.
+!  Most of these variables can be removed, but at a later time.
+      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
+        upstream_bias, wTopEdge, rho0Inv, r
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
+        zMidZLevel, zTopZLevel 
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &amp;
+        tend_h, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
+        MontPot, wTop, divergence, vertViscTopOfEdge
+      type (dm_info) :: dminfo
+
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+      integer, dimension(:,:), pointer :: &amp;
+        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
+        edgesOnEdge, edgesOnVertex
+      real (kind=RKIND) :: u_diffusion
+      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+      call timer_start(&quot;OcnTendH&quot;)
+
+      h           =&gt; s % h % array
+      u           =&gt; s % u % array
+      v           =&gt; s % v % array
+      wTop        =&gt; s % wTop % array
+      h_edge      =&gt; s % h_edge % array
+      circulation =&gt; s % circulation % array
+      vorticity   =&gt; s % vorticity % array
+      divergence  =&gt; s % divergence % array
+      ke          =&gt; s % ke % array
+      ke_edge     =&gt; s % ke_edge % array
+      pv_edge     =&gt; s % pv_edge % array
+      MontPot     =&gt; s % MontPot % array
+      pressure    =&gt; s % pressure % array
+      vertViscTopOfEdge =&gt; d % vertViscTopOfEdge % array
+
+      weightsOnEdge     =&gt; grid % weightsOnEdge % array
+      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      cellsOnVertex     =&gt; grid % cellsOnVertex % array
+      verticesOnEdge    =&gt; grid % verticesOnEdge % array
+      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
+      edgesOnCell       =&gt; grid % edgesOnCell % array
+      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
+      edgesOnEdge       =&gt; grid % edgesOnEdge % array
+      edgesOnVertex     =&gt; grid % edgesOnVertex % array
+      dcEdge            =&gt; grid % dcEdge % array
+      dvEdge            =&gt; grid % dvEdge % array
+      areaCell          =&gt; grid % areaCell % array
+      areaTriangle      =&gt; grid % areaTriangle % array
+      h_s               =&gt; grid % h_s % array
+      fVertex           =&gt; grid % fVertex % array
+      fEdge             =&gt; grid % fEdge % array
+      zMidZLevel        =&gt; grid % zMidZLevel % array
+      zTopZLevel        =&gt; grid % zTopZLevel % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
+      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
+
+      tend_h      =&gt; tend % h % array
+                  
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+      nEdgesSolve = grid % nEdgesSolve
+      nVertices   = grid % nVertices
+      nVertLevels = grid % nVertLevels
+
+      !
+      ! height tendency: start accumulating tendency terms
+      !
+      tend_h = 0.0
+
+      !
+      ! height tendency: horizontal advection term -</font>
<font color="blue">abla\cdot ( hu)
+      !
+      ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3. 
+      ! for explanation of divergence operator.
+      !
+      ! for z-level, only compute height tendency for top layer.
+
+      call timer_start(&quot;OcnTendH-horiz adv&quot;)
+
+      call OcnThickHadvTend(grid, u, h_edge, tend_h, err)
+
+      call timer_stop(&quot;OcnTendH-horiz adv&quot;)
+
+      !
+      ! height tendency: vertical advection term -d/dz(hw)
+      !
+      ! Vertical advection computed for top layer of a z grid only.
+      call timer_start(&quot;OcnTendH-vert adv&quot;)
+
+      call OcnThickVadvTend(grid, wTop, tend_h, err)
+
+      call timer_stop(&quot;OcnTendH-vert adv&quot;)
+      call timer_stop(&quot;OcnTendH&quot;)
+   
+   end subroutine OcnTendH!}}}
+
+!***********************************************************************
+!
+!  routine OcnTendU
+!
+!&gt; \brief   Computes velocity tendency
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes the velocity tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnTendU(tend, s, d, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Compute height and normal wind tendencies, as well as diagnostic variables
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !
+   ! Output: tend - computed tendencies for prognostic variables
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (tend_type), intent(inout) :: tend
+      type (state_type), intent(in) :: s
+      type (diagnostics_type), intent(in) :: d
+      type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into OcnTendU and compute_tend_h.
+!  Some of these variables can be removed, but at a later time.
+      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
+        vertex1, vertex2, eoe, i, j
+
+      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, err
+      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
+        upstream_bias, wTopEdge, rho0Inv, r, visc_vorticity_coef
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
+        zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &amp;
+        tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
+        MontPot, wTop, divergence, vertViscTopOfEdge
+      type (dm_info) :: dminfo
+
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+      integer, dimension(:,:), pointer :: &amp;
+        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
+        edgesOnEdge, edgesOnVertex
+      real (kind=RKIND) :: u_diffusion
+      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+      real (kind=RKIND), dimension(:,:), pointer :: u_src
+      real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+      call timer_start(&quot;OcnTendU&quot;)
+
+      h           =&gt; s % h % array
+      u           =&gt; s % u % array
+      v           =&gt; s % v % array
+      wTop        =&gt; s % wTop % array
+      h_edge      =&gt; s % h_edge % array
+      circulation =&gt; s % circulation % array
+      vorticity   =&gt; s % vorticity % array
+      divergence  =&gt; s % divergence % array
+      ke          =&gt; s % ke % array
+      ke_edge     =&gt; s % ke_edge % array
+      pv_edge     =&gt; s % pv_edge % array
+      MontPot     =&gt; s % MontPot % array
+      pressure    =&gt; s % pressure % array
+      vertViscTopOfEdge =&gt; d % vertViscTopOfEdge % array
+
+      weightsOnEdge     =&gt; grid % weightsOnEdge % array
+      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      cellsOnVertex     =&gt; grid % cellsOnVertex % array
+      verticesOnEdge    =&gt; grid % verticesOnEdge % array
+      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
+      edgesOnCell       =&gt; grid % edgesOnCell % array
+      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
+      edgesOnEdge       =&gt; grid % edgesOnEdge % array
+      edgesOnVertex     =&gt; grid % edgesOnVertex % array
+      dcEdge            =&gt; grid % dcEdge % array
+      dvEdge            =&gt; grid % dvEdge % array
+      areaCell          =&gt; grid % areaCell % array
+      areaTriangle      =&gt; grid % areaTriangle % array
+      h_s               =&gt; grid % h_s % array
+! mrp 110516 cleanup fvertex fedge not used in this subroutine
+      fVertex           =&gt; grid % fVertex % array
+      fEdge             =&gt; grid % fEdge % array
+      zMidZLevel        =&gt; grid % zMidZLevel % array
+      zTopZLevel        =&gt; grid % zTopZLevel % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
+      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
+
+      tend_u      =&gt; tend % u % array
+                  
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+      nEdgesSolve = grid % nEdgesSolve
+      nVertices   = grid % nVertices
+      nVertLevels = grid % nVertLevels
+
+      u_src =&gt; grid % u_src % array
+
+      meshScalingDel2 =&gt; grid % meshScalingDel2 % array
+      meshScalingDel4 =&gt; grid % meshScalingDel4 % array
+
+      !
+      ! velocity tendency: start accumulating tendency terms
+      !
+      ! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
+      tend_u(:,:) = 0.0
+
+      !
+      ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
+      !
+
+      call timer_start(&quot;OcnTendU-coriolis&quot;)
+
+!     call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
+
+      do iEdge=1,grid % nEdgesSolve
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+
+         do k=1,maxLevelEdgeTop(iEdge)
+
+            q = 0.0
+            do j = 1,nEdgesOnEdge(iEdge)
+               eoe = edgesOnEdge(j,iEdge)
+               workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
+               q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe) 
+            end do
+
+           tend_u(k,iEdge) = tend_u(k,iEdge)     &amp;
+                  + q     &amp;
+                  - (   ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
+         end do
+      end do
+
+      call timer_stop(&quot;OcnTendU-coriolis&quot;)
+
+      !
+      ! velocity tendency: vertical advection term -w du/dz
+      !
+      call timer_start(&quot;OcnTendU-vert adv&quot;)
+
+      call OcnVelVadvTend(grid, u, wTop, tend_u, err)
+
+      call timer_stop(&quot;OcnTendU-vert adv&quot;)
+
+      !
+      ! velocity tendency: pressure gradient
+      !
+      call timer_start(&quot;OcnTendU-pressure grad&quot;)
+
+      if (config_vert_grid_type.eq.'isopycnal') then
+          call OcnVelPressureGradTend(grid, MontPot, tend_u, err)
+      elseif (config_vert_grid_type.eq.'zlevel') then
+          call OcnVelPressureGradTend(grid, pressure, tend_u, err)
+      end if
+
+      call timer_stop(&quot;OcnTendU-pressure grad&quot;)
+
+      !
+      ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="blue">abla^2 u
+      !   computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity )
+      !   strictly only valid for config_h_mom_eddy_visc2 == constant
+      !
+      call timer_start(&quot;OcnTendU-horiz mix&quot;)
+
+      call OcnVelHmixTend(grid, divergence, vorticity, tend_u, err)
+
+      call timer_stop(&quot;OcnTendU-horiz mix&quot;)
+
+      !
+      ! velocity tendency: forcing and bottom drag
+      !
+      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
+      ! know the bottom edge with nonzero velocity and place the drag there.
+
+      call timer_start(&quot;OcnTendU-forcings&quot;)
+
+      call OcnVelForcingTend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
+
+      call timer_stop(&quot;OcnTendU-forcings&quot;)
+
+      !
+      ! velocity tendency: vertical mixing d/dz( nu_v du/dz))
+      !
+      if (.not.config_implicit_vertical_mix) then
+          call timer_start(&quot;OcnTendU-explicit vert mix&quot;)
+          call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
+!        allocate(fluxVertTop(nVertLevels+1))
+!        fluxVertTop(1) = 0.0
+!        do iEdge=1,grid % nEdgesSolve
+!        
+!           do k=2,maxLevelEdgeTop(iEdge)
+!             fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &amp;
+!                * ( u(k-1,iEdge) - u(k,iEdge) ) &amp;
+!                * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
+!           enddo
+!           fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+
+!           do k=1,maxLevelEdgeTop(iEdge)
+!             tend_u(k,iEdge) = tend_u(k,iEdge) &amp;
+!               + (fluxVertTop(k) - fluxVertTop(k+1)) &amp;
+!               / h_edge(k,iEdge)
+!           enddo
+
+!        end do
+!        deallocate(fluxVertTop)
+          call timer_stop(&quot;OcnTendU-explicit vert mix&quot;)
+      endif
+      call timer_stop(&quot;OcnTendU&quot;)
+
+   end subroutine OcnTendU!}}}
+
+!***********************************************************************
+!
+!  routine OcnTendSalar
+!
+!&gt; \brief   Computes scalar tendency
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes the scalar (tracer) tendency for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnTendScalar(tend, s, d, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !        note: the variable s % tracers really contains the tracers, 
+   !              not tracers*h
+   !
+   ! Output: tend - computed scalar tendencies
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (tend_type), intent(inout) :: tend
+      type (state_type), intent(in) :: s
+      type (diagnostics_type), intent(in) :: d
+      type (mesh_type), intent(in) :: grid
+
+      integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&amp;
+        nEdges, nCells, nCellsSolve, nVertLevels, num_tracers, err
+      real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+      real (kind=RKIND) :: flux, tracer_edge, r
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        u,h,wTop, h_edge, vertDiffTopOfCell
+      real (kind=RKIND), dimension(:,:,:), pointer :: &amp;
+        tracers, tend_tr
+      integer, dimension(:,:), pointer :: boundaryEdge
+      type (dm_info) :: dminfo
+
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+      integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
+      real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &amp;
+         hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
+      real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &amp;
+            posZTopZLevel
+      real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
+      real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
+
+
+      real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+      real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
+
+      integer :: index_temperature, index_salinity, rrr
+      real (kind=RKIND), dimension(:), pointer :: temperatureRestore, salinityRestore
+
+      call timer_start(&quot;OcnTendScalar&quot;)
+
+      u           =&gt; s % u % array
+      h           =&gt; s % h % array
+      boundaryCell=&gt; grid % boundaryCell % array
+      wTop        =&gt; s % wTop % array
+      tracers     =&gt; s % tracers % array
+      h_edge      =&gt; s % h_edge % array
+      vertDiffTopOfCell =&gt; d % vertDiffTopOfCell % array
+
+      tend_tr     =&gt; tend % tracers % array
+                  
+      areaCell          =&gt; grid % areaCell % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      dvEdge            =&gt; grid % dvEdge % array
+      dcEdge            =&gt; grid % dcEdge % array
+      zTopZLevel        =&gt; grid % zTopZLevel % array
+      zMidZLevel        =&gt; grid % zMidZLevel % array
+      hRatioZLevelK    =&gt; grid % hRatioZLevelK % array
+      hRatioZLevelKm1    =&gt; grid % hRatioZLevelKm1 % array
+      boundaryEdge      =&gt; grid % boundaryEdge % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
+      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
+
+      nEdges      = grid % nEdges
+      nCells      = grid % nCells
+      nCellsSolve = grid % nCellsSolve
+      nVertLevels = grid % nVertLevels
+      num_tracers = s % num_tracers
+
+      meshScalingDel2 =&gt; grid % meshScalingDel2 % array
+      meshScalingDel4 =&gt; grid % meshScalingDel4 % array
+
+
+      deriv_two   =&gt; grid % deriv_two % array
+
+      !
+      ! initialize tracer tendency (RHS of tracer equation) to zero.
+      !
+      tend_tr(:,:,:) = 0.0
+
+      !
+      ! tracer tendency: horizontal advection term -div( h \phi u)
+      !
+      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
+      ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
+      ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
+      ! tracer_edge at the boundary will also need to be defined for flux boundaries.
+
+      call timer_start(&quot;OcnTendScalar-horiz adv&quot;)
+
+      call OcnTracerHadvTend(grid, u, h_edge, tracers, tend_tr, err)
+
+      call timer_stop(&quot;OcnTendScalar-horiz adv&quot;)
+
+
+      !
+      ! tracer tendency: vertical advection term -d/dz( h \phi w)
+      !
+
+      call timer_start(&quot;OcnTendScalar-vert adv&quot;)
+
+      call OcnTracerVadvTend(grid, wTop, tracers, tend_tr, err)
+
+      call timer_stop(&quot;OcnTendScalar-vert adv&quot;)
+
+      !
+      ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
+      !
+      call timer_start(&quot;OcnTendScalar-horiz diff&quot;)
+
+      call OcnTracerHmixTend(grid, h_edge, tracers, tend_tr, err)
+
+      call timer_stop(&quot;OcnTendScalar-horiz diff&quot;)
+
+! mrp 110516 printing
+!print *, 'tend_tr 1',minval(tend_tr(3,1,1:nCells)),&amp;
+!                   maxval(tend_tr(3,1,1:nCells))
+!print *, 'tracer  1',minval(tracers(3,1,1:nCells)),&amp;
+!                   maxval(tracers(3,1,1:nCells))
+! mrp 110516 printing end
+
+      !
+      ! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
+      !
+      if (.not.config_implicit_vertical_mix) then
+         call timer_start(&quot;OcnTendScalar-explicit vert diff&quot;)
+         call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
+!        allocate(fluxVertTop(num_tracers,nVertLevels+1))
+!        fluxVertTop(:,1) = 0.0
+!        do iCell=1,nCellsSolve 
+
+!           do k=2,maxLevelCell(iCell)
+!             do iTracer=1,num_tracers
+!               ! compute \kappa_v d\phi/dz
+!               fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &amp;
+!                  * (   tracers(iTracer,k-1,iCell)    &amp;
+!                      - tracers(iTracer,k  ,iCell) )  &amp;
+!                  * 2 / (h(k-1,iCell) + h(k,iCell))
+
+!             enddo
+!           enddo
+!           fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
+
+!           do k=1,maxLevelCell(iCell)
+!             do iTracer=1,num_tracers
+!               ! This is h d/dz( fluxVertTop) but h and dz cancel, so 
+!               ! reduces to delta( fluxVertTop)
+!               tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &amp;
+!                 + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
+!             enddo
+!           enddo
+!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
+!print '(a,50e12.2)', 'tend_tr    ',tend_tr(3,1,1:maxLevelCell(iCell))
+
+!        enddo ! iCell loop
+!        deallocate(fluxVertTop)
+         call timer_stop(&quot;OcnTendScalar-explicit vert diff&quot;)
+      endif
+
+! mrp 110516 printing
+!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&amp;
+!                   maxval(tend_tr(3,1,1:nCells))
+! mrp 110516 printing end
+
+      !
+      ! add restoring to T and S in top model layer
+      !
+      call timer_start(&quot;OcnTendScalar-restoring&quot;)
+
+      call OcnRestoringTend(grid, h, s%index_temperature, s%index_salinity, tracers, tend_tr, err)
+
+      call timer_stop(&quot;OcnTendScalar-restoring&quot;)
+
+ 10   format(2i8,10e20.10)
+      call timer_stop(&quot;OcnTendScalar&quot;)
+
+   end subroutine OcnTendScalar!}}}
+
+!***********************************************************************
+!
+!  routine OcnDiagnosticSolve
+!
+!&gt; \brief   Computes diagnostic variables
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes the diagnostic variables for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnDiagnosticSolve(dt, s, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Compute diagnostic fields used in the tendency computations
+   !
+   ! Input: grid - grid metadata
+   !
+   ! Output: s - computed diagnostics
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      real (kind=RKIND), intent(in) :: dt
+      type (state_type), intent(inout) :: s
+      type (mesh_type), intent(in) :: grid
+
+
+      integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+
+      integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
+
+
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
+        hZLevel
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&amp;
+        circulation, vorticity, ke, ke_edge, MontPot, wTop, &amp;
+        pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &amp;
+        rho, temperature, salinity
+      real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+      real (kind=RKIND), dimension(:), allocatable:: pTop
+      real (kind=RKIND), dimension(:,:), allocatable:: div_u
+      character :: c1*6
+
+      integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &amp;
+        verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &amp;
+        boundaryEdge, boundaryCell
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &amp;
+        maxLevelVertexBot,  maxLevelVertexTop
+      real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+      real (kind=RKIND) :: coef_3rd_order
+      real (kind=RKIND) :: r, h1, h2
+
+      call timer_start(&quot;OcnDiagnosticSolve&quot;)
+
+      h           =&gt; s % h % array
+      u           =&gt; s % u % array
+      v           =&gt; s % v % array
+      wTop        =&gt; s % wTop % array
+      h_edge      =&gt; s % h_edge % array
+      circulation =&gt; s % circulation % array
+      vorticity   =&gt; s % vorticity % array
+      divergence  =&gt; s % divergence % array
+      ke          =&gt; s % ke % array
+      ke_edge     =&gt; s % ke_edge % array
+      pv_edge     =&gt; s % pv_edge % array
+      pv_vertex   =&gt; s % pv_vertex % array
+      pv_cell     =&gt; s % pv_cell % array
+      gradPVn     =&gt; s % gradPVn % array
+      gradPVt     =&gt; s % gradPVt % array
+      rho         =&gt; s % rho % array
+      tracers     =&gt; s % tracers % array
+      MontPot     =&gt; s % MontPot % array
+      pressure    =&gt; s % pressure % array
+
+      weightsOnEdge     =&gt; grid % weightsOnEdge % array
+      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      cellsOnVertex     =&gt; grid % cellsOnVertex % array
+      verticesOnEdge    =&gt; grid % verticesOnEdge % array
+      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
+      edgesOnCell       =&gt; grid % edgesOnCell % array
+      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
+      edgesOnEdge       =&gt; grid % edgesOnEdge % array
+      edgesOnVertex     =&gt; grid % edgesOnVertex % array
+      dcEdge            =&gt; grid % dcEdge % array
+      dvEdge            =&gt; grid % dvEdge % array
+      areaCell          =&gt; grid % areaCell % array
+      areaTriangle      =&gt; grid % areaTriangle % array
+      h_s               =&gt; grid % h_s % array
+      fVertex           =&gt; grid % fVertex % array
+      fEdge             =&gt; grid % fEdge % array
+      hZLevel           =&gt; grid % hZLevel % array
+      deriv_two         =&gt; grid % deriv_two % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeTop   =&gt; grid % maxLevelEdgeTop % array
+      maxLevelEdgeBot   =&gt; grid % maxLevelEdgeBot % array
+      maxLevelVertexBot =&gt; grid % maxLevelVertexBot % array
+      maxLevelVertexTop =&gt; grid % maxLevelVertexTop % array
+                  
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+      nVertices   = grid % nVertices
+      nVertLevels = grid % nVertLevels
+      vertexDegree = grid % vertexDegree
+
+      boundaryEdge =&gt; grid % boundaryEdge % array
+      boundaryCell =&gt; grid % boundaryCell % array
+
+      !
+      ! Compute height on cell edges at velocity locations
+      !   Namelist options control the order of accuracy of the reconstructed h_edge value
+      !
+      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
+      ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
+
+      ! mrp 110516 efficiency note: For z-level, only do this on level 1.  h_edge for all
+      ! lower levels is defined by hZlevel.
+
+      call timer_start(&quot;OcnDiagnosticSolve-hEdge&quot;)
+
+      coef_3rd_order = 0.
+      if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
+      if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+      if (config_thickness_adv_order == 2) then
+          call timer_start(&quot;OcnDiagnosticSolve-hEdge 2&quot;)
+
+         do iEdge=1,nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+            do k=1,maxLevelEdgeTop(iEdge)
+               h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+            end do
+         end do
+          call timer_stop(&quot;OcnDiagnosticSolve-hEdge 2&quot;)
+
+      else if (config_thickness_adv_order == 3) then
+          call timer_start(&quot;OcnDiagnosticSolve-hEdge 3&quot;)
+
+         do iEdge=1,nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+
+            do k=1,maxLevelEdgeTop(iEdge)
+
+               d2fdx2_cell1 = 0.0
+               d2fdx2_cell2 = 0.0
+
+               !-- if not a boundary cell
+               if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+                  d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+                  d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+                  !-- all edges of cell 1
+                  do i=1, grid % nEdgesOnCell % array (cell1)
+                          d2fdx2_cell1 = d2fdx2_cell1 + &amp;
+                          deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+                  end do
+
+                  !-- all edges of cell 2
+                  do i=1, grid % nEdgesOnCell % array (cell2)
+                          d2fdx2_cell2 = d2fdx2_cell2 + &amp;
+                          deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+                  end do
+
+               endif
+
+               !-- if u &gt; 0:
+               if (u(k,iEdge) &gt; 0) then
+                  h_edge(k,iEdge) =     &amp;
+                       0.5*(h(k,cell1) + h(k,cell2))      &amp;
+                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &amp;
+                       -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+               !-- else u &lt;= 0:
+               else
+                  h_edge(k,iEdge) =     &amp;
+                       0.5*(h(k,cell1) + h(k,cell2))      &amp;
+                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &amp;
+                       +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+               end if
+
+            end do   ! do k
+         end do         ! do iEdge
+
+          call timer_stop(&quot;OcnDiagnosticSolve-hEdge 3&quot;)
+      else  if (config_thickness_adv_order == 4) then
+          call timer_start(&quot;OcnDiagnosticSolve-hEdge 4&quot;)
+
+         do iEdge=1,nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+
+            do k=1,maxLevelEdgeTop(iEdge)
+
+               d2fdx2_cell1 = 0.0
+               d2fdx2_cell2 = 0.0
+
+               !-- if not a boundary cell
+               if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+                  d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+                  d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+                  !-- all edges of cell 1
+                  do i=1, grid % nEdgesOnCell % array (cell1)
+                          d2fdx2_cell1 = d2fdx2_cell1 + &amp;
+                          deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+                  end do
+
+                  !-- all edges of cell 2
+                  do i=1, grid % nEdgesOnCell % array (cell2)
+                          d2fdx2_cell2 = d2fdx2_cell2 + &amp;
+                          deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+                  end do
+
+               endif
+
+               h_edge(k,iEdge) =   &amp;
+                    0.5*(h(k,cell1) + h(k,cell2))      &amp;
+                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+            end do   ! do k
+         end do         ! do iEdge
+
+         call timer_stop(&quot;OcnDiagnosticSolve-hEdge 4&quot;)
+      endif   ! if(config_thickness_adv_order == 2)
+      call timer_stop(&quot;OcnDiagnosticSolve-hEdge&quot;)
+
+      !
+      ! set the velocity and height at dummy address
+      !    used -1e34 so error clearly occurs if these values are used.
+      !
+!mrp 110516 change to zero, change back later:
+      u(:,nEdges+1) = -1e34
+      h(:,nCells+1) = -1e34
+      tracers(s % index_temperature,:,nCells+1) = -1e34
+      tracers(s % index_salinity,:,nCells+1) = -1e34
+
+      !
+      ! Compute circulation and relative vorticity at each vertex
+      !
+      circulation(:,:) = 0.0
+      do iEdge=1,nEdges
+         vertex1 = verticesOnEdge(1,iEdge)
+         vertex2 = verticesOnEdge(2,iEdge)
+         do k=1,maxLevelEdgeBot(iEdge)
+            circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
+            circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
+         end do
+      end do
+      do iVertex=1,nVertices
+         do k=1,maxLevelVertexBot(iVertex)
+            vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+         end do
+      end do
+
+      !
+      ! Compute the divergence at each cell center
+      !
+      divergence(:,:) = 0.0
+      do iEdge=1,nEdges
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         do k=1,maxLevelEdgeBot(iEdge)
+             divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+             divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+         enddo
+      end do
+      do iCell = 1,nCells
+         r = 1.0 / areaCell(iCell)
+         do k = 1,maxLevelCell(iCell)
+            divergence(k,iCell) = divergence(k,iCell) * r
+         enddo
+      enddo
+
+      !
+      ! Compute kinetic energy in each cell
+      !
+      ke(:,:) = 0.0
+      do iEdge=1,nEdges
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         do k=1,maxLevelEdgeBot(iEdge)
+              ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+              ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+         enddo
+      end do
+      do iCell = 1,nCells
+         do k = 1,maxLevelCell(iCell)
+            ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
+         enddo
+      enddo
+
+      !
+      ! Compute v (tangential) velocities
+      !
+      v(:,:) = 0.0
+      do iEdge = 1,nEdges
+         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,maxLevelEdgeTop(iEdge) 
+               v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+            end do
+         end do
+      end do
+
+      !
+      ! Compute ke on cell edges at velocity locations for quadratic bottom drag. 
+      !
+      ! mrp 101025 efficiency note: we could get rid of ke_edge completely by 
+      ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+      ke_edge = 0.0  !mrp remove 0 for efficiency
+      do iEdge=1,nEdges
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         do k=1,maxLevelEdgeTop(iEdge)
+            ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+         end do
+      end do
+
+      !
+      ! Compute height at vertices, pv at vertices, and average pv to edge locations
+      !  ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
+      !
+      if (trim(config_time_integration) == 'RK4') then
+         ! for RK4, PV is really PV = (eta+f)/h
+         fCoef = 1
+      elseif (trim(config_time_integration) == 'split_explicit' &amp;
+          .or.trim(config_time_integration) == 'unsplit_explicit') then
+         ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
+! mrp temp, new should be:
+         fCoef = 0
+! old, for testing:
+!         fCoef = 1
+      end if
+
+      do iVertex = 1,nVertices
+         do k=1,maxLevelVertexBot(iVertex)
+            h_vertex = 0.0
+            do i=1,vertexDegree
+               h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+            end do
+            h_vertex = h_vertex / areaTriangle(iVertex)
+
+            pv_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+         end do
+      end do
+
+      !
+      ! Compute pv at cell centers
+      !    ( this computes pv_cell for all real cells and distance-1 ghost cells )
+      !
+      pv_cell(:,:) = 0.0
+      do iVertex = 1,nVertices
+         do i=1,vertexDegree
+            iCell = cellsOnVertex(i,iVertex)
+            do k = 1,maxLevelCell(iCell)
+               pv_cell(k,iCell) = pv_cell(k,iCell)  &amp;
+                  + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &amp;
+                    / areaCell(iCell)
+            enddo
+         enddo
+      enddo
+
+      !
+      ! Compute pv at the edges
+      !   ( this computes pv_edge at all edges bounding real cells )
+      !
+      pv_edge(:,:) = 0.0
+      do iVertex = 1,nVertices
+         do i=1,vertexDegree
+            iEdge = edgesOnVertex(i,iVertex)
+            do k=1,maxLevelEdgeBot(iEdge)
+               pv_edge(k,iEdge) =  pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+            enddo
+        end do
+      end do
+
+      !
+      ! Compute gradient of PV in normal direction
+      !   ( this computes gradPVn for all edges bounding real cells )
+      !
+      gradPVn(:,:) = 0.0
+      do iEdge = 1,nEdges
+         do k=1,maxLevelEdgeTop(iEdge)
+            gradPVn(k,iEdge) = (  pv_cell(k,cellsOnEdge(2,iEdge)) &amp;
+                                - pv_cell(k,cellsOnEdge(1,iEdge))) &amp;
+                               / dcEdge(iEdge)
+         enddo
+      enddo
+
+      !
+      ! Compute gradient of PV in the tangent direction
+      !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+      !
+      do iEdge = 1,nEdges
+         do k = 1,maxLevelEdgeBot(iEdge)
+           gradPVt(k,iEdge) = (  pv_vertex(k,verticesOnEdge(2,iEdge)) &amp;
+                               - pv_vertex(k,verticesOnEdge(1,iEdge))) &amp;
+                                 /dvEdge(iEdge)
+         enddo
+      enddo
+
+      !
+      ! Modify PV edge with upstream bias.
+      !
+      do iEdge = 1,nEdges
+         do k = 1,maxLevelEdgeBot(iEdge)
+           pv_edge(k,iEdge) = pv_edge(k,iEdge) &amp;
+             - 0.5 * dt* (  u(k,iEdge) * gradPVn(k,iEdge) &amp;
+                          + v(k,iEdge) * gradPVt(k,iEdge) )
+         enddo
+      enddo
+
+      !
+      ! equation of state
+      !
+      ! For an isopycnal model, density should remain constant.
+      ! For zlevel, calculate in-situ density
+      if (config_vert_grid_type.eq.'zlevel') then
+         call equation_of_state(s,grid,0,'relative')
+      ! mrp 110324 In order to visualize rhoDisplaced, include the following
+         call equation_of_state(s, grid, 1,'relative')
+      endif
+
+      !
+      ! Pressure
+      ! This section must be after computing rho
+      !
+      if (config_vert_grid_type.eq.'isopycnal') then
+
+        ! For Isopycnal model.
+        ! Compute pressure at top of each layer, and then
+        ! Montgomery Potential.
+        allocate(pTop(nVertLevels))
+        do iCell=1,nCells
+
+           ! assume atmospheric pressure at the surface is zero for now.
+           pTop(1) = 0.0
+           ! For isopycnal mode, p is the Montgomery Potential.
+           ! At top layer it is g*SSH, where SSH may be off by a 
+           ! constant (ie, h_s can be relative to top or bottom)
+           MontPot(1,iCell) = gravity &amp;
+              * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
+
+           do k=2,nVertLevels
+              pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
+
+              ! from delta M = p delta / rho
+              MontPot(k,iCell) = MontPot(k-1,iCell) &amp;
+                 + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell)) 
+           end do
+
+        end do
+        deallocate(pTop)
+
+      elseif (config_vert_grid_type.eq.'zlevel') then
+
+        ! For z-level model.
+        ! Compute pressure at middle of each level.  
+        ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
+        ! pressure at middle of layer including SSH.
+
+        do iCell=1,nCells
+           ! compute pressure for z-level coordinates
+           ! assume atmospheric pressure at the surface is zero for now.
+
+           pressure(1,iCell) = rho(1,iCell)*gravity &amp;
+              * (h(1,iCell)-0.5*hZLevel(1)) 
+
+           do k=2,maxLevelCell(iCell)
+              pressure(k,iCell) = pressure(k-1,iCell)  &amp;
+                + 0.5*gravity*(  rho(k-1,iCell)*hZLevel(k-1) &amp;
+                               + rho(k  ,iCell)*hZLevel(k  ))
+           end do
+
+        end do
+
+      endif
+
+      call OcnWtop(s,grid)
+
+      call timer_stop(&quot;OcnDiagnosticSolve&quot;)
+
+   end subroutine OcnDiagnosticSolve!}}}
+
+!***********************************************************************
+!
+!  routine OcnWtop
+!
+!&gt; \brief   Computes vertical velocity
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes the vertical velocity in the top layer for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnWtop(s, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Compute diagnostic fields used in the tendency computations
+   !
+   ! Input: grid - grid metadata
+   !
+   ! Output: s - computed diagnostics
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (state_type), intent(inout) :: s
+      type (mesh_type), intent(in) :: grid
+
+      ! mrp 110512 could clean this out, remove pointers?
+      integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+
+      integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
+        hZLevel
+      real (kind=RKIND), dimension(:,:), pointer :: u,wTop
+      real (kind=RKIND), dimension(:,:), allocatable:: div_u
+
+      integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &amp;
+        verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &amp;
+        boundaryEdge, boundaryCell
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &amp;
+        maxLevelVertexBot,  maxLevelVertexTop
+
+        call timer_start(&quot;wTop&quot;)
+
+      u           =&gt; s % u % array
+      wTop        =&gt; s % wTop % array
+
+      areaCell          =&gt; grid % areaCell % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      hZLevel           =&gt; grid % hZLevel % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeBot   =&gt; grid % maxLevelEdgeBot % array
+      dvEdge            =&gt; grid % dvEdge % array
+
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+      nVertLevels = grid % nVertLevels
+
+      !
+      ! vertical velocity through layer interface
+      !
+      if (config_vert_grid_type.eq.'isopycnal') then
+        ! set vertical velocity to zero in isopycnal case
+        wTop=0.0  
+
+      elseif (config_vert_grid_type.eq.'zlevel') then
+
+        !
+        ! Compute div(u) for each cell
+        ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+        !
+        allocate(div_u(nVertLevels,nCells+1))
+        div_u(:,:) = 0.0
+        do iEdge=1,nEdges
+           cell1 = cellsOnEdge(1,iEdge)
+           cell2 = cellsOnEdge(2,iEdge)
+           do k=2,maxLevelEdgeBot(iEdge)
+              flux = u(k,iEdge) * dvEdge(iEdge) 
+              div_u(k,cell1) = div_u(k,cell1) + flux
+              div_u(k,cell2) = div_u(k,cell2) - flux
+           end do 
+        end do 
+
+        do iCell=1,nCells
+           ! Vertical velocity through layer interface at top and 
+           ! bottom is zero.
+           wTop(1,iCell) = 0.0
+           wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+           do k=maxLevelCell(iCell),2,-1
+              wTop(k,iCell) = wTop(k+1,iCell) &amp;
+                 - div_u(k,iCell)/areaCell(iCell)*hZLevel(k)
+           end do
+        end do
+        deallocate(div_u)
+
+      endif
+
+      call timer_stop(&quot;wTop&quot;)
+
+   end subroutine OcnWtop!}}}
+
+!***********************************************************************
+!
+!  routine OcnFUPerp
+!
+!&gt; \brief   Computes f u_perp
+!&gt; \author  Doug Jacobsen
+!&gt; \date    23 September 2011
+!&gt; \version SVN:$Id$
+!&gt; \details 
+!&gt;  This routine computes f u_perp for the ocean
+!
+!-----------------------------------------------------------------------
+
+   subroutine OcnFUPerp(s, grid)!{{{
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Put f*uBcl^{perp} in u as a work variable
+   !
+   ! Input: s - current model state
+   !        grid - grid metadata
+   !
+   ! Output: tend - computed tendencies for prognostic variables
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (state_type), intent(inout) :: s
+      type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+!  Some of these variables can be removed, but at a later time.
+      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
+        vertex1, vertex2, eoe, i, j
+
+      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
+        upstream_bias, wTopEdge, rho0Inv, r
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
+        zMidZLevel, zTopZLevel 
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, uBcl, v, pressure, &amp;
+        tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
+        MontPot, wTop, divergence, vertViscTopOfEdge
+      type (dm_info) :: dminfo
+
+      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
+        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+      integer, dimension(:,:), pointer :: &amp;
+        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
+        edgesOnEdge, edgesOnVertex
+      real (kind=RKIND) :: u_diffusion
+      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+      real (kind=RKIND), dimension(:,:), pointer :: u_src
+      real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+      call timer_start(&quot;OcnFUPerp&quot;)
+
+      h           =&gt; s % h % array
+      u           =&gt; s % u % array
+      uBcl        =&gt; s % uBcl % array
+      v           =&gt; s % v % array
+      wTop        =&gt; s % wTop % array
+      h_edge      =&gt; s % h_edge % array
+      circulation =&gt; s % circulation % array
+      vorticity   =&gt; s % vorticity % array
+      divergence  =&gt; s % divergence % array
+      ke          =&gt; s % ke % array
+      ke_edge     =&gt; s % ke_edge % array
+      pv_edge     =&gt; s % pv_edge % array
+      MontPot     =&gt; s % MontPot % array
+      pressure    =&gt; s % pressure % array
+
+      weightsOnEdge     =&gt; grid % weightsOnEdge % array
+      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
+      cellsOnEdge       =&gt; grid % cellsOnEdge % array
+      cellsOnVertex     =&gt; grid % cellsOnVertex % array
+      verticesOnEdge    =&gt; grid % verticesOnEdge % array
+      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
+      edgesOnCell       =&gt; grid % edgesOnCell % array
+      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
+      edgesOnEdge       =&gt; grid % edgesOnEdge % array
+      edgesOnVertex     =&gt; grid % edgesOnVertex % array
+      dcEdge            =&gt; grid % dcEdge % array
+      dvEdge            =&gt; grid % dvEdge % array
+      areaCell          =&gt; grid % areaCell % array
+      areaTriangle      =&gt; grid % areaTriangle % array
+      h_s               =&gt; grid % h_s % array
+      fVertex           =&gt; grid % fVertex % array
+      fEdge             =&gt; grid % fEdge % array
+      zMidZLevel        =&gt; grid % zMidZLevel % array
+      zTopZLevel        =&gt; grid % zTopZLevel % array
+      maxLevelCell      =&gt; grid % maxLevelCell % array
+      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
+      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
+
+                  
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+      nEdgesSolve = grid % nEdgesSolve
+      nVertices   = grid % nVertices
+      nVertLevels = grid % nVertLevels
+
+      !
+      ! Put f*uBcl^{perp} in u as a work variable
+      !
+      do iEdge=1,grid % nEdgesSolve
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+
+         do k=1,maxLevelEdgeTop(iEdge)
+
+            u(k,iEdge) = 0.0
+            do j = 1,nEdgesOnEdge(iEdge)
+               eoe = edgesOnEdge(j,iEdge)
+               u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe) 
+            end do
+         end do
+      end do
+
+      call timer_stop(&quot;OcnFUPerp&quot;)
+
+   end subroutine OcnFUPerp!}}}
+
+!***********************************************************************
+
+end module OcnTendency
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker

Modified: branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F        2011-09-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/module_mpas_core.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -5,6 +5,8 @@
    use dmpar
    use test_cases
 
+   use OcnTendency
+
    use OcnVelPressureGrad
    use OcnVelVadv
    use OcnVelHmix
@@ -188,7 +190,7 @@
       integer :: i, iEdge, iCell, k
    
    
-      call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+      call OcnDiagnosticSolve(dt, block % state % time_levs(1) % state, mesh)
 
       call compute_mesh_scaling(mesh)
  

Modified: branches/ocean_projects/performance/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/module_time_integration.F        2011-09-23 18:09:34 UTC (rev 1018)
+++ branches/ocean_projects/performance/src/core_ocean/module_time_integration.F        2011-09-23 21:26:14 UTC (rev 1019)
@@ -8,20 +8,8 @@
    use spline_interpolation
    use timer
 
-   use OcnThickHadv
-   use OcnThickVadv
+   use OcnTendency
 
-   use OcnVelCoriolis
-   use OcnVelPressureGrad
-   use OcnVelVadv
-   use OcnVelHmix
-   use OcnVelForcing
-
-   use OcnTracerHadv
-   use OcnTracerVadv
-   use OcnTracerHmix
-   use OcnRestoring
-
    use OcnEquationOfState
    use OcnVmix
 
@@ -182,8 +170,8 @@
            if (.not.config_implicit_vertical_mix) then
               call OcnVmixCoefs(block % mesh, provis, block % diagnostics, err)
            end if
-           call compute_tend_h(block % tend, provis, block % diagnostics, block % mesh)
-           call compute_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+           call OcnTendH(block % tend, provis, block % diagnostics, block % mesh)
+           call OcnTendU(block % tend, provis, block % diagnostics, block % mesh)
 
            ! mrp 110718 filter btr mode out of u_tend
            ! still got h perturbations with just this alone.  Try to set uBtr=0 after full u computation
@@ -191,7 +179,7 @@
                call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
            endif
 
-           call compute_scalar_tend(block % tend, provis, block % diagnostics, block % mesh)
+           call OcnTendScalar(block % tend, provis, block % diagnostics, block % mesh)
            call enforce_boundaryEdge(block % tend, block % mesh)
            block =&gt; block % next
         end do
@@ -241,7 +229,7 @@
                  provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
               end if
 
-              call compute_solve_diagnostics(dt, provis, block % mesh)
+              call OcnDiagnosticSolve(dt, provis, block % mesh)
 
               block =&gt; block % next
            end do
@@ -318,41 +306,7 @@
             !  Implicit vertical solve for momentum
             !
             call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
-!           do iEdge=1,nEdges
-!             if (maxLevelEdgeTop(iEdge).gt.0) then
 
-               ! Compute A(k), C(k) for momentum
-               ! mrp 110315 efficiency note: for z-level, could precompute
-               ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
-               ! h_edge is computed in compute_solve_diag, and is not available yet.
-               ! This could be removed if hZLevel used instead.
-!              cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
-!              cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-!              do k=1,maxLevelEdgeTop(iEdge)
-!                 h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
-!              end do
-
-!              do k=1,maxLevelEdgeTop(iEdge)-1
-!                 A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &amp;
-!                    / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &amp;
-!                    / h_edge(k,iEdge)
-!              enddo
-!              A(maxLevelEdgeTop(iEdge)) = -dt*config_bottom_drag_coeff  &amp;
-!                 *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
-
-!              C(1) = 1 - A(1)
-!              do k=2,maxLevelEdgeTop(iEdge)
-!                 C(k) = 1 - A(k) - A(k-1)
-!              enddo
-
-!              call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
-
-!              u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
-!              u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
-
-!             end if
-!           end do
-
           !  mrp 110718 filter btr mode out of u
            if (config_rk_filter_btr_mode) then
                call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
@@ -364,32 +318,6 @@
             !
 
             call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
-!           do iCell=1,nCells
-               ! Compute A(k), C(k) for tracers
-               ! mrp 110315 efficiency note: for z-level, could precompute
-               ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
-!              do k=1,maxLevelCell(iCell)-1
-!                 A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &amp;
-!                    / (h(k,iCell) + h(k+1,iCell)) &amp;
-!                    / h(k,iCell)
-!              enddo
-
-!              A(maxLevelCell(iCell)) = 0.0
-
-!              C(1) = 1 - A(1)
-!              do k=2,maxLevelCell(iCell)
-!                 C(k) = 1 - A(k) - A(k-1)
-!              enddo
-
-!              call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &amp;
-!                 tracersTemp, maxLevelCell(iCell), &amp;
-!                 nVertLevels,num_tracers)
-
-!              tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
-!              tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
-!           end do
-!           deallocate(A,C,uTemp,tracersTemp)
-!           call timer_stop(&quot;RK4-implicit vert mix&quot;)
          end if
 
          ! mrp 110725 momentum decay term
@@ -419,7 +347,7 @@
             block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
          end if
 
-         call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+         call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
 
          call reconstruct(block % state % time_levs(2) % state, block % mesh)
 
@@ -557,7 +485,7 @@
          if (.not.config_implicit_vertical_mix) then
             call OcnVmixCoefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
          end if
-         call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+         call OcnTendU(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
          call enforce_boundaryEdge(block % tend, block % mesh)
          block =&gt; block % next
       end do
@@ -579,7 +507,7 @@
            allocate(uTemp(block % mesh % nVertLevels))
 
            ! Put f*uBcl^{perp} in uNew as a work variable
-           call compute_f_uperp(block % state % time_levs(2) % state , block % mesh)
+           call OcnFUPerp(block % state % time_levs(2) % state , block % mesh)
 
            do iEdge=1,block % mesh % nEdges
               cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
@@ -1189,7 +1117,7 @@
               block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
           endif
 
-         ! Put new sshEdge values in h_edge array, for the compute_scalar_tend call below.
+         ! Put new sshEdge values in h_edge array, for the OcnTendScalar call below.
              block % state % time_levs(2) % state % h_edge % array(1,iEdge) &amp;
            = sshEdge + block % mesh % hZLevel % array(1)
 
@@ -1200,7 +1128,7 @@
 
           end do ! iEdge
 
-         ! Put new SSH values in h array, for the compute_scalar_tend call below.
+         ! Put new SSH values in h array, for the OcnTendScalar call below.
          do iCell=1,block % mesh % nCells 
              block % state % time_levs(2) % state % h % array(1,iCell) &amp;
            = block % state % time_levs(2) % state % ssh % array(iCell) &amp;
@@ -1231,13 +1159,13 @@
          block =&gt; domain % blocklist
          do while (associated(block))
 
-           call compute_wTop(block % state % time_levs(2) % state, block % mesh)
+           call OcnWtop(block % state % time_levs(2) % state, block % mesh)
 
       if (trim(config_time_integration) == 'unsplit_explicit') then
-           call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+           call OcnTendH(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
       endif
 
-           call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+           call OcnTendScalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
 
            block =&gt; block % next
          end do
@@ -1328,7 +1256,7 @@
           ! compute u*, the velocity for tendency terms.  Put in uNew.
           ! uBclNew is at time n+1/2 here.
           ! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
-          ! The following must occur after  call compute_scalar_tend
+          ! The following must occur after  call OcnTendScalar
            do iEdge=1,block % mesh % nEdges
                block % state % time_levs(2) % state % u    % array(:,iEdge) &amp;
              = block % state % time_levs(2) % state % uBtr % array(iEdge) &amp;
@@ -1337,7 +1265,7 @@
 
          ! mrp 110512  I really only need this to compute h_edge, density, pressure.
          ! I can par this down later.
-         call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+         call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
  
 
          elseif (split_explicit_step == config_n_ts_iter) then
@@ -1430,7 +1358,7 @@
         endif
 
          do iCell=1,block % mesh % nCells
-         ! Put new SSH values in h array, for the compute_scalar_tend call below.
+         ! Put new SSH values in h array, for the OcnTendScalar call below.
              block % state % time_levs(2) % state % h % array(1,iCell) &amp;
            = block % state % time_levs(2) % state % ssh % array(iCell) &amp;
            + block % mesh % hZLevel % array(1)
@@ -1473,69 +1401,10 @@
 
             call OcnVelVmixTendImplicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
 
-!           do iEdge=1,block % mesh % nEdges
-!             if (maxLevelEdgeTop(iEdge).gt.0) then
-
-               ! Compute A(k), C(k) for momentum
-               ! mrp 110315 efficiency note: for z-level, could precompute
-               ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
-               ! h_edge is computed in compute_solve_diag, and is not available yet.
-               ! This could be removed if hZLevel used instead.
-!              cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
-!              cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-!              do k=1,maxLevelEdgeTop(iEdge)
-!                 h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
-!              end do
-
-!              do k=1,maxLevelEdgeTop(iEdge)-1
-!                 A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &amp;
-!                    / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &amp;
-!                    / h_edge(k,iEdge)
-!              enddo
-!              A(maxLevelEdgeTop(iEdge)) = 0.0
-
-!              C(1) = 1 - A(1)
-!              do k=2,maxLevelEdgeTop(iEdge)
-!                 C(k) = 1 - A(k) - A(k-1)
-!              enddo
-
-!              call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
-
-!              u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
-!              u(maxLevelEdgeTop(iEdge)+1:block % mesh % nVertLevels,iEdge) = 0.0
-
-!             end if
-!           end do
-
             !
             !  Implicit vertical solve for tracers
             !
             call OcnTracerVmixTendImplicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
-!           do iCell=1,block % mesh % nCells
-               ! Compute A(k), C(k) for tracers
-               ! mrp 110315 efficiency note: for z-level, could precompute
-               ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
-!              do k=1,maxLevelCell(iCell)-1
-!                 A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &amp;
-!                    / (h(k,iCell) + h(k+1,iCell)) &amp;
-!                    / h(k,iCell)
-!              enddo
-
-!              A(maxLevelCell(iCell)) = 0.0
-
-!              C(1) = 1 - A(1)
-!              do k=2,maxLevelCell(iCell)
-!                 C(k) = 1 - A(k) - A(k-1)
-!              enddo
-
-!              call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &amp;
-!                 tracersTemp, maxLevelCell(iCell), &amp;
-!                 block % mesh % nVertLevels,num_tracers)
-
-!              tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
-!              tracers(:,maxLevelCell(iCell)+1:block % mesh % nVertLevels,iCell) = -1e34
-!           end do
-!           deallocate(A,C,uTemp,tracersTemp)
          end if
 
          ! mrp 110725 adding momentum decay term
@@ -1562,7 +1431,7 @@
             block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
          end if
 
-         call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+         call OcnDiagnosticSolve(dt, block % state % time_levs(2) % state, block % mesh)
 
          call reconstruct(block % state % time_levs(2) % state, block % mesh)
 
@@ -1572,348 +1441,6 @@
 
    end subroutine split_explicit_timestep!}}}
 
-   subroutine compute_tend_h(tend, s, d, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute height and normal wind tendencies, as well as diagnostic variables
-   !
-   ! Input: s - current model state
-   !        grid - grid metadata
-   !
-   ! Output: tend - computed tendencies for prognostic variables
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (tend_type), intent(inout) :: tend
-      type (state_type), intent(in) :: s
-      type (diagnostics_type), intent(in) :: d
-      type (mesh_type), intent(in) :: grid
-
-      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
-        vertex1, vertex2, eoe, i, j, err
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-!  Most of these variables can be removed, but at a later time.
-      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
-      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
-        upstream_bias, wTopEdge, rho0Inv, r
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
-        zMidZLevel, zTopZLevel 
-      real (kind=RKIND), dimension(:,:), pointer :: &amp;
-        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &amp;
-        tend_h, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
-        MontPot, wTop, divergence, vertViscTopOfEdge
-      type (dm_info) :: dminfo
-
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
-      integer, dimension(:,:), pointer :: &amp;
-        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
-        edgesOnEdge, edgesOnVertex
-      real (kind=RKIND) :: u_diffusion
-      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-      call timer_start(&quot;compute_tend_h&quot;)
-
-      h           =&gt; s % h % array
-      u           =&gt; s % u % array
-      v           =&gt; s % v % array
-      wTop        =&gt; s % wTop % array
-      h_edge      =&gt; s % h_edge % array
-      circulation =&gt; s % circulation % array
-      vorticity   =&gt; s % vorticity % array
-      divergence  =&gt; s % divergence % array
-      ke          =&gt; s % ke % array
-      ke_edge     =&gt; s % ke_edge % array
-      pv_edge     =&gt; s % pv_edge % array
-      MontPot     =&gt; s % MontPot % array
-      pressure    =&gt; s % pressure % array
-      vertViscTopOfEdge =&gt; d % vertViscTopOfEdge % array
-
-      weightsOnEdge     =&gt; grid % weightsOnEdge % array
-      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      cellsOnVertex     =&gt; grid % cellsOnVertex % array
-      verticesOnEdge    =&gt; grid % verticesOnEdge % array
-      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
-      edgesOnCell       =&gt; grid % edgesOnCell % array
-      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
-      edgesOnEdge       =&gt; grid % edgesOnEdge % array
-      edgesOnVertex     =&gt; grid % edgesOnVertex % array
-      dcEdge            =&gt; grid % dcEdge % array
-      dvEdge            =&gt; grid % dvEdge % array
-      areaCell          =&gt; grid % areaCell % array
-      areaTriangle      =&gt; grid % areaTriangle % array
-      h_s               =&gt; grid % h_s % array
-      fVertex           =&gt; grid % fVertex % array
-      fEdge             =&gt; grid % fEdge % array
-      zMidZLevel        =&gt; grid % zMidZLevel % array
-      zTopZLevel        =&gt; grid % zTopZLevel % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
-      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
-
-      tend_h      =&gt; tend % h % array
-                  
-      nCells      = grid % nCells
-      nEdges      = grid % nEdges
-      nEdgesSolve = grid % nEdgesSolve
-      nVertices   = grid % nVertices
-      nVertLevels = grid % nVertLevels
-
-      !
-      ! height tendency: start accumulating tendency terms
-      !
-      tend_h = 0.0
-
-      !
-      ! height tendency: horizontal advection term -</font>
<font color="red">abla\cdot ( hu)
-      !
-      ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3. 
-      ! for explanation of divergence operator.
-      !
-      ! for z-level, only compute height tendency for top layer.
-
-      call timer_start(&quot;compute_tend_h-horiz adv&quot;)
-
-      call OcnThickHadvTend(grid, u, h_edge, tend_h, err)
-
-      call timer_stop(&quot;compute_tend_h-horiz adv&quot;)
-
-      !
-      ! height tendency: vertical advection term -d/dz(hw)
-      !
-      ! Vertical advection computed for top layer of a z grid only.
-      call timer_start(&quot;compute_tend_h-vert adv&quot;)
-
-      call OcnThickVadvTend(grid, wTop, tend_h, err)
-
-      call timer_stop(&quot;compute_tend_h-vert adv&quot;)
-      call timer_stop(&quot;compute_tend_h&quot;)
-   
-   end subroutine compute_tend_h!}}}
-
-   subroutine compute_tend_u(tend, s, d, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute height and normal wind tendencies, as well as diagnostic variables
-   !
-   ! Input: s - current model state
-   !        grid - grid metadata
-   !
-   ! Output: tend - computed tendencies for prognostic variables
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (tend_type), intent(inout) :: tend
-      type (state_type), intent(in) :: s
-      type (diagnostics_type), intent(in) :: d
-      type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-!  Some of these variables can be removed, but at a later time.
-      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
-        vertex1, vertex2, eoe, i, j
-
-      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, err
-      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
-        upstream_bias, wTopEdge, rho0Inv, r, visc_vorticity_coef
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
-        zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
-      real (kind=RKIND), dimension(:,:), pointer :: &amp;
-        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &amp;
-        tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
-        MontPot, wTop, divergence, vertViscTopOfEdge
-      type (dm_info) :: dminfo
-
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
-      integer, dimension(:,:), pointer :: &amp;
-        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
-        edgesOnEdge, edgesOnVertex
-      real (kind=RKIND) :: u_diffusion
-      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
-      real (kind=RKIND), dimension(:,:), pointer :: u_src
-      real (kind=RKIND), parameter :: rho_ref = 1000.0
-
-      call timer_start(&quot;compute_tend_u&quot;)
-
-      h           =&gt; s % h % array
-      u           =&gt; s % u % array
-      v           =&gt; s % v % array
-      wTop        =&gt; s % wTop % array
-      h_edge      =&gt; s % h_edge % array
-      circulation =&gt; s % circulation % array
-      vorticity   =&gt; s % vorticity % array
-      divergence  =&gt; s % divergence % array
-      ke          =&gt; s % ke % array
-      ke_edge     =&gt; s % ke_edge % array
-      pv_edge     =&gt; s % pv_edge % array
-      MontPot     =&gt; s % MontPot % array
-      pressure    =&gt; s % pressure % array
-      vertViscTopOfEdge =&gt; d % vertViscTopOfEdge % array
-
-      weightsOnEdge     =&gt; grid % weightsOnEdge % array
-      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      cellsOnVertex     =&gt; grid % cellsOnVertex % array
-      verticesOnEdge    =&gt; grid % verticesOnEdge % array
-      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
-      edgesOnCell       =&gt; grid % edgesOnCell % array
-      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
-      edgesOnEdge       =&gt; grid % edgesOnEdge % array
-      edgesOnVertex     =&gt; grid % edgesOnVertex % array
-      dcEdge            =&gt; grid % dcEdge % array
-      dvEdge            =&gt; grid % dvEdge % array
-      areaCell          =&gt; grid % areaCell % array
-      areaTriangle      =&gt; grid % areaTriangle % array
-      h_s               =&gt; grid % h_s % array
-! mrp 110516 cleanup fvertex fedge not used in this subroutine
-      fVertex           =&gt; grid % fVertex % array
-      fEdge             =&gt; grid % fEdge % array
-      zMidZLevel        =&gt; grid % zMidZLevel % array
-      zTopZLevel        =&gt; grid % zTopZLevel % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
-      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
-
-      tend_u      =&gt; tend % u % array
-                  
-      nCells      = grid % nCells
-      nEdges      = grid % nEdges
-      nEdgesSolve = grid % nEdgesSolve
-      nVertices   = grid % nVertices
-      nVertLevels = grid % nVertLevels
-
-      u_src =&gt; grid % u_src % array
-
-      meshScalingDel2 =&gt; grid % meshScalingDel2 % array
-      meshScalingDel4 =&gt; grid % meshScalingDel4 % array
-
-      !
-      ! velocity tendency: start accumulating tendency terms
-      !
-      ! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
-      tend_u(:,:) = 0.0
-
-      !
-      ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
-      !
-
-      call timer_start(&quot;compute_tend_u-coriolis&quot;)
-
-!     call OcnVelCoriolisTend(grid, pv_edge, h_edge, u, ke, tend_u, err)
-
-      do iEdge=1,grid % nEdgesSolve
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-
-         do k=1,maxLevelEdgeTop(iEdge)
-
-            q = 0.0
-            do j = 1,nEdgesOnEdge(iEdge)
-               eoe = edgesOnEdge(j,iEdge)
-               workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
-               q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe) 
-            end do
-
-           tend_u(k,iEdge) = tend_u(k,iEdge)     &amp;
-                  + q     &amp;
-                  - (   ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
-
-         end do
-      end do
-
-      call timer_stop(&quot;compute_tend_u-coriolis&quot;)
-
-      !
-      ! velocity tendency: vertical advection term -w du/dz
-      !
-      call timer_start(&quot;compute_tend_u-vert adv&quot;)
-
-      call OcnVelVadvTend(grid, u, wTop, tend_u, err)
-
-      call timer_stop(&quot;compute_tend_u-vert adv&quot;)
-
-      !
-      ! velocity tendency: pressure gradient
-      !
-      call timer_start(&quot;compute_tend_u-pressure grad&quot;)
-
-      if (config_vert_grid_type.eq.'isopycnal') then
-          call OcnVelPressureGradTend(grid, MontPot, tend_u, err)
-      elseif (config_vert_grid_type.eq.'zlevel') then
-          call OcnVelPressureGradTend(grid, pressure, tend_u, err)
-      end if
-
-      call timer_stop(&quot;compute_tend_u-pressure grad&quot;)
-
-      !
-      ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="red">abla^2 u
-      !   computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity )
-      !   strictly only valid for config_h_mom_eddy_visc2 == constant
-      !
-      call timer_start(&quot;compute_tend_u-horiz mix&quot;)
-
-      call OcnVelHmixTend(grid, divergence, vorticity, tend_u, err)
-
-      call timer_stop(&quot;compute_tend_u-horiz mix&quot;)
-
-      !
-      ! velocity tendency: forcing and bottom drag
-      !
-      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
-      ! know the bottom edge with nonzero velocity and place the drag there.
-
-      call timer_start(&quot;compute_tend_u-forcings&quot;)
-
-      call OcnVelForcingTend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
-
-      call timer_stop(&quot;compute_tend_u-forcings&quot;)
-
-      !
-      ! velocity tendency: vertical mixing d/dz( nu_v du/dz))
-      !
-      if (.not.config_implicit_vertical_mix) then
-          call timer_start(&quot;compute_tend_u-explicit vert mix&quot;)
-          call OcnVelVmixTendExplicit(grid, u, h_edge, vertViscTopOfEdge, tend_u, err)
-!        allocate(fluxVertTop(nVertLevels+1))
-!        fluxVertTop(1) = 0.0
-!        do iEdge=1,grid % nEdgesSolve
-!        
-!           do k=2,maxLevelEdgeTop(iEdge)
-!             fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &amp;
-!                * ( u(k-1,iEdge) - u(k,iEdge) ) &amp;
-!                * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
-!           enddo
-!           fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
-
-!           do k=1,maxLevelEdgeTop(iEdge)
-!             tend_u(k,iEdge) = tend_u(k,iEdge) &amp;
-!               + (fluxVertTop(k) - fluxVertTop(k+1)) &amp;
-!               / h_edge(k,iEdge)
-!           enddo
-
-!        end do
-!        deallocate(fluxVertTop)
-          call timer_stop(&quot;compute_tend_u-explicit vert mix&quot;)
-      endif
-      call timer_stop(&quot;compute_tend_u&quot;)
-
-   end subroutine compute_tend_u!}}}
-
    subroutine filter_btr_mode_tend_u(tend, s, d, grid)!{{{
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Filter and remove barotropic mode from the tendencies
@@ -2158,887 +1685,6 @@
 
    end subroutine filter_btr_mode_u!}}}
 
-   subroutine compute_f_uperp(s, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Put f*uBcl^{perp} in u as a work variable
-   !
-   ! Input: s - current model state
-   !        grid - grid metadata
-   !
-   ! Output: tend - computed tendencies for prognostic variables
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (state_type), intent(inout) :: s
-      type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-!  Some of these variables can be removed, but at a later time.
-      integer :: iEdge, iCell, iVertex, k, cell1, cell2, &amp;
-        vertex1, vertex2, eoe, i, j
-
-      integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
-      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &amp;
-        upstream_bias, wTopEdge, rho0Inv, r
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
-        zMidZLevel, zTopZLevel 
-      real (kind=RKIND), dimension(:,:), pointer :: &amp;
-        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, uBcl, v, pressure, &amp;
-        tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &amp;
-        MontPot, wTop, divergence, vertViscTopOfEdge
-      type (dm_info) :: dminfo
-
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
-      integer, dimension(:,:), pointer :: &amp;
-        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
-        edgesOnEdge, edgesOnVertex
-      real (kind=RKIND) :: u_diffusion
-      real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
-      real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
-      real (kind=RKIND), dimension(:,:), pointer :: u_src
-      real (kind=RKIND), parameter :: rho_ref = 1000.0
-
-      call timer_start(&quot;compute_f_uperp&quot;)
-
-      h           =&gt; s % h % array
-      u           =&gt; s % u % array
-      uBcl        =&gt; s % uBcl % array
-      v           =&gt; s % v % array
-      wTop        =&gt; s % wTop % array
-      h_edge      =&gt; s % h_edge % array
-      circulation =&gt; s % circulation % array
-      vorticity   =&gt; s % vorticity % array
-      divergence  =&gt; s % divergence % array
-      ke          =&gt; s % ke % array
-      ke_edge     =&gt; s % ke_edge % array
-      pv_edge     =&gt; s % pv_edge % array
-      MontPot     =&gt; s % MontPot % array
-      pressure    =&gt; s % pressure % array
-
-      weightsOnEdge     =&gt; grid % weightsOnEdge % array
-      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      cellsOnVertex     =&gt; grid % cellsOnVertex % array
-      verticesOnEdge    =&gt; grid % verticesOnEdge % array
-      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
-      edgesOnCell       =&gt; grid % edgesOnCell % array
-      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
-      edgesOnEdge       =&gt; grid % edgesOnEdge % array
-      edgesOnVertex     =&gt; grid % edgesOnVertex % array
-      dcEdge            =&gt; grid % dcEdge % array
-      dvEdge            =&gt; grid % dvEdge % array
-      areaCell          =&gt; grid % areaCell % array
-      areaTriangle      =&gt; grid % areaTriangle % array
-      h_s               =&gt; grid % h_s % array
-      fVertex           =&gt; grid % fVertex % array
-      fEdge             =&gt; grid % fEdge % array
-      zMidZLevel        =&gt; grid % zMidZLevel % array
-      zTopZLevel        =&gt; grid % zTopZLevel % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
-      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
-
-                  
-      nCells      = grid % nCells
-      nEdges      = grid % nEdges
-      nEdgesSolve = grid % nEdgesSolve
-      nVertices   = grid % nVertices
-      nVertLevels = grid % nVertLevels
-
-      !
-      ! Put f*uBcl^{perp} in u as a work variable
-      !
-      do iEdge=1,grid % nEdgesSolve
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-
-         do k=1,maxLevelEdgeTop(iEdge)
-
-            u(k,iEdge) = 0.0
-            do j = 1,nEdgesOnEdge(iEdge)
-               eoe = edgesOnEdge(j,iEdge)
-               u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe) 
-            end do
-         end do
-      end do
-
-      call timer_stop(&quot;compute_f_uperp&quot;)
-
-   end subroutine compute_f_uperp!}}}
-
-   subroutine compute_scalar_tend(tend, s, d, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   !
-   ! Input: s - current model state
-   !        grid - grid metadata
-   !        note: the variable s % tracers really contains the tracers, 
-   !              not tracers*h
-   !
-   ! Output: tend - computed scalar tendencies
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (tend_type), intent(inout) :: tend
-      type (state_type), intent(in) :: s
-      type (diagnostics_type), intent(in) :: d
-      type (mesh_type), intent(in) :: grid
-
-      integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&amp;
-        nEdges, nCells, nCellsSolve, nVertLevels, num_tracers, err
-      real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
-      real (kind=RKIND) :: flux, tracer_edge, r
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
-      real (kind=RKIND), dimension(:,:), pointer :: &amp;
-        u,h,wTop, h_edge, vertDiffTopOfCell
-      real (kind=RKIND), dimension(:,:,:), pointer :: &amp;
-        tracers, tend_tr
-      integer, dimension(:,:), pointer :: boundaryEdge
-      type (dm_info) :: dminfo
-
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
-      integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
-      real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &amp;
-         hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
-      real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &amp;
-            posZTopZLevel
-      real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
-      real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
-
-
-      real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
-      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
-      real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
-
-      integer :: index_temperature, index_salinity, rrr
-      real (kind=RKIND), dimension(:), pointer :: temperatureRestore, salinityRestore
-
-      call timer_start(&quot;compute_scalar_tend&quot;)
-
-      u           =&gt; s % u % array
-      h           =&gt; s % h % array
-      boundaryCell=&gt; grid % boundaryCell % array
-      wTop        =&gt; s % wTop % array
-      tracers     =&gt; s % tracers % array
-      h_edge      =&gt; s % h_edge % array
-      vertDiffTopOfCell =&gt; d % vertDiffTopOfCell % array
-
-      tend_tr     =&gt; tend % tracers % array
-                  
-      areaCell          =&gt; grid % areaCell % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      dvEdge            =&gt; grid % dvEdge % array
-      dcEdge            =&gt; grid % dcEdge % array
-      zTopZLevel        =&gt; grid % zTopZLevel % array
-      zMidZLevel        =&gt; grid % zMidZLevel % array
-      hRatioZLevelK    =&gt; grid % hRatioZLevelK % array
-      hRatioZLevelKm1    =&gt; grid % hRatioZLevelKm1 % array
-      boundaryEdge      =&gt; grid % boundaryEdge % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeTop      =&gt; grid % maxLevelEdgeTop % array
-      maxLevelVertexBot    =&gt; grid % maxLevelVertexBot % array
-
-      nEdges      = grid % nEdges
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nVertLevels = grid % nVertLevels
-      num_tracers = s % num_tracers
-
-      meshScalingDel2 =&gt; grid % meshScalingDel2 % array
-      meshScalingDel4 =&gt; grid % meshScalingDel4 % array
-
-
-      deriv_two   =&gt; grid % deriv_two % array
-
-      !
-      ! initialize tracer tendency (RHS of tracer equation) to zero.
-      !
-      tend_tr(:,:,:) = 0.0
-
-      !
-      ! tracer tendency: horizontal advection term -div( h \phi u)
-      !
-      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
-      ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
-      ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
-      ! tracer_edge at the boundary will also need to be defined for flux boundaries.
-
-      call timer_start(&quot;compute_scalar_tend-horiz adv&quot;)
-
-      call OcnTracerHadvTend(grid, u, h_edge, tracers, tend_tr, err)
-
-      call timer_stop(&quot;compute_scalar_tend-horiz adv&quot;)
-
-
-      !
-      ! tracer tendency: vertical advection term -d/dz( h \phi w)
-      !
-
-      call timer_start(&quot;compute_scalar_tend-vert adv&quot;)
-
-      call OcnTracerVadvTend(grid, wTop, tracers, tend_tr, err)
-
-      call timer_stop(&quot;compute_scalar_tend-vert adv&quot;)
-
-      !
-      ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
-      !
-      call timer_start(&quot;compute_scalar_tend-horiz diff&quot;)
-
-      call OcnTracerHmixTend(grid, h_edge, tracers, tend_tr, err)
-
-      call timer_stop(&quot;compute_scalar_tend-horiz diff&quot;)
-
-! mrp 110516 printing
-!print *, 'tend_tr 1',minval(tend_tr(3,1,1:nCells)),&amp;
-!                   maxval(tend_tr(3,1,1:nCells))
-!print *, 'tracer  1',minval(tracers(3,1,1:nCells)),&amp;
-!                   maxval(tracers(3,1,1:nCells))
-! mrp 110516 printing end
-
-      !
-      ! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
-      !
-      if (.not.config_implicit_vertical_mix) then
-         call timer_start(&quot;compute_scalar_tend-explicit vert diff&quot;)
-         call OcnTracerVmixTendExplicit(grid, h, vertDiffTopOfCell, tracers, tend_tr, err)
-!        allocate(fluxVertTop(num_tracers,nVertLevels+1))
-!        fluxVertTop(:,1) = 0.0
-!        do iCell=1,nCellsSolve 
-
-!           do k=2,maxLevelCell(iCell)
-!             do iTracer=1,num_tracers
-!               ! compute \kappa_v d\phi/dz
-!               fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &amp;
-!                  * (   tracers(iTracer,k-1,iCell)    &amp;
-!                      - tracers(iTracer,k  ,iCell) )  &amp;
-!                  * 2 / (h(k-1,iCell) + h(k,iCell))
-
-!             enddo
-!           enddo
-!           fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
-!           do k=1,maxLevelCell(iCell)
-!             do iTracer=1,num_tracers
-!               ! This is h d/dz( fluxVertTop) but h and dz cancel, so 
-!               ! reduces to delta( fluxVertTop)
-!               tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &amp;
-!                 + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
-!             enddo
-!           enddo
-!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
-!print '(a,50e12.2)', 'tend_tr    ',tend_tr(3,1,1:maxLevelCell(iCell))
-
-!        enddo ! iCell loop
-!        deallocate(fluxVertTop)
-         call timer_stop(&quot;compute_scalar_tend-explicit vert diff&quot;)
-      endif
-
-! mrp 110516 printing
-!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&amp;
-!                   maxval(tend_tr(3,1,1:nCells))
-! mrp 110516 printing end
-
-      !
-      ! add restoring to T and S in top model layer
-      !
-      call timer_start(&quot;compute_scalar_tend-restoring&quot;)
-
-      call OcnRestoringTend(grid, h, s%index_temperature, s%index_salinity, tracers, tend_tr, err)
-
-      call timer_stop(&quot;compute_scalar_tend-restoring&quot;)
-
- 10   format(2i8,10e20.10)
-      call timer_stop(&quot;compute_scalar_tend&quot;)
-
-   end subroutine compute_scalar_tend!}}}
-
-   subroutine compute_solve_diagnostics(dt, s, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute diagnostic fields used in the tendency computations
-   !
-   ! Input: grid - grid metadata
-   !
-   ! Output: s - computed diagnostics
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      real (kind=RKIND), intent(in) :: dt
-      type (state_type), intent(inout) :: s
-      type (mesh_type), intent(in) :: grid
-
-
-      integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
-      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
-
-      integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
-
-
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
-        hZLevel
-      real (kind=RKIND), dimension(:,:), pointer :: &amp;
-        weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&amp;
-        circulation, vorticity, ke, ke_edge, MontPot, wTop, &amp;
-        pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &amp;
-        rho, temperature, salinity
-      real (kind=RKIND), dimension(:,:,:), pointer :: tracers
-      real (kind=RKIND), dimension(:), allocatable:: pTop
-      real (kind=RKIND), dimension(:,:), allocatable:: div_u
-      character :: c1*6
-
-      integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &amp;
-        verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &amp;
-        boundaryEdge, boundaryCell
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &amp;
-        maxLevelVertexBot,  maxLevelVertexTop
-      real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
-      real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
-      real (kind=RKIND) :: coef_3rd_order
-      real (kind=RKIND) :: r, h1, h2
-
-      call timer_start(&quot;compute_solve_diagnostics&quot;)
-
-      h           =&gt; s % h % array
-      u           =&gt; s % u % array
-      v           =&gt; s % v % array
-      wTop        =&gt; s % wTop % array
-      h_edge      =&gt; s % h_edge % array
-      circulation =&gt; s % circulation % array
-      vorticity   =&gt; s % vorticity % array
-      divergence  =&gt; s % divergence % array
-      ke          =&gt; s % ke % array
-      ke_edge     =&gt; s % ke_edge % array
-      pv_edge     =&gt; s % pv_edge % array
-      pv_vertex   =&gt; s % pv_vertex % array
-      pv_cell     =&gt; s % pv_cell % array
-      gradPVn     =&gt; s % gradPVn % array
-      gradPVt     =&gt; s % gradPVt % array
-      rho         =&gt; s % rho % array
-      tracers     =&gt; s % tracers % array
-      MontPot     =&gt; s % MontPot % array
-      pressure    =&gt; s % pressure % array
-
-      weightsOnEdge     =&gt; grid % weightsOnEdge % array
-      kiteAreasOnVertex =&gt; grid % kiteAreasOnVertex % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      cellsOnVertex     =&gt; grid % cellsOnVertex % array
-      verticesOnEdge    =&gt; grid % verticesOnEdge % array
-      nEdgesOnCell      =&gt; grid % nEdgesOnCell % array
-      edgesOnCell       =&gt; grid % edgesOnCell % array
-      nEdgesOnEdge      =&gt; grid % nEdgesOnEdge % array
-      edgesOnEdge       =&gt; grid % edgesOnEdge % array
-      edgesOnVertex     =&gt; grid % edgesOnVertex % array
-      dcEdge            =&gt; grid % dcEdge % array
-      dvEdge            =&gt; grid % dvEdge % array
-      areaCell          =&gt; grid % areaCell % array
-      areaTriangle      =&gt; grid % areaTriangle % array
-      h_s               =&gt; grid % h_s % array
-      fVertex           =&gt; grid % fVertex % array
-      fEdge             =&gt; grid % fEdge % array
-      hZLevel           =&gt; grid % hZLevel % array
-      deriv_two         =&gt; grid % deriv_two % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeTop   =&gt; grid % maxLevelEdgeTop % array
-      maxLevelEdgeBot   =&gt; grid % maxLevelEdgeBot % array
-      maxLevelVertexBot =&gt; grid % maxLevelVertexBot % array
-      maxLevelVertexTop =&gt; grid % maxLevelVertexTop % array
-                  
-      nCells      = grid % nCells
-      nEdges      = grid % nEdges
-      nVertices   = grid % nVertices
-      nVertLevels = grid % nVertLevels
-      vertexDegree = grid % vertexDegree
-
-      boundaryEdge =&gt; grid % boundaryEdge % array
-      boundaryCell =&gt; grid % boundaryCell % array
-
-      !
-      ! Compute height on cell edges at velocity locations
-      !   Namelist options control the order of accuracy of the reconstructed h_edge value
-      !
-      ! mrp 101115 note: in order to include flux boundary conditions, we will need to 
-      ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
-
-      ! mrp 110516 efficiency note: For z-level, only do this on level 1.  h_edge for all
-      ! lower levels is defined by hZlevel.
-
-      call timer_start(&quot;compute_solve_diagnostics-hEdge&quot;)
-
-      coef_3rd_order = 0.
-      if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
-      if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
-      if (config_thickness_adv_order == 2) then
-          call timer_start(&quot;compute_solve_diagnostics-hEdge 2&quot;)
-
-         do iEdge=1,nEdges
-            cell1 = cellsOnEdge(1,iEdge)
-            cell2 = cellsOnEdge(2,iEdge)
-            do k=1,maxLevelEdgeTop(iEdge)
-               h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
-            end do
-         end do
-          call timer_stop(&quot;compute_solve_diagnostics-hEdge 2&quot;)
-
-      else if (config_thickness_adv_order == 3) then
-          call timer_start(&quot;compute_solve_diagnostics-hEdge 3&quot;)
-
-         do iEdge=1,nEdges
-            cell1 = cellsOnEdge(1,iEdge)
-            cell2 = cellsOnEdge(2,iEdge)
-
-            do k=1,maxLevelEdgeTop(iEdge)
-
-               d2fdx2_cell1 = 0.0
-               d2fdx2_cell2 = 0.0
-
-               !-- if not a boundary cell
-               if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
-                  d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
-                  d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
-                  !-- all edges of cell 1
-                  do i=1, grid % nEdgesOnCell % array (cell1)
-                          d2fdx2_cell1 = d2fdx2_cell1 + &amp;
-                          deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
-                  end do
-
-                  !-- all edges of cell 2
-                  do i=1, grid % nEdgesOnCell % array (cell2)
-                          d2fdx2_cell2 = d2fdx2_cell2 + &amp;
-                          deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
-                  end do
-
-               endif
-
-               !-- if u &gt; 0:
-               if (u(k,iEdge) &gt; 0) then
-                  h_edge(k,iEdge) =     &amp;
-                       0.5*(h(k,cell1) + h(k,cell2))      &amp;
-                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &amp;
-                       -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
-               !-- else u &lt;= 0:
-               else
-                  h_edge(k,iEdge) =     &amp;
-                       0.5*(h(k,cell1) + h(k,cell2))      &amp;
-                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.          &amp;
-                       +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
-               end if
-
-            end do   ! do k
-         end do         ! do iEdge
-
-          call timer_stop(&quot;compute_solve_diagnostics-hEdge 3&quot;)
-      else  if (config_thickness_adv_order == 4) then
-          call timer_start(&quot;compute_solve_diagnostics-hEdge 4&quot;)
-
-         do iEdge=1,nEdges
-            cell1 = cellsOnEdge(1,iEdge)
-            cell2 = cellsOnEdge(2,iEdge)
-
-            do k=1,maxLevelEdgeTop(iEdge)
-
-               d2fdx2_cell1 = 0.0
-               d2fdx2_cell2 = 0.0
-
-               !-- if not a boundary cell
-               if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
-                  d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
-                  d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
-                  !-- all edges of cell 1
-                  do i=1, grid % nEdgesOnCell % array (cell1)
-                          d2fdx2_cell1 = d2fdx2_cell1 + &amp;
-                          deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
-                  end do
-
-                  !-- all edges of cell 2
-                  do i=1, grid % nEdgesOnCell % array (cell2)
-                          d2fdx2_cell2 = d2fdx2_cell2 + &amp;
-                          deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
-                  end do
-
-               endif
-
-               h_edge(k,iEdge) =   &amp;
-                    0.5*(h(k,cell1) + h(k,cell2))      &amp;
-                       -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
-            end do   ! do k
-         end do         ! do iEdge
-
-         call timer_stop(&quot;compute_solve_diagnostics-hEdge 4&quot;)
-      endif   ! if(config_thickness_adv_order == 2)
-      call timer_stop(&quot;compute_solve_diagnostics-hEdge&quot;)
-
-      !
-      ! set the velocity and height at dummy address
-      !    used -1e34 so error clearly occurs if these values are used.
-      !
-!mrp 110516 change to zero, change back later:
-      u(:,nEdges+1) = -1e34
-      h(:,nCells+1) = -1e34
-      tracers(s % index_temperature,:,nCells+1) = -1e34
-      tracers(s % index_salinity,:,nCells+1) = -1e34
-
-      !
-      ! Compute circulation and relative vorticity at each vertex
-      !
-      circulation(:,:) = 0.0
-      do iEdge=1,nEdges
-         vertex1 = verticesOnEdge(1,iEdge)
-         vertex2 = verticesOnEdge(2,iEdge)
-         do k=1,maxLevelEdgeBot(iEdge)
-            circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
-            circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
-         end do
-      end do
-      do iVertex=1,nVertices
-         do k=1,maxLevelVertexBot(iVertex)
-            vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
-         end do
-      end do
-
-      !
-      ! Compute the divergence at each cell center
-      !
-      divergence(:,:) = 0.0
-      do iEdge=1,nEdges
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-         do k=1,maxLevelEdgeBot(iEdge)
-             divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
-             divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
-         enddo
-      end do
-      do iCell = 1,nCells
-         r = 1.0 / areaCell(iCell)
-         do k = 1,maxLevelCell(iCell)
-            divergence(k,iCell) = divergence(k,iCell) * r
-         enddo
-      enddo
-
-      !
-      ! Compute kinetic energy in each cell
-      !
-      ke(:,:) = 0.0
-      do iEdge=1,nEdges
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-         do k=1,maxLevelEdgeBot(iEdge)
-              ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
-              ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
-         enddo
-      end do
-      do iCell = 1,nCells
-         do k = 1,maxLevelCell(iCell)
-            ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
-         enddo
-      enddo
-
-      !
-      ! Compute v (tangential) velocities
-      !
-      v(:,:) = 0.0
-      do iEdge = 1,nEdges
-         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,maxLevelEdgeTop(iEdge) 
-               v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
-            end do
-         end do
-      end do
-
-      !
-      ! Compute ke on cell edges at velocity locations for quadratic bottom drag. 
-      !
-      ! mrp 101025 efficiency note: we could get rid of ke_edge completely by 
-      ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
-      ke_edge = 0.0  !mrp remove 0 for efficiency
-      do iEdge=1,nEdges
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-         do k=1,maxLevelEdgeTop(iEdge)
-            ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
-         end do
-      end do
-
-      !
-      ! Compute height at vertices, pv at vertices, and average pv to edge locations
-      !  ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
-      !
-      if (trim(config_time_integration) == 'RK4') then
-         ! for RK4, PV is really PV = (eta+f)/h
-         fCoef = 1
-      elseif (trim(config_time_integration) == 'split_explicit' &amp;
-          .or.trim(config_time_integration) == 'unsplit_explicit') then
-         ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
-! mrp temp, new should be:
-         fCoef = 0
-! old, for testing:
-!         fCoef = 1
-      end if
-
-      do iVertex = 1,nVertices
-         do k=1,maxLevelVertexBot(iVertex)
-            h_vertex = 0.0
-            do i=1,vertexDegree
-               h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
-            end do
-            h_vertex = h_vertex / areaTriangle(iVertex)
-
-            pv_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
-         end do
-      end do
-
-      !
-      ! Compute pv at cell centers
-      !    ( this computes pv_cell for all real cells and distance-1 ghost cells )
-      !
-      pv_cell(:,:) = 0.0
-      do iVertex = 1,nVertices
-         do i=1,vertexDegree
-            iCell = cellsOnVertex(i,iVertex)
-            do k = 1,maxLevelCell(iCell)
-               pv_cell(k,iCell) = pv_cell(k,iCell)  &amp;
-                  + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &amp;
-                    / areaCell(iCell)
-            enddo
-         enddo
-      enddo
-
-      !
-      ! Compute pv at the edges
-      !   ( this computes pv_edge at all edges bounding real cells )
-      !
-      pv_edge(:,:) = 0.0
-      do iVertex = 1,nVertices
-         do i=1,vertexDegree
-            iEdge = edgesOnVertex(i,iVertex)
-            do k=1,maxLevelEdgeBot(iEdge)
-               pv_edge(k,iEdge) =  pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
-            enddo
-        end do
-      end do
-
-      !
-      ! Compute gradient of PV in normal direction
-      !   ( this computes gradPVn for all edges bounding real cells )
-      !
-      gradPVn(:,:) = 0.0
-      do iEdge = 1,nEdges
-         do k=1,maxLevelEdgeTop(iEdge)
-            gradPVn(k,iEdge) = (  pv_cell(k,cellsOnEdge(2,iEdge)) &amp;
-                                - pv_cell(k,cellsOnEdge(1,iEdge))) &amp;
-                               / dcEdge(iEdge)
-         enddo
-      enddo
-
-      !
-      ! Compute gradient of PV in the tangent direction
-      !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
-      !
-      do iEdge = 1,nEdges
-         do k = 1,maxLevelEdgeBot(iEdge)
-           gradPVt(k,iEdge) = (  pv_vertex(k,verticesOnEdge(2,iEdge)) &amp;
-                               - pv_vertex(k,verticesOnEdge(1,iEdge))) &amp;
-                                 /dvEdge(iEdge)
-         enddo
-      enddo
-
-      !
-      ! Modify PV edge with upstream bias.
-      !
-      do iEdge = 1,nEdges
-         do k = 1,maxLevelEdgeBot(iEdge)
-           pv_edge(k,iEdge) = pv_edge(k,iEdge) &amp;
-             - 0.5 * dt* (  u(k,iEdge) * gradPVn(k,iEdge) &amp;
-                          + v(k,iEdge) * gradPVt(k,iEdge) )
-         enddo
-      enddo
-
-      !
-      ! equation of state
-      !
-      ! For an isopycnal model, density should remain constant.
-      ! For zlevel, calculate in-situ density
-      if (config_vert_grid_type.eq.'zlevel') then
-         call equation_of_state(s,grid,0,'relative')
-      ! mrp 110324 In order to visualize rhoDisplaced, include the following
-         call equation_of_state(s, grid, 1,'relative')
-      endif
-
-      !
-      ! Pressure
-      ! This section must be after computing rho
-      !
-      if (config_vert_grid_type.eq.'isopycnal') then
-
-        ! For Isopycnal model.
-        ! Compute pressure at top of each layer, and then
-        ! Montgomery Potential.
-        allocate(pTop(nVertLevels))
-        do iCell=1,nCells
-
-           ! assume atmospheric pressure at the surface is zero for now.
-           pTop(1) = 0.0
-           ! For isopycnal mode, p is the Montgomery Potential.
-           ! At top layer it is g*SSH, where SSH may be off by a 
-           ! constant (ie, h_s can be relative to top or bottom)
-           MontPot(1,iCell) = gravity &amp;
-              * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
-
-           do k=2,nVertLevels
-              pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
-
-              ! from delta M = p delta / rho
-              MontPot(k,iCell) = MontPot(k-1,iCell) &amp;
-                 + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell)) 
-           end do
-
-        end do
-        deallocate(pTop)
-
-      elseif (config_vert_grid_type.eq.'zlevel') then
-
-        ! For z-level model.
-        ! Compute pressure at middle of each level.  
-        ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
-        ! pressure at middle of layer including SSH.
-
-        do iCell=1,nCells
-           ! compute pressure for z-level coordinates
-           ! assume atmospheric pressure at the surface is zero for now.
-
-           pressure(1,iCell) = rho(1,iCell)*gravity &amp;
-              * (h(1,iCell)-0.5*hZLevel(1)) 
-
-           do k=2,maxLevelCell(iCell)
-              pressure(k,iCell) = pressure(k-1,iCell)  &amp;
-                + 0.5*gravity*(  rho(k-1,iCell)*hZLevel(k-1) &amp;
-                               + rho(k  ,iCell)*hZLevel(k  ))
-           end do
-
-        end do
-
-      endif
-
-      call compute_wtop(s,grid)
-
-      call timer_stop(&quot;compute_solve_diagnostics&quot;)
-
-   end subroutine compute_solve_diagnostics!}}}
-
-   subroutine compute_wtop(s, grid)!{{{
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute diagnostic fields used in the tendency computations
-   !
-   ! Input: grid - grid metadata
-   !
-   ! Output: s - computed diagnostics
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (state_type), intent(inout) :: s
-      type (mesh_type), intent(in) :: grid
-
-      ! mrp 110512 could clean this out, remove pointers?
-      integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
-      real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
-
-      integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
-
-
-      real (kind=RKIND), dimension(:), pointer :: &amp;
-        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &amp;
-        hZLevel
-      real (kind=RKIND), dimension(:,:), pointer :: u,wTop
-      real (kind=RKIND), dimension(:,:), allocatable:: div_u
-
-      integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &amp;
-        verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &amp;
-        boundaryEdge, boundaryCell
-      integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &amp;
-        maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &amp;
-        maxLevelVertexBot,  maxLevelVertexTop
-
-        call timer_start(&quot;wTop&quot;)
-
-      u           =&gt; s % u % array
-      wTop        =&gt; s % wTop % array
-
-      areaCell          =&gt; grid % areaCell % array
-      cellsOnEdge       =&gt; grid % cellsOnEdge % array
-      hZLevel           =&gt; grid % hZLevel % array
-      maxLevelCell      =&gt; grid % maxLevelCell % array
-      maxLevelEdgeBot   =&gt; grid % maxLevelEdgeBot % array
-      dvEdge            =&gt; grid % dvEdge % array
-
-      nCells      = grid % nCells
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-
-      !
-      ! vertical velocity through layer interface
-      !
-      if (config_vert_grid_type.eq.'isopycnal') then
-        ! set vertical velocity to zero in isopycnal case
-        wTop=0.0  
-
-      elseif (config_vert_grid_type.eq.'zlevel') then
-
-        !
-        ! Compute div(u) for each cell
-        ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
-        !
-        allocate(div_u(nVertLevels,nCells+1))
-        div_u(:,:) = 0.0
-        do iEdge=1,nEdges
-           cell1 = cellsOnEdge(1,iEdge)
-           cell2 = cellsOnEdge(2,iEdge)
-           do k=2,maxLevelEdgeBot(iEdge)
-              flux = u(k,iEdge) * dvEdge(iEdge) 
-              div_u(k,cell1) = div_u(k,cell1) + flux
-              div_u(k,cell2) = div_u(k,cell2) - flux
-           end do 
-        end do 
-
-        do iCell=1,nCells
-           ! Vertical velocity through layer interface at top and 
-           ! bottom is zero.
-           wTop(1,iCell) = 0.0
-           wTop(maxLevelCell(iCell)+1,iCell) = 0.0
-           do k=maxLevelCell(iCell),2,-1
-              wTop(k,iCell) = wTop(k+1,iCell) &amp;
-                 - div_u(k,iCell)/areaCell(iCell)*hZLevel(k)
-           end do
-        end do
-        deallocate(div_u)
-
-      endif
-
-      call timer_stop(&quot;wTop&quot;)
-
-   end subroutine compute_wtop!}}}
-
    subroutine enforce_boundaryEdge(tend, grid)!{{{
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Enforce any boundary conditions on the normal velocity at each edge

</font>
</pre>