<p><b>mpetersen@lanl.gov</b> 2010-04-09 11:14:34 -0600 (Fri, 09 Apr 2010)</p><p>Add computation of horizontal thickness flux to compute_tend.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry        2010-04-08 19:38:18 UTC (rev 186)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry        2010-04-09 17:14:34 UTC (rev 187)
@@ -112,6 +112,7 @@
 var real    pmid ( nVertLevels nCells Time ) o pmid - -
 var real    pbot ( nVertLevels nCells Time ) o pbot - -
 var real    MontPot ( nVertLevels nCells Time ) o MontPot - -
+var real    w ( nVertLevels nCells Time ) o w - -
 # mrp 100112 end
 
 

Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F        2010-04-08 19:38:18 UTC (rev 186)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F        2010-04-09 17:14:34 UTC (rev 187)
@@ -256,17 +256,18 @@
       real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
 
       integer :: nCells, nEdges, nVertices, nVertLevels
-      real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
-      real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &amp;
-                                                    circulation, vorticity, ke, pv_edge, divergence
-      integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+      real (kind=RKIND), dimension(:), pointer :: &amp;
+        h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+      real (kind=RKIND), dimension(:,:), pointer :: &amp;
+        vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, &amp;
+        tend_h, tend_u, Fv, circulation, vorticity, ke, pv_edge, &amp;
+        divergence, MontPot
       integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+      integer, dimension(:,:), pointer :: &amp;
+        cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &amp;
+        edgesOnEdge, edgesOnVertex
       real (kind=RKIND) :: u_diffusion, visc
 
-      !mrp 100112:
-      real (kind=RKIND), dimension(:,:), pointer :: MontPot
-      !mrp 100112 end
-
 !ocean
       real (kind=RKIND), dimension(:,:), pointer :: u_src
       real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -308,6 +309,7 @@
 
       tend_h      =&gt; tend % h % array
       tend_u      =&gt; tend % u % array
+      Fv          =&gt; tend % w % array
                   
       nCells      = grid % nCells
       nEdges      = grid % nEdges
@@ -321,6 +323,7 @@
 
       !
       ! Compute height tendency for each cell
+      ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
       !
       tend_h(:,:) = 0.0
       do iEdge=1,nEdges
@@ -339,11 +342,51 @@
                end do 
             end if
       end do 
+
+      ! mrp100409 computation of h tendency was, only had div.(hu) term
+      !do iCell=1,grid % nCellsSolve
+      !   do k=1,nVertLevels
+      !      tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+      !   end do
+      !end do
+      ! mrp100409 replace with
       do iCell=1,grid % nCellsSolve
          do k=1,nVertLevels
+            ! tend_h is just the -div.(hu) term at this point:
             tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
          end do
+
+         k=nVertLevels
+            Fv(k,iCell) = tend_h(k,iCell)*h(k,icell)
+   
+            ! note if you substitute line above this is just tend_h=0.
+            ! so this line can be changed to tend_h=0 in the future.
+            tend_h(k,iCell) =   tend_h(k,iCell) &amp;
+                              - (Fv(k,iCell))/h(k,icell)
+
+         do k=nVertLevels-1,-1,2
+
+            ! F^v_{k-1/2} = 
+            ! F^v_{k+1/2} - </font>
<font color="blue">abla_h \cdot \left( F^h \right) h_k
+            ! note +tend_h because tend_h is -div.(hu)
+            Fv(k,iCell) = Fv(k+1,iCell) + tend_h(k,iCell)*h(k,icell)
+    
+            ! now tend_h is div.(hu) + d/dz(hw):
+            ! - </font>
<font color="gray">abla_h \cdot F^h_k 
+            ! - \frac{F^v_{k-1/2} - F^v_{k+1/2}}{ h^n_k} 
+            ! note if you substitute line above this is just tend_h=0.
+            ! so this line can be changed to tend_h=0 in the future.
+            tend_h(k,iCell) =   tend_h(k,iCell) &amp;
+                              - (Fv(k,iCell) - Fv(k+1,iCell))/h(k,icell)
+         end do
+
+         k=1
+            Fv(k,iCell) = 0.0
+            tend_h(k,iCell) =   tend_h(k,iCell) &amp;
+                              - ( - Fv(k+1,iCell))/h(k,icell)
+
       end do
+      ! mrp100409 end
 
 #ifdef LANL_FORMULATION
       !
@@ -454,10 +497,10 @@
                      flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) * tracer_edge
                      tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
                      tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
-                  end do 
-               end do 
+                  end do
+               end do
             end if
-      end do 
+      end do
 
       do iCell=1,grid % nCellsSolve
          do k=1,grid % nVertLevelsSolve

</font>
</pre>