diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 909244fb..c5f4cc3d 100644 --- a/libglide/glide_setup.F90 +++ b/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) diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index 147608e9..d55a727e 100644 --- a/libglide/glide_types.F90 +++ b/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) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 4093a137..efa2d483 100644 --- a/libglissade/glissade.F90 +++ b/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 !------------------------------------------------------------------------- diff --git a/libglissade/glissade_therm.F90 b/libglissade/glissade_therm.F90 index f6364650..f155521b 100644 --- a/libglissade/glissade_therm.F90 +++ b/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) diff --git a/libglissade/glissade_transport.F90 b/libglissade/glissade_transport.F90 index a1c57219..e0974b96 100644 --- a/libglissade/glissade_transport.F90 +++ b/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 diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index e9845437..77b2f9f5 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/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.