Minor code changes to support the slab test

I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error.
This is done only when the user does *not* specify adaptive subcycling.
The clean abort allows the new slabStability script to keep going, launching a new run
with a shorter timestep.

In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option
(option 3 of whichflwa).
This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa
in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1.
We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid
breaking config files in test cases.

Note: This option was added by Matt Hoffman in 2014.  I am unaware of test cases that
use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them
by switching to whichflwa = 0.

In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold
for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2.
For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors.

In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort
for large CFL violations (advective CFL > 10).
Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled.
This prevents excessive subcycling for egregious CFL violations.
If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled,
then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1.

I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2.
With the default value, simple prognostic tests like the dome are not mass-conserving.
Not changing just yet because answers will change for the dome test.

libglide/glide_setup.F90

@@ -883,11 +883,10 @@ subroutine print_options(model)
          'advective-diffusive balance ',&
          'temp from external file     ' /)
 
-    character(len=*), dimension(0:3), parameter :: flow_law = (/ &
-         'const 1e-16 Pa^-n a^-1      ', &
+    character(len=*), dimension(0:2), parameter :: flow_law = (/ &
+         'uniform factor flwa         ', &
          'Paterson and Budd (T = -5 C)', &
-         'Paterson and Budd           ', &
-         'read flwa/flwastag from file' /)
+         'Paterson and Budd           ' /)
 
     !TODO - Rename slip_coeff to which_btrc?
     character(len=*), dimension(0:5), parameter :: slip_coeff = (/ &
@@ -2034,9 +2033,9 @@ subroutine handle_parameters(section, model)
     call GetValue(section,'pmp_offset',         model%temper%pmp_offset)
     call GetValue(section,'pmp_threshold',      model%temper%pmp_threshold)
     call GetValue(section,'geothermal',         model%paramets%geot)
-    !TODO - Change default_flwa to flwa_constant?  Would have to change config files.
     call GetValue(section,'flow_factor',        model%paramets%flow_enhancement_factor)
     call GetValue(section,'flow_factor_float',  model%paramets%flow_enhancement_factor_float)
+    !TODO - Change default_flwa to flwa_constant?  Would have to change config files.
     call GetValue(section,'default_flwa',       model%paramets%default_flwa)
     call GetValue(section,'efvs_constant',      model%paramets%efvs_constant)
     call GetValue(section,'effstrain_min',      model%paramets%effstrain_min)

libglide/glide_types.F90

@@ -104,7 +104,6 @@ module glide_types
   integer, parameter :: FLWA_CONST_FLWA = 0
   integer, parameter :: FLWA_PATERSON_BUDD_CONST_TEMP = 1
   integer, parameter :: FLWA_PATERSON_BUDD = 2
-  integer, parameter :: FLWA_INPUT = 3
 
   integer, parameter :: BTRC_ZERO = 0
   integer, parameter :: BTRC_CONSTANT = 1
@@ -470,7 +469,6 @@ module glide_types
     !> \item[1] \emph{Paterson and Budd} relationship, 
     !> with temperature set to $-5^{\circ}\mathrm{C}$ 
     !> \item[2] \emph{Paterson and Budd} relationship
-    !> \item[3] Read flwa/flwastag from file
     !> \end{description}
 
     integer :: whichbtrc = 0
@@ -2148,6 +2146,7 @@ module glide_types
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   !TODO - Move these parameters to types associated with a certain kind of physics
+  !TODO - Set default geot = 0, so that idealized tests by default have no mass loss
   type glide_paramets
     real(dp),dimension(5) :: bpar = (/ 0.2d0, 0.5d0, 0.0d0 ,1.0d-2, 1.0d0/)
     real(dp) :: btrac_const = 0.d0     ! m yr^{-1} Pa^{-1} (gets scaled during init)

libglissade/glissade.F90

@@ -1899,7 +1899,7 @@ subroutine glissade_thermal_solve(model, dt)
                                 model%temper%btemp_ground,                                    & ! deg C
                                 model%temper%btemp_float,                                     & ! deg C
                                 bmlt_ground_unscaled)                                           ! m/s
-                                     
+
     ! Update basal hydrology, if needed
     ! Note: glissade_calcbwat uses SI units
 
@@ -1977,6 +1977,7 @@ subroutine glissade_thickness_tracer_solve(model)
     use glissade_bmlt_float, only: verbose_bmlt_float
     use glissade_calving, only: verbose_calving
     use glissade_grid_operators, only: glissade_vertical_interpolate
+    use glide_stop, only: glide_finalise
 
     implicit none
 
@@ -2165,21 +2166,25 @@ subroutine glissade_thickness_tracer_solve(model)
                                model%geomderv%dusrfdew*thk0/len0, model%geomderv%dusrfdns*thk0/len0,           &
                                model%velocity%uvel * scyr * vel0, model%velocity%vvel * scyr * vel0,           &
                                model%numerics%dt_transport * tim0 / scyr,                                      &
+                               model%numerics%adaptive_cfl_threshold,                                          &
                                model%numerics%adv_cfl_dt,         model%numerics%diff_cfl_dt)
 
        ! Set the transport timestep.
        ! The timestep is model%numerics%dt by default, but optionally can be reduced for subcycling
 
+       !WHL - debug
+!      if (main_task) then
+!         print*, 'Checked advective CFL threshold'
+!         print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr
+!         print*, 'adv_cfl_dt    =', model%numerics%adv_cfl_dt
+!      endif
+
+       advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt
+
        if (model%numerics%adaptive_cfl_threshold > 0.0d0) then
 
-          !WHL - debug
-!          if (main_task) then
-!             print*, 'Check advective CFL threshold'
-!             print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr
-!             print*, 'adv_cfl_dt    =', model%numerics%adv_cfl_dt
-!          endif
+          ! subcycle the transport when advective_cfl exceeds the threshold
 
-          advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt
           if (advective_cfl > model%numerics%adaptive_cfl_threshold) then
 
              ! compute the number of subcycles
@@ -2192,14 +2197,29 @@ subroutine glissade_thickness_tracer_solve(model)
                 print*, 'Ratio =', advective_cfl / model%numerics%adaptive_cfl_threshold
                 print*, 'nsubcyc =', nsubcyc
              endif
+
           else
              nsubcyc = 1
           endif
           dt_transport = model%numerics%dt * tim0 / real(nsubcyc,dp)   ! convert to s
 
        else  ! no adaptive subcycling
-          nsubcyc = model%numerics%subcyc
-          dt_transport = model%numerics%dt_transport * tim0  ! convert to s
+
+          advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt
+
+          ! If advective_cfl exceeds 1.0, then abort cleanly.  Otherwise, set dt_transport and proceed.
+          ! Note: Usually, it would be enough to write a fatal abort message.
+          !       The call to glide_finalise was added to allow CISM to finish cleanly when running
+          !        a suite of automated stability tests, e.g. with the stabilitySlab.py script.
+          if (advective_cfl > 1.0d0) then
+             if (main_task) print*, 'advective CFL violation; call glide_finalise and exit cleanly'
+             call glide_finalise(model, crash=.true.)
+             stop
+          else
+             nsubcyc = model%numerics%subcyc
+             dt_transport = model%numerics%dt_transport * tim0  ! convert to s
+          endif
+
        endif
 
        !-------------------------------------------------------------------------

libglissade/glissade_therm.F90

@@ -1115,11 +1115,14 @@ subroutine glissade_therm_driver(whichtemp,                         &
              dissipcol(ew,ns) = dissipcol(ew,ns) * thck(ew,ns)*rhoi*shci
 
              ! Verify that the net input of energy into the column is equal to the change in
-             ! internal energy.  
+             ! internal energy.
 
              delta_e = (ucondflx(ew,ns) - lcondflx(ew,ns) + dissipcol(ew,ns)) * dttem
 
-             if (abs((efinal-einit-delta_e)/dttem) > 1.0d-8) then
+             ! Note: For very small dttem (e.g., 1.0d-6 year or less), this error can be triggered
+             !       by roundoff error.  In that case, the user may need to increase the threshold.
+             ! July 2021: Increased from 1.0d-8 to 1.0d-7 to allow smaller dttem.
+             if (abs((efinal-einit-delta_e)/dttem) > 1.0d-7) then
 
                 if (verbose_column) then
                    print*, 'Ice thickness:', thck(ew,ns)
@@ -2416,7 +2419,7 @@ subroutine glissade_flow_factor(whichflwa,               whichtemp,  &
 
     integer,                   intent(in)    :: whichflwa !> which method of calculating A
     integer,                   intent(in)    :: whichtemp !> which method of calculating temperature;
-                                                           !> include waterfrac in calculation if using enthalpy method
+                                                          !> include waterfrac in calculation if using enthalpy method
     real(dp),dimension(:),     intent(in)    :: stagsigma !> vertical coordinate at layer midpoints
     real(dp),dimension(:,:),   intent(in)    :: thck      !> ice thickness (m)
     real(dp),dimension(:,:,:), intent(in)    :: temp      !> 3D temperature field (deg C)
@@ -2490,17 +2493,16 @@ subroutine glissade_flow_factor(whichflwa,               whichtemp,  &
     endif
 
     ! Multiply the default rate factor by the enhancement factor if applicable
-    ! Note: Here, default_flwa is assumed to have units of Pa^{-n} s^{-1},
+    ! Note: Here, the input default_flwa is assumed to have units of Pa^{-n} s^{-1},
     !       whereas model%paramets%default_flwa has units of Pa^{-n} yr^{-1}.
 
     ! initialize
-    if (whichflwa /= FLWA_INPUT) then
-       do ns = 1, nsn
-          do ew = 1, ewn
-             flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa
-          enddo
+    !TODO - Move the next few lines inside the select case construct.
+    do ns = 1, nsn
+       do ew = 1, ewn
+          flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa
        enddo
-    endif
+    enddo
 
     select case(whichflwa)
 
@@ -2560,21 +2562,12 @@ subroutine glissade_flow_factor(whichflwa,               whichtemp,  &
       end do
 
     case(FLWA_CONST_FLWA)
-
-       ! do nothing (flwa is initialized to default_flwa above)
-
-    case(FLWA_INPUT)
-      ! do nothing - use flwa from input or forcing file
-      print *, 'FLWA', minval(flwa), maxval(flwa)
+       ! do nothing (flwa is set above, with units Pa^{-n} s^{-1})
 
     end select
 
-    ! This logic assumes that the input flwa is already in dimensionless model units.
-    ! TODO: Make a different assumption about input units?
-    if (whichflwa /= FLWA_INPUT) then
-       ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1})
-       flwa(:,:,:) = flwa(:,:,:) / vis0
-    endif
+    ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1})
+    flwa(:,:,:) = flwa(:,:,:) / vis0
 
     deallocate(enhancement_factor)
 

libglissade/glissade_transport.F90

@@ -979,6 +979,7 @@ subroutine glissade_check_cfl(ewn,     nsn,      nlyr,          &
                                   parallel,                         &
                                   stagthk, dusrfdew, dusrfdns,      &
                                   uvel,    vvel,     deltat,        &
+                                  adaptive_cfl_threshold,           &
                                   allowable_dt_adv,  allowable_dt_diff)
 
       ! Calculate maximum allowable time step based on both 
@@ -1015,6 +1016,10 @@ subroutine glissade_check_cfl(ewn,     nsn,      nlyr,          &
       real(dp), intent(in) :: &
          deltat      ! model deltat (yrs)
 
+      real(dp), intent(in) :: &
+         adaptive_cfl_threshold  ! threshold for adaptive subcycling
+                                 ! if = 0, there is no adaptive subcycling; code aborts when CFL > 1
+
       real(dp), intent(out) :: &
          allowable_dt_adv     ! maximum allowable dt (yrs) based on advective CFL 
 
@@ -1159,13 +1164,16 @@ subroutine glissade_check_cfl(ewn,     nsn,      nlyr,          &
           write(message,*) 'Advective CFL violation!  Maximum allowable time step for advective CFL condition is ' &
                // trim(adjustl(dt_string)) // ' yr, limited by global position i=' &
                // trim(adjustl(xpos_string)) // ' j=' //trim(adjustl(ypos_string))
+          call write_log(trim(message),GM_WARNING)
 
-          ! If the violation is egregious (defined as deltat > 10 * allowable_dt_adv), then abort.
-          ! Otherwise, write a warning and proceed.
-          if (deltat > 10.d0 * allowable_dt_adv) then
-             call write_log(trim(message),GM_FATAL)
-          else
-             call write_log(trim(message),GM_WARNING)
+          ! If adaptive subcyling is allowed, then make the code abort for egregious CFL violations,
+          ! (defined as deltat > 10 * allowable_dt_adv), to prevent excessive subcycling.
+
+          if (main_task .and. adaptive_cfl_threshold > 0.0d0) then
+             if (deltat > 10.d0 * allowable_dt_adv) then
+                print*, 'deltat, allowable_dt_adv, ratio =', deltat, allowable_dt_adv, deltat/allowable_dt_adv
+                call write_log('Aborting with CFL violation', GM_FATAL)
+             endif
           endif
 
       endif

libglissade/glissade_velo_higher.F90

@@ -200,8 +200,8 @@ module glissade_velo_higher
 !    logical :: verbose = .true.  
     logical :: verbose_init = .false.   
 !    logical :: verbose_init = .true.   
-!    logical :: verbose_solver = .false.
-    logical :: verbose_solver = .true.
+    logical :: verbose_solver = .false.
+!    logical :: verbose_solver = .true.
     logical :: verbose_Jac = .false.
 !    logical :: verbose_Jac = .true.
     logical :: verbose_residual = .false.