<p><b>mpetersen@lanl.gov</b> 2011-12-06 14:35:39 -0700 (Tue, 06 Dec 2011)</p><p>Continue with changes required for ALE in vertical:<br>
<br>
4. Change pressure gradient calculation to not involve hZLevel.<br>
<br>
5. add - rho g /rho_0 grad z_k^{mid} term to the pressure gradient of the momentum equation.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/Registry 2011-12-06 21:00:00 UTC (rev 1239)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/Registry 2011-12-06 21:35:39 UTC (rev 1240)
@@ -216,6 +216,7 @@
var persistent real u_diffusionBtr ( nEdges Time ) 1 - u_diffusionBtr state - -
# Diagnostic fields: only written to output
+var persistent real zMid ( nVertLevels nCells Time ) 2 io zMid state - -
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F 2011-12-06 21:00:00 UTC (rev 1239)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F 2011-12-06 21:35:39 UTC (rev 1240)
@@ -253,7 +253,7 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, rho, zMid, pressure, &
tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
MontPot, wTop, divergence, vertViscTopOfEdge
type (dm_info) :: dminfo
@@ -279,7 +279,9 @@
h => s % h % array
u => s % u % array
v => s % v % array
+ rho => s % rho % array
wTop => s % wTop % array
+ zMid => s % zMid % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -359,9 +361,9 @@
call mpas_timer_start("ocn_tend_u-pressure grad")
if (config_vert_grid_type.eq.'isopycnal') then
- call ocn_vel_pressure_grad_tend(grid, MontPot, tend_u, err)
+ call ocn_vel_pressure_grad_tend(grid, MontPot, zMid, rho, tend_u, err)
elseif (config_vert_grid_type.eq.'zlevel') then
- call ocn_vel_pressure_grad_tend(grid, pressure, tend_u, err)
+ call ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend_u, err)
end if
call mpas_timer_stop("ocn_tend_u-pressure grad")
@@ -441,7 +443,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge, vertDiffTopOfCell
+ u,h,rho, wTop, h_edge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:,:), pointer :: boundaryEdge
@@ -469,6 +471,7 @@
u => s % u % array
h => s % h % array
+ rho => s % rho % array
boundaryCell=> grid % boundaryCell % array
wTop => s % wTop % array
tracers => s % tracers % array
@@ -615,10 +618,10 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel
+ hZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
- circulation, vorticity, ke, ke_edge, MontPot, wTop, &
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
rho, temperature, salinity, kev, kevc
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -660,6 +663,7 @@
tracers => s % tracers % array
MontPot => s % MontPot % array
pressure => s % pressure % array
+ zMid => s % zMid % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -679,6 +683,7 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
@@ -1092,15 +1097,36 @@
! compute pressure for z-level coordinates
! assume atmospheric pressure at the surface is zero for now.
+! mrp 111206 pressure for z-level, old
+! pressure(1,iCell) = rho(1,iCell)*gravity &
+! * (h(1,iCell)-0.5*hZLevel(1))
+! do k=2,maxLevelCell(iCell)
+! pressure(k,iCell) = pressure(k-1,iCell) &
+! + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+! + rho(k ,iCell)*hZLevel(k ))
+! end do
+
+ ! pressure for generalized coordinates
pressure(1,iCell) = rho(1,iCell)*gravity &
- * (h(1,iCell)-0.5*hZLevel(1))
+ * 0.5*h(1,iCell)
do k=2,maxLevelCell(iCell)
pressure(k,iCell) = pressure(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
- + rho(k ,iCell)*hZLevel(k ))
+ + 0.5*gravity*( rho(k-1,iCell)*h(k-1,iCell) &
+ + rho(k ,iCell)*h(k ,iCell))
end do
+ ! Compute zMid, the z-coordinate of the middle of the layer.
+ k = maxLevelCell(iCell)
+ zMid(k:nVertLevels,iCell) = zTopZLevel(k+1) + 0.5*h(k,iCell)
+
+ do k=maxLevelCell(iCell)-1, 1, -1
+ zMid(k,iCell) = zMid(k+1,iCell) &
+ + 0.5*( h(k+1,iCell) &
+ + h(k ,iCell))
+ end do
+
+
end do
endif
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2011-12-06 21:00:00 UTC (rev 1239)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2011-12-06 21:35:39 UTC (rev 1240)
@@ -17,6 +17,7 @@
use mpas_grid_types
use mpas_configure
+ use mpas_constants
implicit none
private
@@ -43,7 +44,7 @@
!
!--------------------------------------------------------------------
- real (kind=RKIND) :: rho0Inv
+ real (kind=RKIND) :: rho0Inv, grho0Inv
!***********************************************************************
@@ -64,7 +65,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_pressure_grad_tend(grid, pressure, tend, err)!{{{
+ subroutine ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -73,7 +74,9 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
- pressure !< Input: Pressure field or Mongomery potential
+ pressure, & !< Input: Pressure field or Mongomery potential
+ zMid, & !< Input: z-coordinate at mid-depth of layer
+ rho !< Input: density
type (mesh_type), intent(in) :: &
grid !< Input: grid information
@@ -123,15 +126,33 @@
- (pressure(k,cell2) - pressure(k,cell1))/dcEdge(iEdge)
end do
enddo
- elseif (config_vert_grid_type.eq.'zlevel') then
+ else
+! mrp 111206 pressure for z-level, old, delete before trunk merge.
+! do iEdge=1,nEdgesSolve
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
+! tend(k,iEdge) = tend(k,iEdge) &
+! - rho0Inv*( pressure(k,cell2) &
+! - pressure(k,cell1) )/dcEdge(iEdge)
+! end do
+! enddo
+
+ ! pressure for generalized coordinates
+ ! -1/rho_0 (grad p_k + rho g grad z_k^{mid})
+ grho0Inv = gravity*rho0Inv
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
-
tend(k,iEdge) = tend(k,iEdge) &
- rho0Inv*( pressure(k,cell2) &
- - pressure(k,cell1) )/dcEdge(iEdge)
+ - pressure(k,cell1) )/dcEdge(iEdge) &
+ - grho0Inv* 0.5*(rho(k,cell1)+rho(k,cell2)) &
+ *( zMid(k,cell2) &
+ - zMid(k,cell1) )/dcEdge(iEdge)
+
end do
enddo
@@ -178,11 +199,8 @@
err = 0
- if (config_vert_grid_type.eq.'isopycnal') then
- rho0Inv = 1.0
- elseif (config_vert_grid_type.eq.'zlevel') then
- rho0Inv = 1.0/config_rho0
- end if
+ rho0Inv = 1.0/config_rho0
+ grho0Inv = gravity/config_rho0
!--------------------------------------------------------------------
</font>
</pre>