Added a hybrid restart option

Until now, the config option 'restart' has had two possible values:
* restart = 0 if not a restart (i.e., read in a CF input file and initialize the ice state)
* restart = 1 if a restart (i.e., read in a CF restart file that includes the full ice state)

This commit adds a new option, restart = 2, called a 'hybrid restart' because of its
similarity to a CESM hybrid run. We now refer to option 1 as a 'standard restart'.

A hybrid restart works as follows:

- The run is initialized from a file in the [CF input] section, as for restart = 0.
  However, this file has 'restart' or '.r.' in its name and includes the full ice state
  as needed for exact restart. Typically, it is the final restart time slice from a spin-up.
- The initial time (model%numerics%time) is set to 'tstart' as specified in the config file.
  This differs from a standard restart, which takes its initial time from the restart file.
- The initial tstep_count = 0.
  This differs from a standard restart, which takes tstep_count from the restart file.

For glaciers, we can use the hybrid restart option for commitment runs and other
forward runs starting from a spun-up state. For instance, say we want to do a 2000-year
spin-up followed by a forward run from 2003–2100. The workflow is as follows:
- Do the spin-up and write a final restart file.
- Set up a directory for the forward run with the required input and forcing files,
  including the restart file from the spin-up.
- In the config file:
  * Set tstart = 2003, tend = 2100., and restart = 2.
  * In the [CF input] section, set 'name' to the name of the restart file from the spin-up.
    This should be different from the name of the [CF restart] file for the forward run.
    E.g., one file could be spinup.restart.nc, and the other could be forward.restart.nc.
  * Change other options changes as needed, e.g. change the inversion options.
- Launch the forward run.

It is no longer necessary to modify the model time or tstep_count by hand in the restart file
from the spin-up. Also, it is not necessary to use the same name for (1) the restart file
from the spin-up and (2) the restart file for the forward run.

If the [CF output] file for the forward run is configured with 'write_init = .true.',
then this output file will include the ice state at the start of the forward run.
For a standard restart, CISM does not write to the output file at the start of the forward run.

Testing the new restart option, I confirmed that a hybrid restart is exact, as expected,
apart from the new values of the model time and tstep_count.

I also confirmed that the standard restart (restart = 1) works as before.

Author: whlipscomb

cism_driver/cism_front_end.F90

@@ -229,8 +229,10 @@ subroutine cism_init_dycore(model)
 
 
   ! --- Output the initial state -------------
+  ! Note: For a standard restart, the initial state is not output, because this state
+  !       should already have been written to the output file when the previous run ended.
 
-  if (model%options%is_restart == RESTART_FALSE .or. model%options%forcewrite_restart) then
+  if (model%options%is_restart == NO_RESTART .or. model%options%is_restart == HYBRID_RESTART) then
      call t_startf('initial_io_writeall')
      call glide_io_writeall(model, model, time=time)          ! MJH The optional time argument needs to be supplied 
                                                               !     since we have not yet set model%numerics%time

libglide/glide.F90

@@ -463,7 +463,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
        l_evolve_ice = .true.
     end if
 
-    if (model%options%is_restart == RESTART_TRUE) then
+    if (model%options%is_restart == STANDARD_RESTART .or. &
+        model%options%is_restart == HYBRID_RESTART) then
        ! On a restart, just assign the basal velocity from uvel/vvel (which are restart variables)
        ! to ubas/vbas which are used by the temperature solver to calculate basal heating.
        ! During time stepping ubas/vbas are calculated by slipvelo during thickness evolution or below on a cold start.

libglide/glide_diagnostics.F90

@@ -1165,7 +1165,7 @@ subroutine glide_write_diag (model,       time)
        write(message,'(a35,i14)') 'Diagnostic glacier index (CISM)    ', ng
        call write_log(trim(message), type = GM_DIAGNOSTIC)
 
-       write(message,'(a35,f14.6)') 'Glacier area_init(km^2)            ', &
+       write(message,'(a35,f14.6)') 'Glacier area_init (km^2)           ', &
             model%glacier%area_init(ng) / 1.0d6
        call write_log(trim(message), type = GM_DIAGNOSTIC)
 

libglide/glide_lithot.F90

@@ -82,7 +82,7 @@ subroutine init_lithot(model)
     !TODO - Make sure the sign is correct for the geothermal flux.
     !NOTE: CISM convention is that geot is positive down, so geot < 0 for upward geothermal flux
 
-    if (model%options%is_restart == RESTART_FALSE) then
+    if (model%options%is_restart == NO_RESTART) then
        ! set initial temp distribution to thermal gradient
        factor = model%paramets%geot / model%lithot%con_r
        do k=1,model%lithot%nlayer
@@ -112,7 +112,7 @@ subroutine spinup_lithot(model)
 
     integer t
 
-    if (model%options%is_restart == RESTART_FALSE .and. model%lithot%numt > 0) then
+    if (model%options%is_restart == NO_RESTART .and. model%lithot%numt > 0) then
        call write_log('Spinning up GTHF calculations',type=GM_INFO)
        call not_parallel(__FILE__,__LINE__)
        do t=1,model%lithot%numt

libglide/glide_setup.F90

@@ -766,7 +766,6 @@ subroutine handle_options(section, model)
     !      Going forward, only 'restart' is supported.
     call GetValue(section,'restart',model%options%is_restart)
     call GetValue(section,'restart_extend_velo',model%options%restart_extend_velo)
-    call GetValue(section,'forcewrite_restart',model%options%forcewrite_restart)
 
   end subroutine handle_options
 
@@ -1653,17 +1652,18 @@ subroutine print_options(model)
        call write_log('  Slightly cheated with how temperature is implemented.',GM_WARNING)
     end if
 
-    if (model%options%is_restart == RESTART_TRUE) then
+    if (model%options%is_restart == STANDARD_RESTART) then
        call write_log('Restarting model from a previous run')
        if (model%options%restart_extend_velo == RESTART_EXTEND_VELO_TRUE) then
           call write_log('Using extended velocity fields for restart')
        endif
+    elseif (model%options%is_restart == HYBRID_RESTART) then
+       call write_log('Hybrid restart from a previous run')
+       if (model%options%restart_extend_velo == RESTART_EXTEND_VELO_TRUE) then
+          call write_log('Using extended velocity fields for restart')
+       endif
     end if
 
-    if (model%options%forcewrite_restart) then
-       call write_log('Will write to output files on restart')
-    endif
-
     !HO options
 
     if (model%options%whichdycore /= DYCORE_GLIDE) then   ! glissade higher-order

libglide/glide_temp.F90

@@ -199,7 +199,7 @@ subroutine glide_init_temp(model)
     !TODO - Make sure cells in the Glide temperature halo are initialized to reasonable values
     !       (not unphys_val), e.g. if reading temps from input or restart file.
 
-    if (model%options%is_restart == RESTART_TRUE) then
+    if (model%options%is_restart == STANDARD_RESTART .or. model%options%is_restart == HYBRID_RESTART) then
 
        ! Temperature has already been initialized from a restart file. 
        ! (Temperature is always a restart variable.)
@@ -291,7 +291,7 @@ subroutine glide_init_temp(model)
 
     ! ====== Calculate initial value of flwa ==================
 
-    if (model%options%is_restart == RESTART_FALSE) then
+    if (model%options%is_restart == NO_RESTART) then
        call write_log("Calculating initial flwa from temp and thk fields")
 
        ! Calculate Glen's A --------------------------------------------------------   

libglide/glide_types.F90

@@ -214,8 +214,9 @@ module glide_types
   integer, parameter :: SIGMA_COMPUTE_EVEN = 3
   integer, parameter :: SIGMA_COMPUTE_PATTYN = 4
 
-  integer, parameter :: RESTART_FALSE = 0
-  integer, parameter :: RESTART_TRUE = 1
+  integer, parameter :: NO_RESTART = 0
+  integer, parameter :: STANDARD_RESTART = 1
+  integer, parameter :: HYBRID_RESTART = 2
 
   integer, parameter :: RESTART_EXTEND_VELO_FALSE = 0
   integer, parameter :: RESTART_EXTEND_VELO_TRUE = 1
@@ -753,12 +754,12 @@ module glide_types
     !> \item[4] compute Pattyn sigma coordinates
     !> \end{description}
 
-    !TODO - Make is_restart a logical variable?
     integer :: is_restart = 0
     !> if the run is a restart of a previous run
     !> \begin{description}
     !> \item[0] normal start-up (take init fields from .nc input file OR if absent, use default options)
     !> \item[1] restart model from previous run (do not calc. temp, rate factor, or vel)
+    !> \item[2] hybrid restart; use restart from previous run as the input file, and reset the time
     !> \end{description}
 
     integer :: restart_extend_velo = 0
@@ -769,9 +770,6 @@ module glide_types
     !>          (required if restart velocities are nonzero on global boundaries)
     !> \end{description}
 
-    logical :: forcewrite_restart = .false.
-    !> flag that indicates whether to force writing of output on restart
-
     ! This is a Glimmer serial option
     ! The parallel code enforces periodic EW and NS boundary conditions by default
     logical :: periodic_ew = .false.

libglimmer/glimmer_ncio.F90

@@ -84,7 +84,7 @@ subroutine openall_out(model,outfiles)
 
           call glimmer_nc_openappend(oc,model)
 
-       elseif (model%options%is_restart == RESTART_TRUE) then   ! reopen the file if it exists
+       elseif (model%options%is_restart == STANDARD_RESTART) then   ! reopen the file if it exists
 
           status = parallel_open(process_path(oc%nc%filename),NF90_WRITE,oc%nc%id)
 
@@ -100,6 +100,7 @@ subroutine openall_out(model,outfiles)
           endif
 
        else  ! assume the file does not exist; create it
+             ! Note: For hybrid restarts, the file is created at initialization
 
           call glimmer_nc_createfile(oc, model)
 
@@ -714,7 +715,7 @@ subroutine glimmer_nc_checkread(infile,model,time)
 
     implicit none
 
-    type(glimmer_nc_input), pointer :: infile  !> structure containg output netCDF descriptor
+    type(glimmer_nc_input), pointer :: infile  !> structure containing output netCDF descriptor
     type(glide_global_type) :: model    !> the model instance
     real(dp),optional :: time           !> Optional alternative time
 
@@ -749,21 +750,36 @@ subroutine glimmer_nc_checkread(infile,model,time)
 
           if (pos /= 0 .or. pos_cesm /= 0) then   ! get the start time based on the current time slice
 
-             restart_time = infile%times(infile%current_time)      ! years
-             model%numerics%tstart = restart_time
-             model%numerics%time = restart_time
-
-             if (infile%tstep_counts_read) then
-                model%numerics%tstep_count = infile%tstep_counts(infile%current_time)
-             else
-                ! BACKWARDS_COMPATIBILITY(wjs, 2017-05-17) Older files may not have
-                ! 'tstep_count', so compute it ourselves here. We don't want to use this
-                ! formulation in general because it is prone to roundoff errors.
-                model%numerics%tstep_count = nint(model%numerics%time/model%numerics%tinc)
-             end if
-
-             write(message,*) 'Restart: New tstart, tstep_count =', model%numerics%tstart, model%numerics%tstep_count
-             call write_log(message)
+             if (model%options%is_restart == STANDARD_RESTART) then
+
+                restart_time = infile%times(infile%current_time)      ! years
+                model%numerics%tstart = restart_time
+                model%numerics%time = restart_time
+
+                if (infile%tstep_counts_read) then
+                   model%numerics%tstep_count = infile%tstep_counts(infile%current_time)
+                else
+                   ! BACKWARDS_COMPATIBILITY(wjs, 2017-05-17) Older files may not have
+                   ! 'tstep_count', so compute it ourselves here. We don't want to use this
+                   ! formulation in general because it is prone to roundoff errors.
+                   model%numerics%tstep_count = nint(model%numerics%time/model%numerics%tinc)
+                end if
+
+                write(message,*) 'Standard restart: New tstart, tstep_count =', &
+                     model%numerics%tstart, model%numerics%tstep_count
+                call write_log(message)
+
+             elseif (model%options%is_restart == HYBRID_RESTART) then
+
+                ! Use tstart from the config file, not the time from the restart file
+                model%numerics%time = model%numerics%tstart  ! years
+                model%numerics%tstep_count = 0
+
+                write(message,*) 'Hybrid restart: New tstart, tstep_count =', &
+                     model%numerics%tstart, model%numerics%tstep_count
+                call write_log(message)
+
+             endif   ! is_restart
 
           endif  ! pos/=0 or pos_cesm/=0
 

libglimmer/glimmer_ncparams.F90

@@ -75,17 +75,25 @@ subroutine glimmer_nc_readparams(model,config)
     character(len=fname_length) :: restart_filename
     character(len=256) :: message
 
-    ! Note on restart files:
+    ! Notes on restart files:
     ! If a file is listed in the 'CF restart' section, then it is added to the glimmer_nc_output data structure
     !  and written at the specified frequency.
-    ! If model%options%is_restart = RESTART_TRUE, then the file listed in 'CF restart' (provided it exists) 
+    !
+    ! If model%options%is_restart = STANDARD_RESTART, then the file listed in 'CF restart' (provided it exists)
     !  is added to the glimmer_nc_input data structure, overriding any file listed in the 'CF input' section.
     !  The latest time slice will be read in.
-    ! Thus when restarting the model, it is only necessary to set restart = RESTART_TRUE (i.e, restart = 1)
-    !  in the config file; it is not necesssary to change filenames in 'CF input' or 'CF restart'.
-    ! At most one file should be listed in the 'CF restart' section, and it should contain the string 'restart'
-    ! If model%options%is_restart = RESTART_TRUE and there is no 'CF restart' section, then the model will restart
+    ! Thus when restarting the model, it is only necessary to set restart = 1 (i.e., STANDARD_RESTART)
+    !  in the config file; it is not necesssary to change the filenames in 'CF input' or 'CF restart'.
+    ! At most one file should be listed in the 'CF restart' section, and it should contain the string 'restart' or '.r.'
+    ! If model%options%is_restart = STANDARD_RESTART and there is no 'CF restart' section, then the model will restart
     !  from the file and time slice specified in the 'CF input' section. (This is the old Glimmer behavior.)
+    !
+    ! If model%options%is_restart = HYBRID_RESTART, then the file listed in 'CF input' is used to initialize the model.
+    ! This file should be a restart file from a previous run (e.g., a long ice-sheet spin-up),
+    !  which provides the initial ice state for the hybrid run.
+    ! The differences from STANDARD_RESTART (besides the config section where the filename is given) are
+    ! (1) tstep_count is set to 0, replacing the value in the CF input file.
+    ! (2) model%numerics%time is set to tstart from the config file, replacing the value in the CF input file.
 
     ! get config string
     call ConfigAsString(config,configstring)
@@ -135,7 +143,7 @@ subroutine glimmer_nc_readparams(model,config)
     end do
 
     ! set up restart input
-    if (model%options%is_restart == RESTART_TRUE) then
+    if (model%options%is_restart == STANDARD_RESTART) then
 
        ! If there is a 'CF restart' section, the model will restart from the file listed there (if it exists).
        ! Else the model will start from the input file in the 'CF input' section.
@@ -187,7 +195,7 @@ subroutine glimmer_nc_readparams(model,config)
 
        endif   ! associated(section)
 
-    endif  ! model%options%is_restart = RESTART_TRUE
+    endif  ! model%options%is_restart
 
     ! setup forcings
     call GetSection(config,section,'CF forcing')