Updated the Dukowicz slab test case

libglide/felix_dycore_interface.F90

@@ -146,7 +146,7 @@ end subroutine felix_velo_init
    subroutine felix_velo_driver(model)
 
       use glimmer_global, only : dp
-      use glimmer_physcon, only: gn, scyr
+      use glimmer_physcon, only: scyr
       use glimmer_paramets, only: thk0, len0, vel0, vis0
       use glimmer_log
       use glide_types

libglide/glide_setup.F90

@@ -160,7 +160,7 @@ end subroutine glide_printconfig
   subroutine glide_scale_params(model)
     !> scale parameters
     use glide_types
-    use glimmer_physcon,  only: scyr, gn
+    use glimmer_physcon,  only: scyr
     use glimmer_paramets, only: thk0, tim0, len0, vel0, vis0, acc0, tau0
 
     implicit none
@@ -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 = (/ &
@@ -1996,7 +1995,7 @@ subroutine handle_parameters(section, model)
     use glimmer_config
     use glide_types
     use glimmer_log
-    use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt
+    use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt, n_glen
 
     implicit none
     type(ConfigSection), pointer :: section
@@ -2021,6 +2020,7 @@ subroutine handle_parameters(section, model)
     call GetValue(section,'lhci', lhci)
     call GetValue(section,'trpt', trpt)
 #endif
+    call GetValue(section,'n_glen', n_glen)
 
     loglevel = GM_levels-GM_ERROR
     call GetValue(section,'log_level',loglevel)
@@ -2033,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)
@@ -2206,7 +2206,7 @@ end subroutine handle_parameters
 
   subroutine print_parameters(model)
 
-    use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav
+    use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav, n_glen
     use glide_types
     use glimmer_log
     implicit none
@@ -2371,6 +2371,9 @@ subroutine print_parameters(model)
     write(message,*) 'triple point of water (K)     : ', trpt
     call write_log(message)
 
+    write(message,*) 'Glen flow law exponent        : ', n_glen
+    call write_log(message)
+
     write(message,*) 'geothermal flux  (W/m^2)      : ', model%paramets%geot
     call write_log(message)
 

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)

libglimmer/glimmer_paramets.F90

@@ -33,7 +33,7 @@
 module glimmer_paramets
 
   use glimmer_global, only : dp
-  use glimmer_physcon, only : scyr, gn
+  use glimmer_physcon, only : scyr
 
   implicit none
   save
@@ -118,6 +118,7 @@ module glimmer_paramets
   real(dp), parameter :: grav_glam = 9.81d0           ! m s^{-2}
 
   ! GLAM scaling parameters; units are correct if thk0 has units of meters
+  integer, parameter :: gn = 3                              ! Glen flow exponent; fixed at 3 for purposes of setting vis0
   real(dp), parameter :: tau0 = rhoi_glam*grav_glam*thk0    ! stress scale in GLAM ( Pa )  
   real(dp), parameter :: evs0 = tau0 / (vel0/len0)          ! eff. visc. scale in GLAM ( Pa s )
   real(dp), parameter :: vis0 = tau0**(-gn) * (vel0/len0)   ! rate factor scale in GLAM ( Pa^-3 s^-1 )

libglimmer/glimmer_physcon.F90

@@ -79,11 +79,17 @@ module glimmer_physcon
 !  real(dp) :: trpt = 273.16d0                     !< Triple point of water (K)
 #endif
 
+  ! The Glen flow-law exponent, n_glen, is used in Glissade.
+  ! It is not a parameter, because the default can be overridden in the config file.
+  ! TODO: Allow setting n_glen independently for each ice sheet instance?
+  ! Note: Earlier code used an integer parameter, gn = 3, for all flow-law calculations.
+  !       For backward compatiblity, gn = 3 is retained for Glide.
+  real(dp) :: n_glen = 3.0d0                     !< Exponent in Glen's flow law; user-configurable real(dp) in Glissade
+  integer, parameter :: gn = 3                   !< Exponent in Glen's flow law; fixed integer parameter in Glide
   real(dp),parameter :: celsius_to_kelvin = 273.15d0  !< Note: Not quite equal to trpt
   real(dp),parameter :: scyr = 31536000.d0       !< Number of seconds in a year of exactly 365 days
   real(dp),parameter :: rhom = 3300.0d0          !< The density of magma(?) (kg m<SUP>-3</SUP>) 
   real(dp),parameter :: rhos = 2600.0d0          !< The density of solid till (kg m$^{-3}$) 
-  integer, parameter :: gn = 3                   !< The power dependency of Glen's flow law.
   real(dp),parameter :: actenh = 139.0d3         !< Activation energy in Glen's flow law for \f$T^{*}\geq263\f$K. (J mol<SUP>-1</SUP>)
   real(dp),parameter :: actenl = 60.0d3          !< Activation energy in Glen's flow law for \f$T^{*}<263\f$K. (J mol<SUP>-1</SUP>)  
   real(dp),parameter :: arrmlh = 1.733d3         !< Constant of proportionality in Arrhenius relation

libglimmer/glimmer_scales.F90

@@ -45,7 +45,7 @@ subroutine glimmer_init_scales
 
     ! set scale factors for I/O (can't have non-integer powers)
 
-    use glimmer_physcon, only : scyr, gn
+    use glimmer_physcon, only : scyr
     use glimmer_paramets, only : thk0, tim0, vel0, vis0, len0, acc0, tau0, evs0
     implicit none
 

libglissade/glissade.F90

@@ -457,7 +457,11 @@ subroutine glissade_initialise(model, evolve_ice)
     ! handle relaxed/equilibrium topo
     ! Initialise isostasy first
 
-    call init_isostasy(model)
+    if (model%options%isostasy == ISOSTASY_COMPUTE) then
+
+       call init_isostasy(model)
+
+    endif
 
     select case(model%isostasy%whichrelaxed)
 
@@ -1895,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
 
@@ -1973,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
 
@@ -2161,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
@@ -2188,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_basal_traction.F90

@@ -440,9 +440,12 @@ subroutine calcbeta (whichbabc,                    &
        !       If this factor is not present in the input file, it is set to 1 everywhere.
 
        ! Compute beta
-       ! gn = Glen's n from physcon module
-       beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/gn - 1.0d0) * &
-            (speed(:,:) + basal_physics%effecpress_stag(:,:)**gn * big_lambda)**(-1.0d0/gn)
+       ! Note: Where this equation has powerlaw_m, we used to have Glen's flow exponent n,
+       !       following the notation of Leguy et al. (2014).
+       !       Changed to powerlaw_m to be consistent with the Schoof and Tsai laws.
+       m = basal_physics%powerlaw_m
+       beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/m - 1.0d0) * &
+            (speed(:,:) + basal_physics%effecpress_stag(:,:)**m * big_lambda)**(-1.0d0/m)
 
        ! If c_space_factor /= 1.0 everywhere, then multiply beta by c_space_factor
        if (maxval(abs(basal_physics%c_space_factor_stag(:,:) - 1.0d0)) > tiny(0.0d0)) then

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)
@@ -1640,10 +1643,11 @@ subroutine glissade_enthalpy_matrix_elements(dttem,                       &
     ! At each temperature point, compute the temperature part of the enthalpy.
     ! enth_T = enth for cold ice, enth_T < enth for temperate ice
 
-    enth_T(0) = rhoi*shci*temp(0)  !WHL - not sure enth_T(0) is needed
-    do up = 1, upn
+    do up = 1, upn-1
        enth_T(up) = (1.d0 - waterfrac(up)) * rhoi*shci*temp(up)
     enddo
+    enth_T(0) = rhoi*shci*temp(0)
+    enth_T(up) = rhoi*shci*temp(up)
 
 !WHL - debug
     if (verbose_column) then
@@ -2250,8 +2254,7 @@ subroutine glissade_interior_dissipation_sia(ewn,       nsn,       &
     ! Compute the dissipation source term associated with strain heating,
     ! based on the shallow-ice approximation.
 
-    use glimmer_physcon, only : gn   ! Glen's n
-    use glimmer_physcon, only: rhoi, shci, grav
+    use glimmer_physcon, only: rhoi, shci, grav, n_glen
 
     integer, intent(in) :: ewn, nsn, upn   ! grid dimensions
 
@@ -2267,12 +2270,14 @@ subroutine glissade_interior_dissipation_sia(ewn,       nsn,       &
     real(dp), dimension(:,:,:), intent(out) ::  &
          dissip       ! interior heat dissipation (deg/s)
 
-    integer, parameter :: p1 = gn + 1  
-
     integer :: ew, ns
     real(dp), dimension(upn-1) :: sia_dissip_fact  ! factor in SIA dissipation calculation
     real(dp) :: geom_fact         ! geometric factor
 
+    real(dp) :: p1                ! exponent = n_glen + 1
+
+    p1 = n_glen + 1.0d0
+
     ! Two methods of doing this calculation: 
     ! 1. find dissipation at u-pts and then average
     ! 2. find dissipation at H-pts by averaging quantities from u-pts
@@ -2414,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)
@@ -2488,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)
 
@@ -2558,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)