<p><b>mpetersen@lanl.gov</b> 2012-04-11 09:22:26 -0600 (Wed, 11 Apr 2012)</p><p>Merging zstar_restart_new branch onto trunk.<br>
<br>
Added the ability to restart using zstar when last run used zlevel. The following changes were made:<br>
<br>
1. config_enforce_zstar_at_restart flag, must be true for thickness change to occur.<br>
2. added subroutine ocn_init_h_zstar where thickness are recalculated.<br>
3. Had to split ocn_init_z_level into two subroutines:<br>
ocn_init_z_level called before ocn_init_h_zstar<br>
ocn_init_split_timestep requires h field, so called after ocn_init_h_zstar<br>
4. Changed zstar weights to be a new Registry variable, zstarWeight<br>
5. Consolidated wTop computation so it is identical for all z-type coordinates,<br>
using only zstarWeight to differentiate among grid types.<br>
<br>
Tested using 120km global grid for one day, bit-for-bit matches between trunk and revised branch for zlevel and zstarWeights. zstar agrees to 7 digits, because weights changed from h to hZLevel.<br>
<br>
Ran 120km for 1 day with zlevel, then restarted using zstar. There is a slight (1e-5) adjustment visible in global KE, but smooths out after 0.5 days (p87k-n).<br>
<br>
Reviewed briefly by Todd using a 15km global grid, restarting from zlevel, testing both zlevel and zstar.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+ /branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/ocean_projects/zstar_restart_new:1762-1770
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2012-04-11 00:51:17 UTC (rev 1770)
+++ trunk/mpas/src/core_ocean/Registry 2012-04-11 15:22:26 UTC (rev 1771)
@@ -24,6 +24,7 @@
namelist character grid config_vert_grid_type isopycnal
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
+namelist logical grid config_enforce_zstar_at_restart false
namelist integer split_explicit_ts config_n_ts_iter 2
namelist integer split_explicit_ts config_n_bcl_iter_beg 2
namelist integer split_explicit_ts config_n_bcl_iter_mid 2
@@ -185,6 +186,7 @@
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
+var persistent real zstarWeight ( nVertLevels ) 0 - zstarWeight mesh - -
% Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-04-11 00:51:17 UTC (rev 1770)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-04-11 15:22:26 UTC (rev 1771)
@@ -105,6 +105,12 @@
call ocn_init_z_level(domain)
+ if (config_enforce_zstar_at_restart) then
+ call ocn_init_h_zstar(domain)
+ endif
+
+ call ocn_init_split_timestep(domain)
+
print *, ' Vertical grid type is: ',config_vert_grid_type
if (config_vert_grid_type.ne.'isopycnal'.and. &
@@ -499,7 +505,7 @@
end subroutine mpas_timestep!}}}
subroutine ocn_init_z_level(domain)!{{{
- ! Initialize maxLevel and bouncary grid variables.
+ ! Initialize zlevel-type variables
use mpas_grid_types
use mpas_configure
@@ -513,7 +519,8 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, referenceBottomDepthTopOfCell
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, &
+ referenceBottomDepthTopOfCell, zstarWeight, hZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -525,6 +532,8 @@
h => block % state % time_levs(1) % state % h % array
referenceBottomDepth => block % mesh % referenceBottomDepth % array
referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
+ zstarWeight => block % mesh % zstarWeight % array
+ hZLevel => block % mesh % hZLevel % array
nVertLevels = block % mesh % nVertLevels
! mrp 120208 right now hZLevel is in the grid.nc file.
@@ -532,9 +541,9 @@
! as the defining variable instead, and will transition soon.
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
- referenceBottomDepth(1) = block % mesh % hZLevel % array(1)
+ referenceBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
+ referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
end do
! TopOfCell needed where zero depth for the very top may be referenced.
@@ -543,6 +552,65 @@
referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
end do
+ ! Initialization of zstarWeights. This determines how SSH perturbations
+ ! are distributed throughout the column.
+ if (config_vert_grid_type.eq.'zlevel') then
+
+ zstarWeight = 0.0
+ zstarWeight(1) = 1.0
+
+ elseif (config_vert_grid_type.eq.'zstar') then
+
+ do k = 1,nVertLevels
+ zstarWeight(k) = hZLevel(k)
+ enddo
+
+ elseif (config_vert_grid_type.eq.'zstarWeights') then
+
+ ! This is a test with other weights, just to make sure zstar functions
+ ! using variable weights.
+
+ zstarWeight = 0.0
+ zstarWeight(1:5) = 1.0
+ do k=1,10
+ zstarWeight(5+k) = 1.0-k*0.1
+ end do
+
+ endif
+
+ block => block % next
+ end do
+
+ end subroutine ocn_init_z_level!}}}
+
+ subroutine ocn_init_split_timestep(domain)!{{{
+ ! Initialize splitting variables
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell, cell1, cell2
+ real (kind=RKIND) :: uhSum, hSum, hEdge1
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ nVertLevels = block % mesh % nVertLevels
+
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
if (trim(config_time_integration) == 'unsplit_explicit') then
@@ -618,8 +686,66 @@
block => block % next
end do
- end subroutine ocn_init_z_level!}}}
+ end subroutine ocn_init_split_timestep!}}}
+ subroutine ocn_init_h_zstar(domain)!{{{
+ ! If changing from zlevel to zstar, compute h based on zstar weights,
+ ! where SSH is distributed through the layers. We only change h.
+ ! We do not remap the tracer variables, so this breaks total global
+ ! conservation.
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+
+ integer :: i, iCell, iEdge, iVertex, k, nVertLevels
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: hSum, sumZstarWeights
+ real (kind=RKIND), dimension(:), pointer :: hZLevel, zstarWeight, &
+ referenceBottomDepth
+ real (kind=RKIND), dimension(:,:), pointer :: h
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ nVertLevels = block % mesh % nVertLevels
+ hZLevel => block % mesh % hZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+ zstarWeight => block % mesh % zstarWeight % array
+ referenceBottomDepth => block % mesh % referenceBottomDepth % array
+
+ do iCell=1,block % mesh % nCells
+ ! Compute the total column thickness, hSum, and the sum of zstar weights.
+ hSum = 0.0
+ sumZstarWeights = 0.0
+ do k = 1,maxLevelCell(iCell)
+ hSum = hSum + h(k,iCell)
+ sumZstarWeights = sumZstarWeights + zstarWeight(k)
+ enddo
+
+ ! h_k = h_k^{zlevel} + zeta * W_k/sum(W_k)
+ ! where zeta is SSH and W_k are weights
+ do k = 1,maxLevelCell(iCell)
+ h(k,iCell) = hZLevel(k) &
+ + (hSum - referenceBottomDepth(maxLevelCell(iCell))) &
+ * zstarWeight(k)/sumZstarWeights
+ enddo
+
+ enddo
+
+ block => block % next
+ end do
+
+ end subroutine ocn_init_h_zstar!}}}
+
subroutine ocn_compute_max_level(domain)!{{{
! Initialize maxLevel and bouncary grid variables.
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-04-11 00:51:17 UTC (rev 1770)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-04-11 15:22:26 UTC (rev 1771)
@@ -876,10 +876,10 @@
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zstarWeight
real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
real (kind=RKIND), dimension(:,:), allocatable:: div_hu
- real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col, h_weights
+ real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
@@ -898,13 +898,14 @@
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
+ zstarWeight => grid % zstarWeight % array
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels), h_weights(nVertLevels))
+ h_tend_col(nVertLevels))
!
! Compute div(h^{edge} u) for each cell
@@ -937,51 +938,14 @@
! set vertical velocity to zero in isopycnal case
wTop=0.0
- elseif (config_vert_grid_type.eq.'zlevel') then
+ else ! zlevel or zstar type vertical grid
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_hu(k,iCell) - h_tend_col(k)
- end do
- end do
-
- 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)
+ h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
+ hSum = hSum + zstarWeight(k)*h(k,iCell)
end do
h_tend_col = h_tend_col / hSum
@@ -994,37 +958,9 @@
end do
end do
- elseif (config_vert_grid_type.eq.'zstarWeights') then
-
- ! This is a test with other weights, not meant to be permanent.
-
- h_weights = 0.0
- h_weights(1:5) = 1.0
- do k=1,10
- h_weights(5+k) = 1.0-k*0.1
- end do
-
- do iCell=1,nCells
-
- hSum = 0.0
- do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - h_weights(k)*h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + h_weights(k)*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, h_weights)
+ deallocate(div_hu, div_hu_btr, h_tend_col)
end subroutine ocn_wtop!}}}
</font>
</pre>