<p><b>mpetersen@lanl.gov</b> 2010-04-27 09:43:25 -0600 (Tue, 27 Apr 2010)</p><p>Z-level branch now has a config_vert_grid_type=isopycnal or zlevel namelist flag. When zlevel is chosen, Eulerian grid vertical fluxes are computed, and the horizontal pressure gradient grad p/rho0 is used. When isopycnal is chosen, vertical fluxes are zero, and the Montgomery Potential is used.<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-26 22:58:02 UTC (rev 213)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-04-27 15:43:25 UTC (rev 214)
@@ -6,9 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
-# mrp 100120:
namelist integer sw_model config_stats_interval 100
-# mrp 100120 end
namelist real sw_model config_visc 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -16,6 +14,10 @@
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
+# config_vert_grid_type: zlevel or isopycnal
+namelist character grid config_vert_grid_type isopycnal
+# mrp 100426 added rho0, only used in zlevel pgrad term.
+namelist real grid config_rho0 1028
#
# dim type name_in_file name_in_code
@@ -82,6 +84,13 @@
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+# Arrays for z-level version of mpas-ocean
+var integer kmaxCell ( nCells ) iro kmaxCell - -
+var integer kmaxEdge ( nEdges ) iro kmaxEdge - -
+var real hZLevel ( nVertLevels ) iro hZLevel - -
+var real zmidZLevel ( nVertLevels ) iro zmidZLevel - -
+var real zbotZLevel ( nVertLevels ) iro zbotZLevel - -
+
# Boundary conditions: read from input, saved in restart and written to output
var integer uBC ( nVertLevels nEdges ) iro uBC - -
var real u_src ( nVertLevels nEdges ) iro u_src - -
@@ -89,7 +98,7 @@
# Prognostic variables: read from input, saved in restart, and written to output
var real u ( nVertLevels nEdges Time ) iro u - -
var real h ( nVertLevels nCells Time ) iro h - -
-var real rho ( nVertLevels nCells Time ) io rho - -
+var real rho ( nVertLevels nCells Time ) iro rho - -
var real temperature ( nVertLevels nCells Time ) iro temperature tracers dynamics
var real salinity ( nVertLevels nCells Time ) iro salinity tracers dynamics
var real tracer1 ( nVertLevels nCells Time ) iro tracer1 tracers testing
@@ -113,6 +122,7 @@
var real zSurface ( nCells Time ) o zSurface - -
var real pmid ( nVertLevels nCells Time ) o pmid - -
var real pbot ( nVertLevels nCells Time ) o pbot - -
+var real pmidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
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_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-04-26 22:58:02 UTC (rev 213)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-04-27 15:43:25 UTC (rev 214)
@@ -23,7 +23,8 @@
integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
- real (kind=RKIND), dimension(:), pointer :: xCell,yCell
+ real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
+ hZLevel, zmidZLevel, zbotZLevel
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
@@ -79,7 +80,7 @@
end if
! mrp 100121:
- !
+ !
! Initialize u_src, rho, alpha
! This is a temporary fix until everything is specified correctly
! in an initial conditions file.
@@ -95,6 +96,10 @@
xCell => block_ptr % mesh % xCell % array
yCell => block_ptr % mesh % yCell % array
+ hZLevel => block_ptr % mesh % hZLevel % array
+ zmidZLevel => block_ptr % mesh % zmidZLevel % array
+ zbotZLevel => block_ptr % mesh % zbotZLevel % array
+
nCells = block_ptr % mesh % nCells
nEdges = block_ptr % mesh % nEdges
nVertices = block_ptr % mesh % nVertices
@@ -118,12 +123,14 @@
!tracers(2,iLevel,iCell) = 1.0
!tracers(3,iLevel,iCell) = xCell(iCell)
!if (xCell(iCell)<500*1e3) tracers(4,iLevel,iCell) = 1.0
-! tracers(2,iLevel,iCell) = 5.0 ! temperature
-! tracers(3,iLevel,iCell) = 35.0 ! salinity
- tracers(index_temperature,iLevel,iCell) = &
- 10.0* xCell(iCell)/2500.e3
- tracers(index_salinity,iLevel,iCell) = &
- 34.0 + 2* yCell(iCell)/4000.e3
+! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
+! tracers(index_temperature,iLevel,iCell) = 6.0-iLevel ! temperature
+ tracers(index_temperature,iLevel,iCell) = 6.0-iLevel+yCell(iCell)/4000.e3 ! temperature
+ tracers(index_salinity,iLevel,iCell) = 35.0 ! salinity
+! tracers(index_temperature,iLevel,iCell) = &
+! 10.0* xCell(iCell)/2500.e3
+! tracers(index_salinity,iLevel,iCell) = &
+! 34.0 + 2* yCell(iCell)/4000.e3
tracers(index_tracer1,iLevel,iCell) = 1.0
tracers(index_tracer2,iLevel,iCell) = &
(yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
@@ -156,6 +163,29 @@
print '(10a)', 'EXPAND_LEVELS compiler flag is off.'
#endif
+
+ ! mrp 100426: initialize z-level variables.
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ hZLevel = h(:,1)
+ zmidZLevel(1) = -0.5*hZLevel(1)
+ zbotZLevel(1) = -hZLevel(1)
+ do iLevel = 2,nVertLevels
+ zmidZLevel(iLevel) = zbotZLevel(iLevel-1)-0.5*hZLevel(iLevel)
+ zbotZLevel(iLevel) = zbotZLevel(iLevel-1)- hZLevel(iLevel)
+ enddo
+ if (config_vert_grid_type.eq.'isopycnal') then
+ print *, ' Using isopycnal coordinates'
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ print *, ' Using z-level coordinates'
+ else
+ print *, ' Incorrect choice of config_vert_grid_type:',&
+ config_vert_grid_type
+ stop
+ endif
+
+ ! mrp 100426 end: initialize z-level variables.
+
do i=2,nTimeLevs
call copy_state(block_ptr % time_levs(1) % state, &
block_ptr % time_levs(i) % state)
@@ -165,13 +195,20 @@
print '(10a)', 'ilevel',&
' rho ',&
' min u max u ',&
+ ' min T max T ',&
+ ' min S max S ',&
+ ' min u_src max u_src ', &
' min h max h ',&
- ' min u_src max u_src '
+ ' hZLevel zmidZlevel', &
+ ' zbotZlevel'
do iLevel = 1,nVertLevels
print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
minval(u(iLevel,:)), maxval(u(iLevel,:)), &
+ minval(tracers(index_temperature,iLevel,:)), maxval(tracers(index_temperature,iLevel,:)), &
+ minval(tracers(index_salinity,iLevel,:)), maxval(tracers(index_salinity,iLevel,:)), &
+ minval(u_src(iLevel,:)), maxval(u_src(iLevel,:)), &
minval(h(iLevel,:)), maxval(h(iLevel,:)), &
- minval(u_src(iLevel,:)), maxval(u_src(iLevel,:))
+ hZLevel(iLevel),zmidZlevel(iLevel),zbotZlevel(iLevel)
enddo
block_ptr => block_ptr % next
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-26 22:58:02 UTC (rev 213)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-04-27 15:43:25 UTC (rev 214)
@@ -283,14 +283,14 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wHat, areaSumInv
+ upstream_bias, wHat, areaSumInv, rho0Inv
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudz
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
+ divergence, MontPot, pmidZLevel
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
@@ -317,6 +317,7 @@
vh => s % vh % array
!mrp 100112:
MontPot => s % MontPot % array
+ pmidZLevel => s % pmidZLevel % array
!mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
@@ -380,47 +381,48 @@
!end do
! mrp 100409 replace with
- if (1==2) then ! isopycnal coordinates
- 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
+ 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
- else ! 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
+ 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
- ! change nvertlevels to kmt later
- ! this next line can be set permanently somewhere upon startup.
- Fv(nVertLevels,iCell) = 0.0
+ ! change nvertlevels to kmt later
+ ! this next line can be set permanently somewhere upon startup.
+ Fv(nVertLevels,iCell) = 0.0
- do k=nVertLevels,2,-1
+ do k=nVertLevels,2,-1
- ! F^v_{k-1/2} =
- ! F^v_{k+1/2} - </font>
<font color="red">abla_h \cdot \left( F^h \right) h_k
- ! note +tend_h because tend_h is -div.(hu)
- Fv(k-1,iCell) = Fv(k,iCell) + tend_h(k,iCell)*h(k,iCell)
+ ! F^v_{k-1/2} =
+ ! F^v_{k+1/2} - </font>
<font color="red">abla_h \cdot \left( F^h \right) h_k
+ ! note +tend_h because tend_h is -div.(hu)
+ Fv(k-1,iCell) = Fv(k,iCell) + tend_h(k,iCell)*h(k,iCell)
- ! now tend_h is div.(hu) + d/dz(hw):
- ! - </font>
<font color="red">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) &
+ ! 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-1,iCell) - Fv(k,iCell))/h(k,iCell)
- end do
+ end do
- k=1
- ! Surface tracer Fv(k,iCell) = 0.0 is not used
- tend_h(k,iCell) = tend_h(k,iCell) &
- - ( - Fv(k,iCell))/h(k,iCell)
+ k=1
+ ! Surface tracer Fv(k,iCell) = 0.0 is not used
+ tend_h(k,iCell) = tend_h(k,iCell) &
+ - ( - Fv(k,iCell))/h(k,iCell)
- end do
+ end do
endif ! coordinate type
! mrp 100409 end
@@ -438,6 +440,7 @@
! Compute u (normal) velocity tendency for each edge (cell face)
!
tend_u(:,:) = 0.0
+ rho0Inv = 1.0/config_rho0
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -465,6 +468,21 @@
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
+ 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
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - rho0Inv*( pmidZLevel(k,cell2) &
+ - pmidZLevel(k,cell1) )/dcEdge(iEdge)
+ end do
+ endif
+ ! mrp 100426 end: add pressure gradient
+
do k=1,nVertLevels
!
@@ -481,21 +499,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 100112, this is the original shallow water formulation with grad H:
- !tend_u(k,iEdge) = &
- ! q &
- ! + u_diffusion &
- ! - ( ke(k,cell2) - ke(k,cell1) &
- ! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ! ) / dcEdge(iEdge)
- ! mrp 100112, changed to grad Montgomery potential:
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 100112 end
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
!ocean
tend_u(k,iEdge) = tend_u(k,iEdge) + u_src(k,iEdge)/rho_ref/h_edge(k,iEdge)
@@ -670,17 +677,18 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, pbot_temp
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zSurface
+ zSurface, hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, &
tend_h, tend_u, circulation, vorticity, ke, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- zmid, zbot, pmid, pbot, MontPot, rho, temperature, salinity
+ zmid, zbot, pmid, pbot, MontPot, rho, temperature, salinity, pmidZLevel
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
character :: c1*6
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
@@ -707,11 +715,11 @@
gradPVt => s % gradPVt % array
!mrp 100112:
rho => s % rho % array
- temperature => s % tracers % array(index_temperature,:,:)
- salinity => s % tracers % array(index_salinity,:,:)
+ tracers => s % tracers % array
zmid => s % zmid % array
zbot => s % zbot % array
pmid => s % pmid % array
+ pmidZLevel => s % pmidZLevel % array
pbot => s % pbot % array
MontPot => s % MontPot % array
zSurface => s % zSurface % array
@@ -734,6 +742,7 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
+ hZLevel => grid % hZLevel % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -988,19 +997,20 @@
!
! equation of state
!
-! do iCell=1,nCells
-! do k=1,nVertLevels
-! ! Linear equation of state, for the time being
-! rho(k,iCell) = 1000.0*( 1.0 &
-! - 2.5e-4*temperature(k,iCell) &
-! + 7.6e-4*salinity(k,iCell))
-! end do
-! end do
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ ! Linear equation of state, for the time being
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(index_temperature,k,iCell) &
+ + 7.6e-4*tracers(index_salinity,k,iCell))
+ end do
+ end do
!mrp 100324 end
!mrp 100112:
!
+ ! For Isopycnal model.
! Compute mid- and bottom-depth of each layer, from bottom up.
! Then compute mid- and bottom-pressure of each layer, and
! Montgomery Potential, from top down
@@ -1019,9 +1029,9 @@
! rather than using zbot(0,iCell), I am adding another 2D variable.
zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
- ! assume pressure at the surface is zero for now.
- pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
- pbot(1,iCell) = rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
+ ! assume atmospheric pressure at the surface is zero for now.
+ pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell)
+ pbot(1,iCell) = rho(1,iCell)*gravity* h(1,iCell)
MontPot(1,iCell) = gravity * zSurface(iCell)
do k=2,nVertLevels
delta_p = rho(k,iCell)*gravity* h(k,iCell)
@@ -1032,9 +1042,34 @@
MontPot(k,iCell) = MontPot(k-1,iCell) &
+ pbot(k-1,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
+
end do
!mrp 100112 end
+
+ !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
+ ! pressure at middle of layer including SSH, and pmidZLevel is
+ ! pressure at a depth of hZLevel(1)/2.
+ !
+ do iCell=1,nCells
+ ! compute pressure for z-level coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
+ pmidZLevel(1,iCell) = rho(1,iCell)*gravity &
+ * (h(1,iCell)-0.5*hZLevel(1))
+ pbot_temp = rho(1,iCell)*gravity*h(1,iCell)
+ do k=2,nVertLevels
+ delta_p = rho(k,iCell)*gravity*hZLevel(k)
+ pmidZLevel(k,iCell) = pmidZLevel(k-1,iCell) + 0.5*delta_p
+ pbot_temp = pbot_temp + delta_p
+ end do
+
+ end do
+ !mrp 100426 end
+
end subroutine compute_solve_diagnostics
</font>
</pre>