<p><b>mpetersen@lanl.gov</b> 2012-10-26 11:05:47 -0600 (Fri, 26 Oct 2012)</p><p>merge branch to trunk: restart_reproducibility<br>
<br>
The following three changes were required to have bit-for-bit restartibility:<br>
<br>
1. Call diagnostics before implicit vertical mixing (IVM)<br>
<br>
Variables like rho, h_edge, and ke_edge must be updated at the end of<br>
a timestep but before IVM. Right now we just added an additional call<br>
to ocn_diagnostics_solve just before IVM. After kpp is added, we<br>
could improve performance by only computing the diagnostics required<br>
for that particular IVM scheme beforehand. Diagnostics must be<br>
recomputed after IVM as well, since IVM alters u and tracers.<br>
<br>
This item prevented bit restartibility for RK4 and split explicit.<br>
<br>
2. Compute wTop with the correct advection velocity<br>
<br>
There are two velocities for advection:<br>
u, used by momentum equation<br>
uTransport = u + uCorr + uBolus, used by h and tracer equation<br>
<br>
wTop in the vertical advection must be computed with the correct<br>
horizontal velocity for each conservation equation. In this revision,<br>
I changed wTop to accept u as an input variable. wTop is called with<br>
the correct advection velocity just before computing the tendancies.<br>
<br>
This item prevented bit restartibility for split explicit because wTop<br>
was always computed with uTransport, except upon restart when it used<br>
u.<br>
<br>
3. Add uBtr to restart file.<br>
<br>
In split explicit timestepping, the velocity u is not resplit at the<br>
beginning of each timestep. Instead, we keep uBtr from the previous<br>
timestep. Note that uBcl may now include a barotropic component,<br>
because the weights h have changed. That is OK, because the<br>
GBtrForcing variable subtracts out the barotropic component from the<br>
baroclinic.<br>
<br>
This item prevented bit restartibility for split explicit because uBtr<br>
was computed upon restart. In could be solved by freshly splitting u<br>
at the beginning of each timestep. However, Stage 1 is set up to take<br>
care of that with the G term, and the additional splitting is not<br>
needed. Instead, we chose to save uBtr to the restart file for bit<br>
restartability, since this is only the addition of a 2D field.<br>
</p><hr noshade><pre><font color="gray">Index: trunk/mpas
===================================================================
--- trunk/mpas 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas 2012-10-26 17:05:47 UTC (rev 2273)
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
## -11,6 +11,7 ##
/branches/ocean_projects/monthly_forcing:1810-1867
/branches/ocean_projects/option3_b4b_test:2201-2231
/branches/ocean_projects/partial_bottom_cells:2172-2226
+/branches/ocean_projects/restart_reproducibility:2239-2272
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
\ No newline at end of property
Index: trunk/mpas/src/core_ocean
===================================================================
--- trunk/mpas/src/core_ocean 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean 2012-10-26 17:05:47 UTC (rev 2273)
Property changes on: trunk/mpas/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -11,6 +11,7 ##
/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
/branches/ocean_projects/partial_bottom_cells/src/core_ocean:2172-2226
+/branches/ocean_projects/restart_reproducibility/src/core_ocean:2239-2272
/branches/ocean_projects/split_explicit_mrp/src/core_ocean:1134-1138
/branches/ocean_projects/split_explicit_timestepping/src/core_ocean:1044-1097
/branches/ocean_projects/vert_adv_mrp/src/core_ocean:704-745
\ No newline at end of property
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean/Registry 2012-10-26 17:05:47 UTC (rev 2273)
@@ -206,21 +206,21 @@
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 - -
-var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
-var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+% Boundary conditions and masks
+var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 - boundaryEdge mesh - -
+var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 - boundaryVertex mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 - boundaryCell mesh - -
var persistent integer edgeMask ( nVertLevels nEdges ) 0 o edgeMask mesh - -
var persistent integer vertexMask ( nVertLevels nVertices ) 0 o vertexMask mesh - -
var persistent integer cellMask ( nVertLevels nCells ) 0 o cellMask mesh - -
-var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
-var persistent real temperatureRestore ( nCells ) 0 iro temperatureRestore mesh - -
-var persistent real salinityRestore ( nCells ) 0 iro salinityRestore mesh - -
-% Monthly forcing variables. Change to iro when monthly forcing is used.
-var persistent real windStressMonthly ( nMonths nEdges ) 0 - windStressMonthly mesh - -
-var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 - temperatureRestoreMonthly mesh - -
-var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 - salinityRestoreMonthly mesh - -
+% Forcing variables.
+var persistent real u_src ( nVertLevels nEdges ) 0 ir u_src mesh - -
+var persistent real temperatureRestore ( nCells ) 0 ir temperatureRestore mesh - -
+var persistent real salinityRestore ( nCells ) 0 ir salinityRestore mesh - -
+var persistent real windStressMonthly ( nMonths nEdges ) 0 ir windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 ir temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 ir salinityRestoreMonthly mesh - -
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
@@ -239,7 +239,7 @@
var persistent real tend_tracer1 ( nVertLevels nCells Time ) 1 - tracer1 tend tracers testing
% state variables for Split Explicit timesplitting
-var persistent real uBtr ( nEdges Time ) 2 - uBtr state - -
+var persistent real uBtr ( nEdges Time ) 2 ir uBtr state - -
var persistent real ssh ( nCells Time ) 2 o ssh state - -
var persistent real uBtrSubcycle ( nEdges Time ) 2 - uBtrSubcycle state - -
var persistent real sshSubcycle ( nCells Time ) 2 - sshSubcycle state - -
@@ -248,7 +248,7 @@
var persistent real uBcl ( nVertLevels nEdges Time ) 2 - uBcl state - -
% Diagnostic fields: only written to output
-var persistent real zMid ( nVertLevels nCells Time ) 2 io zMid state - -
+var persistent real zMid ( nVertLevels nCells Time ) 2 - zMid state - -
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real uTransport ( nVertLevels nEdges Time ) 2 - uTransport state - -
var persistent real uBolusGM ( nVertLevels nEdges Time ) 2 - uBolusGM state - -
Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-26 17:05:47 UTC (rev 2273)
@@ -114,7 +114,7 @@
call ocn_init_h_zstar(domain)
endif
- call ocn_init_split_timestep(domain)
+ if (.not.config_do_restart) call ocn_init_split_timestep(domain)
write (0,'(a,a10)'), ' Vertical grid type is: ',config_vert_grid_type
@@ -265,8 +265,6 @@
= block % state % time_levs(1) % state % u % array(:,:) &
+ block % state % time_levs(1) % state % uBolusGM % array(:,:)
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh)
-
call ocn_compute_mesh_scaling(mesh)
call mpas_rbf_interp_initialize(mesh)
@@ -715,7 +713,7 @@
enddo
if (.not. consistentSSH) then
- write(0,*) 'Warning: SSH is not consisntent'
+ write(0,*) 'Warning: SSH is not consistent. Most likely, initial h does not match bottomDepth.'
end if
endif
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-10-26 17:05:47 UTC (rev 2273)
@@ -854,13 +854,43 @@
!> This routine computes the vertical velocity in the top layer for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s1,s2, grid)!{{{
- implicit none
+ subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
- type (state_type), intent(inout) :: s1 !< Input/Output: State 1 information
- type (state_type), intent(inout) :: s2 !< Input/Output: State 2 information
- type (mesh_type), intent(in) :: grid !< Input: Grid information
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: h interpolated to an edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: velocity
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ wTop !< Output: vertical transport at top edge
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum, invAreaCell
@@ -869,7 +899,6 @@
real (kind=RKIND), dimension(:), pointer :: &
dvEdge, areaCell, zstarWeight
- real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
real (kind=RKIND), dimension(:), allocatable:: div_hu, h_tend_col
real (kind=RKIND) :: div_hu_btr
@@ -880,10 +909,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
- h => s1 % h % array
- h_edge => s1 % h_edge % array
- uTransport => s2 % uTransport % array
- wTop => s2 % wTop % array
+ err = 0
nEdgesOnCell => grid % nEdgesOnCell % array
areaCell => grid % areaCell % array
@@ -923,7 +949,7 @@
iEdge = edgesOnCell(i, iCell)
do k = 1, maxLevelEdgeBot(iEdge)
- flux = uTransport(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
div_hu(k) = div_hu(k) - flux
div_hu_btr = div_hu_btr - flux
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-10-26 17:05:47 UTC (rev 2273)
@@ -149,17 +149,19 @@
block => domain % blocklist
do while (associated(block))
- ! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(block % provis, block % provis, block % mesh)
-
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
end if
- call ocn_tend_h(block % tend, block % provis, block % mesh)
+
+ ! advection of u uses u, while advection of h and tracers use uTransport.
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % u % array, block % provis % wTop % array, err)
call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
- ! mrp 110718 filter btr mode out of u_tend
- ! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % uTransport % array, block % provis % wTop % array, err)
+ call ocn_tend_h(block % tend, block % provis, block % mesh)
+
if (config_rk_filter_btr_mode) then
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
@@ -273,6 +275,15 @@
call mpas_timer_start("RK4-implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
block => block % next
end do
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-25 23:36:34 UTC (rev 2272)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-26 17:05:47 UTC (rev 2273)
@@ -86,9 +86,9 @@
type (dm_info) :: dminfo
integer :: iCell, i,k,j, iEdge, cell1, cell2, split_explicit_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
- n_bcl_iter(config_n_ts_iter)
+ n_bcl_iter(config_n_ts_iter), stage1_tend_time
type (block_type), pointer :: block
- real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, &
+ real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, hEdge1, &
CoriolisTerm, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -117,6 +117,9 @@
! The baroclinic velocity needs be recomputed at the beginning of a
! timestep because the implicit vertical mixing is conducted on the
! total u. We keep uBtr from the previous timestep.
+ ! Note that uBcl may now include a barotropic component, because the
+ ! weights h have changed. That is OK, because the GBtrForcing variable
+ ! subtracts out the barotropic component from the baroclinic.
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
= block % state % time_levs(1) % state % u % array(k,iEdge) &
- block % state % time_levs(1) % state % uBtr % array( iEdge)
@@ -181,10 +184,22 @@
block => domain % blocklist
do while (associated(block))
+
+ stage1_tend_time = min(split_explicit_step,2)
+
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
end if
- call ocn_tend_u(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+
+ ! compute wTop. Use u (rather than uTransport) for momentum advection.
+ ! Use the most recent time level available.
+ call ocn_wtop(block % mesh, block % state % time_levs(stage1_tend_time) % state % h % array, &
+ block % state % time_levs(stage1_tend_time) % state % h_edge % array, &
+ block % state % time_levs(stage1_tend_time) % state % u % array, &
+ block % state % time_levs(stage1_tend_time) % state % wTop % array, err)
+
+ call ocn_tend_u(block % tend, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % mesh)
+
block => block % next
end do
@@ -246,6 +261,7 @@
= 0.5*( &
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+
enddo
enddo ! iEdge
@@ -767,8 +783,14 @@
! dwj: 02/22/12 splitting thickness and tracer tendency computations and halo updates to allow monotonic advection.
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
+ ! compute wTop. Use uTransport for advection of h and tracers.
+ ! Use time level 1 values of h and h_edge because h has not yet been computed for time level 2.
+ call ocn_wtop(block % mesh, block % state % time_levs(1) % state % h % array, &
+ block % state % time_levs(1) % state % h_edge % array, &
+ block % state % time_levs(2) % state % uTransport % array, &
+ block % state % time_levs(2) % state % wTop % array, err)
+
call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -921,7 +943,17 @@
call mpas_timer_start("se implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+
block => block % next
end do
@@ -959,6 +991,12 @@
= block % state % time_levs(2) % state % u % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
+!!!!!!!!!!!!!!!!!!!!!!!!!! mrp
+ ! Recompute wTop with all updated information. The wTop in stage 3 used
+ ! the previous h and u before implicit vertical mixing.
+!!!!!!
+! remove: call ocn_wtop(block % state % time_levs(2) % state,block % state % time_levs(2) % state, block % mesh)
+
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
</font>
</pre>