<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, &
- circulation, vorticity, ke, pv_edge, divergence
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+ 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, Fv, circulation, vorticity, ke, pv_edge, &
+ divergence, MontPot
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ 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 => tend % h % array
tend_u => tend % u % array
+ Fv => 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) &
+ - (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) &
+ - (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) &
+ - ( - 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>