<p><b>mpetersen@lanl.gov</b> 2010-05-12 14:15:56 -0600 (Wed, 12 May 2010)</p><p>Compute wBot in diagnostics. Found bug: Normalization tend_h(k,iCell) / areaCell(iCell) should loop to nCells, not nCellsSolve, because I use that to compute advective terms before a boundary update. Also updated some variable names.<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-05-11 21:20:41 UTC (rev 265)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-12 20:15:56 UTC (rev 266)
@@ -132,6 +132,7 @@
var real pMidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
var real w ( nVertLevels nCells Time ) o w - -
+var real wBot ( nVertLevels nCells Time ) o wBot - -
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
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-05-11 21:20:41 UTC (rev 265)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-12 20:15:56 UTC (rev 266)
@@ -284,12 +284,13 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wHat, rho0Inv
+ upstream_bias, wBotEdge, rho0Inv
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zSurfaceEdge
- real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudz
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zSurfaceEdge, &
+ zMidZLevel
+ real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudzBotEdge
real (kind=RKIND), dimension(:,:), pointer :: &
- vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, &
+ vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wBot, &
tend_h, tend_u, Fv, circulation, vorticity, ke, pv_edge, &
divergence, MontPot, pmidZLevel, zMidEdge, zBotEdge
@@ -308,6 +309,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
+ wBot => s % wBot % array
Fv => tend % w % array
h_edge => s % h_edge % array
circulation => s % circulation % array
@@ -339,6 +341,7 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
+ zmidZLevel => grid % zmidZLevel % array
tend_h => tend % h % array
tend_u => tend % u % array
@@ -350,7 +353,6 @@
u_src => grid % u_src % array
-
!
! Compute height tendency for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
@@ -373,29 +375,17 @@
end if
end do
- ! mrp 100409 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
- ! mrp 100409 replace with
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
if (config_vert_grid_type.eq.'isopycnal') then
Fv=0.0 ! vertical fluxes are zero
- do iCell=1,grid % nCellsSolve
- do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
- end do
- end do
elseif (config_vert_grid_type.eq.'zlevel') then
- ! z-level with vertical advection
- 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
+ do iCell=1,nCells
! change nvertlevels to kmt later
! this next line can be set permanently somewhere upon startup.
@@ -448,33 +438,34 @@
vertex2 = verticesOnEdge(2,iEdge)
! mrp 100422 add -w*dudz term to tendancy
- ! For isopycnal mode, Fv should be zero.
+ ! For isopycnal mode, wBot should be zero.
! change nvertlevels to kmt later
do k=1,nVertLevels-1
- ! w =Fv/h. Also average w from cell center to edge
- wHat = ( Fv(k,cell1)/h(k+1,cell1) &
- + Fv(k,cell2)/h(k+1,cell2))/2.0
+ ! Average w from cell center to edge
+ wBotEdge = 0.5*(wBot(k,cell1)+wBot(k,cell2))
+
! compute dudz at vertical interface with first order derivative.
- w_dudz(k) = wHat * (u(k,iEdge)-u(k+1,iEdge)) &
+! mrp 100511 new:
+! w_dudzBotEdge(k) = wBotEdge * (u(k,iEdge)-u(k+1,iEdge)) &
+! / (zMidZLevel(k) - zMidZLevel(k+1))
+! mrp 100511 old:
+ w_dudzBotEdge(k) = wBotEdge * (u(k,iEdge)-u(k+1,iEdge)) &
* 2.0 / (h_edge(k,iEdge) + h_edge(k+1,iEdge))
end do
- w_dudz(nVertLevels) = 0.0
+ w_dudzBotEdge(nVertLevels) = 0.0
- ! Average w*du/dz from vertical edges to vertical middle of cell
-!mrp 100510 remove zlevel additions
-! tend_u(1,iEdge) = - 0.5 * w_dudz(1)
+ ! Average w*du/dz from vertical interface to vertical middle of cell
+ tend_u(1,iEdge) = - 0.5 * w_dudzBotEdge(1)
do k=2,nVertLevels
-!mrp 100510 remove zlevel additions
-! tend_u(k,iEdge) = - 0.5 * (w_dudz(k-1) + w_dudz(k))
+ tend_u(k,iEdge) = - 0.5 * (w_dudzBotEdge(k-1) + w_dudzBotEdge(k))
enddo
! mrp 100422 end: add -w*dudz term to tendancy
! mrp 100426: add pressure gradient
if (config_vert_grid_type.eq.'isopycnal') then
do k=1,nVertLevels
-!mrp 100510 remove zlevel additions
-! tend_u(k,iEdge) = tend_u(k,iEdge) &
-! - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
elseif (config_vert_grid_type.eq.'zlevel') then
do k=1,nVertLevels
@@ -501,19 +492,10 @@
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
-!mrp 100510 remove zlevel additions
- ! tend_u(k,iEdge) = tend_u(k,iEdge) &
- ! + q &
- ! + u_diffusion &
- ! - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
-!mrp 100510 remove zlevel additions replace with
- tend_u(k,iEdge) = &
- q &
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + q &
+ u_diffusion &
- - ( ke(k,cell2) - ke(k,cell1) &
- + MontPot(k,cell2) - MontPot(k,cell1) &
- ) / dcEdge(iEdge)
-!mrp 100510 remove zlevel additions end
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
end do
end do
@@ -548,7 +530,6 @@
dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
endif
-!mrp 100510 remove zlevel additions
tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
enddo
@@ -903,11 +884,12 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zSurface, hZLevel, zSurfaceEdge
real (kind=RKIND), dimension(:,:), pointer :: &
- vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, &
- tend_h, tend_u, circulation, vorticity, ke, &
+ vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wBot, &
+ circulation, vorticity, ke, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zmid, zbot, zmidEdge, zbotEdge, pmid, pbot, MontPot, rho, temperature, salinity, pmidZLevel
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
@@ -919,10 +901,9 @@
u => s % u % array
v => s % v % array
w => s % w % array
+ wBot => s % wBot % array
vh => s % vh % array
h_edge => s % h_edge % array
- tend_h => s % h % array
- tend_u => s % u % array
circulation => s % circulation % array
vorticity => s % vorticity % array
divergence => s % divergence % array
@@ -1064,18 +1045,6 @@
end do
!
- ! Compute vertical velocity in each cell
- ! w(k,iCell) is w_{k-1/2}, the velocity through the interface above cell k
- ! Upon entry, w(k,iCell) is F^v_{k-1/2}, the vertical thickness flux,
- ! computed in compute_tend.
- !
-! print '(a,i5,f10.5)', &
-! 'k, min(w(k,2:)),max(w(k,2:))'
-! do k=1,nVertLevels
-! print '(i5,10es10.2)', &
-! k,minval(w(k,1:grid % ncellssolve)),maxval(w(k,1:grid % ncellssolve))
-! enddo
- !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
@@ -1269,8 +1238,6 @@
enddo
- !mrp 100426:
- !
! For z-level model.
! Compute pressure at middle of each level.
! pmidZLevel and pmid should only differ at k=1, where pmid is
@@ -1290,8 +1257,62 @@
end do
end do
- !mrp 100426 end
+ ! compute vertical velocity
+ if (config_vert_grid_type.eq.'isopycnal') then
+ ! set vertical velocity to zero in isopycnal case
+ wBot=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))
+ div_u(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells) then
+ do k=1,nVertLevels
+! should be this
+ flux = u(k,iEdge) * dvEdge(iEdge)
+! testing with this
+! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ div_u(k,cell1) = div_u(k,cell1) + flux
+ end do
+ end if
+ if (cell2 <= nCells) then
+ do k=1,nVertLevels
+! should be this
+ flux = u(k,iEdge) * dvEdge(iEdge)
+! testing with this
+! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ div_u(k,cell2) = div_u(k,cell2) - flux
+ end do
+ end if
+ end do
+
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
+ end do
+
+ ! change nvertlevels to kmt later
+ ! this next line can be set permanently somewhere upon startup.
+ wBot(nVertLevels,iCell) = 0.0
+
+ do k=nVertLevels,2,-1
+ wBot(k-1,iCell) = wBot(k,iCell) - div_u(k,iCell)*h(k,iCell)
+ end do
+
+ end do
+
+ endif
+ deallocate(div_u)
+
+
end subroutine compute_solve_diagnostics
</font>
</pre>