<p><b>mpetersen@lanl.gov</b> 2011-12-07 08:34:55 -0700 (Wed, 07 Dec 2011)</p><p>I now have a working zstar setting in this branch. Changes in this commit:<br>
<br>
6. Change h_tend to always solve for full equation, regardless of vertical coord.<br>
7. change vel_vadv to not use zlevel variables<br>
8. change<br>
elseif (config_vert_grid_type.eq.'zlevel') then<br>
to<br>
else<br>
to accomodate zstar, ztilde, etc.<br>
9. Add branch for config_vert_grid_type.eq.'zstar' in wtop computation.<br>
Weights are proportional to layer thickness.<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 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/Registry 2011-12-07 15:34:55 UTC (rev 1242)
@@ -174,9 +174,6 @@
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
-var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
-var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
-var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
@@ -188,7 +185,7 @@
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
-var persistent real h ( nVertLevels nCells Time ) 2 ir h state - -
+var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
var persistent real rho ( nVertLevels nCells Time ) 2 iro rho state - -
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
var persistent real salinity ( nVertLevels nCells Time ) 2 iro salinity state tracers dynamics
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_equation_of_state.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_equation_of_state.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_equation_of_state.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -159,7 +159,7 @@
linearEos = .false.
jmEos = .false.
- if(config_vert_grid_type.eq.'zlevel') then
+ if(config_vert_grid_type.ne.'isopycnal') then
eosON = .true.
if (config_eos_type.eq.'linear') then
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -88,10 +88,17 @@
call compute_maxLevel(domain)
if (config_vert_grid_type.eq.'isopycnal') then
- print *, ' Using isopycnal coordinates'
+ print *, ' Using isopycnal vertical coordinates'
elseif (config_vert_grid_type.eq.'zlevel') then
- print *, ' Using z-level coordinates'
+ print *, ' Using z-level vertical coordinates'
call init_ZLevel(domain)
+ elseif (config_vert_grid_type.eq.'zstar1') then
+ print *, ' Using z-star vertical coordinates,', &
+ ' all convergence in top layer'
+ call init_ZLevel(domain)
+ elseif (config_vert_grid_type.eq.'zstar') then
+ print *, ' Using z-star vertical coordinates'
+ call init_ZLevel(domain)
else
print *, ' Incorrect choice of config_vert_grid_type:',&
config_vert_grid_type
@@ -474,8 +481,7 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
real (kind=RKIND), dimension(:), pointer :: &
- hZLevel, zMidZLevel, zTopZLevel, &
- hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
+ hZLevel, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -488,9 +494,6 @@
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
nVertLevels = block % mesh % nVertLevels
- hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
- hRatioZLevelK => block % mesh % hRatioZLevelK % array
- hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
! These should eventually be in an input file. For now
! I just read them in from h(:,1).
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 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -351,7 +351,7 @@
!
call mpas_timer_start("ocn_tend_u-vert adv")
- call ocn_vel_vadv_tend(grid, u, wtop, tend_u, err)
+ call ocn_vel_vadv_tend(grid, u, h_edge, wtop, tend_u, err)
call mpas_timer_stop("ocn_tend_u-vert adv")
@@ -362,7 +362,7 @@
if (config_vert_grid_type.eq.'isopycnal') then
call ocn_vel_pressure_grad_tend(grid, MontPot, zMid, rho, tend_u, err)
- elseif (config_vert_grid_type.eq.'zlevel') then
+ else
call ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend_u, err)
end if
@@ -1043,7 +1043,7 @@
!
! For an isopycnal model, density should remain constant.
! For zlevel, calculate in-situ density
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.ne.'isopycnal') then
call ocn_equation_of_state_rho(s, grid, 0, 'relative', err)
! mrp 110324 In order to visualize rhoDisplaced, include the following
call ocn_equation_of_state_rho(s, grid, 1, 'relative', err)
@@ -1146,15 +1146,16 @@
! mrp 110512 could clean this out, remove pointers?
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: u,wTop, h_edge
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
+ real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge
+ real (kind=RKIND), dimension(:,:), allocatable:: div_hu
+ real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
@@ -1166,6 +1167,7 @@
call mpas_timer_start("wTop")
u => s % u % array
+ h => s % h % array
wTop => s % wTop % array
h_edge => s % h_edge % array
@@ -1179,7 +1181,33 @@
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
+ allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
+ h_tend_col(nVertLevels))
+
!
+ ! Compute div(h^{edge} u) for each cell
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ !
+ div_hu(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,maxLevelEdgeBot(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ div_hu(k,cell1) = div_hu(k,cell1) + flux
+ div_hu(k,cell2) = div_hu(k,cell2) - flux
+ end do
+ end do
+
+ do iCell=1,nCells
+ div_hu_btr(iCell) = 0.0
+ do k=1,maxLevelCell(iCell)
+ div_hu(k,iCell) = div_hu(k,iCell) / areaCell(iCell)
+ div_hu_btr(iCell) = div_hu_btr(iCell) + div_hu(k,iCell)
+ end do
+ end do
+
+ !
! vertical velocity through layer interface
!
if (config_vert_grid_type.eq.'isopycnal') then
@@ -1188,36 +1216,64 @@
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+1))
- div_u(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- div_u(k,cell1) = div_u(k,cell1) + flux
- div_u(k,cell2) = div_u(k,cell2) - flux
- end do
- end do
+ do iCell=1,nCells
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell)
+ end do
+ end do
+ elseif (config_vert_grid_type.eq.'zstar1') then
+
+ ! This is a testing setting. The computation is similar to zstar,
+ ! but the weights are all in the top layer, so is a bit-for-bit
+ ! match with zlevel.
+
do iCell=1,nCells
+
+ h_tend_col = 0.0
+ h_tend_col(1) = - div_hu_btr(iCell)
+
! Vertical velocity through layer interface at top and
! bottom is zero.
wTop(1,iCell) = 0.0
wTop(maxLevelCell(iCell)+1,iCell) = 0.0
do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) &
- - div_u(k,iCell)/areaCell(iCell)
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
end do
end do
- deallocate(div_u)
+ elseif (config_vert_grid_type.eq.'zstar') then
+
+ ! Distribute the change in total column height due to the external
+ ! mode, div_hu_btr, among all the layers. Distribute in proportion
+ ! to the layer thickness.
+
+ do iCell=1,nCells
+
+ hSum = 0.0
+ do k=1,maxLevelCell(iCell)
+ h_tend_col(k) = - h(k,iCell)*div_hu_btr(iCell)
+ hSum = hSum + h(k,iCell)
+ end do
+ h_tend_col = h_tend_col / hSum
+
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
+ end do
+ end do
+
endif
+ deallocate(div_hu, div_hu_btr, h_tend_col)
+
call mpas_timer_stop("wTop")
end subroutine ocn_wtop!}}}
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_hadv.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_hadv.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -104,7 +104,7 @@
integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
integer :: iCell, nCells
- integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND) :: flux
@@ -124,46 +124,27 @@
nCells = grid % nCells
nVertLevels = grid % nVertLevels
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
cellsOnEdge => grid % cellsOnEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- if (config_vert_grid_type.eq.'isopycnal') then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux
- tend(k,cell2) = tend(k,cell2) + flux
- end do
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend(k,cell1) = tend(k,cell1) - flux
+ tend(k,cell2) = tend(k,cell2) + flux
end do
- do iCell=1,nCells
- do k=1,nVertLevels
- tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
- end do
+ end do
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
end do
+ end do
- elseif (config_vert_grid_type.eq.'zlevel') then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,min(1,maxLevelEdgeTop(iEdge))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux
- tend(k,cell2) = tend(k,cell2) + flux
- end do
- end do
- do iCell=1,nCells
- tend(1,iCell) = tend(1,iCell) / areaCell(iCell)
- end do
-
- endif ! config_vert_grid_type
-
-
!--------------------------------------------------------------------
end subroutine ocn_thick_hadv_tend!}}}
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_vadv.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_vadv.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_thick_vadv.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -98,7 +98,8 @@
!
!-----------------------------------------------------------------
- integer :: iCell, nCells
+ integer :: iCell, nCells, nVertLevels, k
+ integer, dimension(:), pointer :: MaxLevelCell
!-----------------------------------------------------------------
!
@@ -110,13 +111,17 @@
err = 0
+ maxLevelCell => grid % maxLevelCell % array
+
nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
do iCell=1,nCells
- tend(1,iCell) = tend(1,iCell) + wTop(2,iCell)
+ do k=1,maxLevelCell(iCell)
+ tend(k,iCell) = tend(k,iCell) + wTop(k+1,iCell) - wTop(k,iCell)
+ end do
end do
-
!--------------------------------------------------------------------
end subroutine ocn_thick_vadv_tend!}}}
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -167,7 +167,7 @@
err = 0
vadvOn = .false.
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.ne.'isopycnal') then
vadvOn = .true.
call ocn_tracer_vadv_stencil_init(err1)
call ocn_tracer_vadv_spline_init(err2)
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_vadv.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_vadv.F 2011-12-06 22:52:00 UTC (rev 1241)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_vadv.F 2011-12-07 15:34:55 UTC (rev 1242)
@@ -64,7 +64,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_vadv_tend(grid, u, wTop, tend, err)!{{{
+ subroutine ocn_vel_vadv_tend(grid, u, h_edge, wTop, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -75,6 +75,7 @@
real (kind=RKIND), dimension(:,:), intent(in) :: &
u !< Input: Horizontal velocity
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge,&!< Input: thickness at edge
wTop !< Input: Vertical velocity on top layer
type (mesh_type), intent(in) :: &
@@ -134,7 +135,7 @@
! compute dudz at vertical interface with first order derivative.
w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
- / (zMidZLevel(k-1) - zMidZLevel(k))
+ / (0.5*(h_edge(k-1,iEdge) + h_edge(k,iEdge)))
end do
w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
! Average w*du/dz from vertical interface to vertical middle of cell
@@ -179,7 +180,7 @@
err = 0
velVadvOn = .false.
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.ne.'isopycnal') then
velVadvOn = .true.
end if
</font>
</pre>