MODULE stopping_crit_core !----------------------------------------------------------------------- ! NAME ! stopping_crit_core ! ! DESCRIPTION ! Stopping criteria for PEM. ! ! AUTHORS & DATE ! JB Clement, 09/2023 ! ! NOTES ! The surface ice is monitored by two criteria: ! - the disappearance of the patch that was sublimating at ! initialization; ! - the appearance of ice outside of what the PCM announced, that ! is to say outside of the cells that already held ice or that ! were condensing at initialization. Such reservoirs can only be ! built by a process internal to the PEM, like the glacier flow. ! For H2O, only the ice thick enough to be a perennial reservoir ! for the PCM counts as a new reservoir. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4 ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- real(dp), private :: tol_h2oice_surf ! Tolerance for the relative change of the H2O ice surface before stopping the PEM real(dp), private :: tol_co2ice_surf ! Tolerance for the relative change of the CO2 ice surface before stopping the PEM real(dp), private :: tol_ice_planetfrac ! Accepted change of a surface ice coverage, as a fraction of the total planet surface real(dp), private :: tol_ps_avg ! Tolerance for the relative change of averaged surface pressure before stopping the PEM ! Minimum change of ice-covered surface, in number of ice-bearing cells, before stopping the PEM (see the function "eff_tol_area" below) real(dp), parameter, private :: tol_ice_ncell = 0._dp ! Default disables the mechanism, = 1. neutralizes a single cell switch ! Upper bound of the per-cell floor of the accepted change of an ice surface (see the function "eff_tol_area" below) real(dp), parameter, private :: tol_coverage_max = 0.5_dp ! Terms setting the accepted change of an ice surface (see the function "eff_tol_area" below) integer(di), parameter, private :: TOL_FROM_USER = 1 ! The relative tolerance asked by the user integer(di), parameter, private :: TOL_FROM_PLANET = 2 ! The fraction of the total planet surface integer(di), parameter, private :: TOL_FROM_GRID = 3 ! The per-cell floor integer(di), parameter, private :: TOL_FROM_CAP = 4 ! The per-cell floor limited by its upper bound ! Stopping codes integer(di), parameter, private :: CODE_H2OICE_SURF = 1 integer(di), parameter, private :: CODE_H2O_FLUX = 2 integer(di), parameter, private :: CODE_CO2ICE_SURF = 3 integer(di), parameter, private :: CODE_PRESSURE = 4 integer(di), parameter, private :: CODE_NMAX_YR_RUN = 5 integer(di), parameter, private :: CODE_NMAX_YR_RUNORB = 6 integer(di), parameter, private :: CODE_NMAX_YR_SIM = 7 integer(di), parameter, private :: CODE_TIME_LIMIT = 8 integer(di), parameter, private :: CODE_MASS_BALANCE = 9 integer(di), parameter, private :: CODE_H2OICE_NEWRESVR = 10 integer(di), parameter, private :: CODE_CO2ICE_NEWRESVR = 11 ! DATA STRUCTURES ! --------------- type :: stopFlags logical(k4) :: surface_h2oice_change = .false. logical(k4) :: newresvr_h2oice_change = .false. logical(k4) :: h2o_flux_balance_unclosed = .false. logical(k4) :: surface_co2ice_change = .false. logical(k4) :: newresvr_co2ice_change = .false. logical(k4) :: pressure_change = .false. logical(k4) :: nmax_yr_run_reached = .false. logical(k4) :: nmax_yr_runorb_reached = .false. logical(k4) :: nmax_yr_sim_reached = .false. logical(k4) :: time_limit_reached = .false. logical(k4) :: mass_balance_exceeded = .false. contains procedure :: is_any_set procedure :: stop_code procedure :: stop_message end type ! VISIBILITY ! ---------- private :: is_any_set, stop_code, stop_message private :: eff_tol_area, report_tol_source contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE set_stopping_crit_config(tol_h2oice_surf_in,tol_co2ice_surf_in,tol_ice_planetfrac_in,tol_ps_avg_in) !----------------------------------------------------------------------- ! NAME ! set_stopping_crit_config ! ! DESCRIPTION ! Setter for "stopping_crit_core" configuration parameters. ! ! AUTHORS & DATE ! JB Clement, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: real2str use display, only: print_msg, LVL_NFO use stoppage, only: stop_clean ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: tol_h2oice_surf_in, tol_co2ice_surf_in, tol_ice_planetfrac_in, tol_ps_avg_in ! CODE ! ---- tol_h2oice_surf = tol_h2oice_surf_in tol_co2ice_surf = tol_co2ice_surf_in tol_ice_planetfrac = tol_ice_planetfrac_in tol_ps_avg = tol_ps_avg_in call print_msg('tol_h2oice_surf = '//real2str(tol_h2oice_surf),LVL_NFO) call print_msg('tol_co2ice_surf = '//real2str(tol_co2ice_surf),LVL_NFO) call print_msg('tol_ice_planetfrac = '//real2str(tol_ice_planetfrac),LVL_NFO) call print_msg('tol_ps_avg = '//real2str(tol_ps_avg),LVL_NFO) if (tol_h2oice_surf <= 0._dp .or. tol_h2oice_surf > 1._dp) call stop_clean(__FILE__,__LINE__,'''tol_h2oice_surf'' outside admissible range ]0.,1.]!',1) if (tol_co2ice_surf <= 0._dp .or. tol_co2ice_surf > 1._dp) call stop_clean(__FILE__,__LINE__,'''tol_co2ice_surf'' outside admissible range ]0.,1.]!',1) if (tol_ice_planetfrac < 0._dp .or. tol_ice_planetfrac > 1._dp) call stop_clean(__FILE__,__LINE__,'''tol_ice_planetfrac'' outside admissible range [0.,1.]!',1) if (tol_ps_avg <= 0._dp .or. tol_ps_avg > 1._dp) call stop_clean(__FILE__,__LINE__,'''tol_ps_avg'' outside admissible range ]0.,1.]!',1) END SUBROUTINE set_stopping_crit_config !======================================================================= !======================================================================= FUNCTION is_any_set(this) RESULT(stopPEM) !----------------------------------------------------------------------- ! NAME ! is_any_set ! ! DESCRIPTION ! Detect if any stopping flag is set to true. ! ! AUTHORS & DATE ! JB Clement, 11/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- class(stopFlags), intent(in) :: this ! Stopping flags object ! RESULT ! ------ logical(k4) :: stopPEM ! True if any flag is set ! CODE ! ---- stopPEM = this%surface_h2oice_change .or. & this%newresvr_h2oice_change .or. & this%surface_co2ice_change .or. & this%newresvr_co2ice_change .or. & this%pressure_change .or. & this%h2o_flux_balance_unclosed .or. & this%nmax_yr_run_reached .or. & this%nmax_yr_runorb_reached .or. & this%nmax_yr_sim_reached .or. & this%time_limit_reached .or. & this%mass_balance_exceeded END FUNCTION is_any_set !======================================================================= !======================================================================= FUNCTION stop_code(this) result(code) !----------------------------------------------------------------------- ! NAME ! stop_code ! ! DESCRIPTION ! Return digit code corresponding to the active stopping flag. ! ! AUTHORS & DATE ! JB Clement, 11/2025 ! ! NOTES ! The physical criteria are tested before the bookkeeping ones so ! that a run whose climate has changed too much is reported as a ! climate stop even if it reaches a year or time limit. !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- class(StopFlags), intent(in) :: this ! Stopping flags object ! RESULT ! ------ integer(di) :: code ! Code corresponding to active flag (0 if none) ! CODE ! ---- code = 0 if (this%surface_h2oice_change) then; code = CODE_H2OICE_SURF else if (this%newresvr_h2oice_change) then; code = CODE_H2OICE_NEWRESVR else if (this%h2o_flux_balance_unclosed) then; code = CODE_H2O_FLUX else if (this%surface_co2ice_change) then; code = CODE_CO2ICE_SURF else if (this%newresvr_co2ice_change) then; code = CODE_CO2ICE_NEWRESVR else if (this%pressure_change) then; code = CODE_PRESSURE else if (this%nmax_yr_run_reached) then; code = CODE_NMAX_YR_RUN else if (this%nmax_yr_runorb_reached) then; code = CODE_NMAX_YR_RUNORB else if (this%nmax_yr_sim_reached) then; code = CODE_NMAX_YR_SIM else if (this%time_limit_reached) then; code = CODE_TIME_LIMIT else if (this%mass_balance_exceeded) then; code = CODE_MASS_BALANCE end if END FUNCTION stop_code !======================================================================= !======================================================================= FUNCTION stop_message(this) result(msg) !----------------------------------------------------------------------- ! NAME ! stop_message ! ! DESCRIPTION ! Return descriptive message corresponding to the active stopping flag. ! ! AUTHORS & DATE ! JB Clement, 11/2025 ! ! NOTES ! The same order as in "stop_code" must be respected. !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- class(StopFlags), intent(in) :: this ! Stopping flags object ! RESULT ! ------ character(:), allocatable :: msg ! Descriptive message for active flag ! CODE ! ---- if (this%surface_h2oice_change) then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)." else if (this%newresvr_h2oice_change) then; msg = "**** STOPPING because new H2O ice reservoirs have appeared (see message above)." else if (this%h2o_flux_balance_unclosed) then; msg = "**** STOPPING because H2O flux imbalance cannot be closed under physical bounds (see message above)." else if (this%surface_co2ice_change) then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)." else if (this%newresvr_co2ice_change) then; msg = "**** STOPPING because new CO2 ice reservoirs have appeared (see message above)." else if (this%pressure_change) then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)." else if (this%nmax_yr_run_reached) then; msg = "**** STOPPING because maximum number of years is reached." else if (this%nmax_yr_runorb_reached) then; msg = "**** STOPPING because maximum number of years due to orbital parameters is reached." else if (this%nmax_yr_sim_reached) then; msg = "**** STOPPING because the number of years to be simulated is reached." else if (this%time_limit_reached) then; msg = "**** STOPPING because the job time limit is going to be reached." else if (this%mass_balance_exceeded) then; msg = "**** STOPPING because a mass balance is not conserved (see message above)." else msg = "**** STOPPING for unknown reason (this should not happen!)." end if END FUNCTION stop_message !======================================================================= !======================================================================= FUNCTION eff_tol_area(tol,cell_area_ini,coverage_ref,tol_planetfrac,win_tol) RESULT(tol_area) !----------------------------------------------------------------------- ! NAME ! eff_tol_area ! ! DESCRIPTION ! Give the change of an ice-covered area that is accepted before stopping the PEM. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! The accepted change is the loosest of three terms: ! - the relative tolerance asked by the user, worth a fraction of ! the reference area; ! - the absolute allowance 'tol_planetfrac', worth a fraction of ! the total planet surface, which gives the criterion a meaning ! that does not depend on the size of the ice reservoir. ! Be carefull: if too loose, the initial patch can disappear ! entirely without stopping the PEM; ! - the granularity that the grid can resolve, which requires the ! the change to be worth at least 'tol_ice_ncell' cells. It ! prevents the PEM to stop when a single cell is switching on/off ! with the ice spanning only a few cells on a coarse grid. That ! per-cell field is floor shrinks as the grid refines and is inert ! wherever the ice well resolved. ! Putting 'tol_ice_ncell = 0.' disables it entirely. ! This per-cell floor is safeguarded from above by 'tol_coverage_max': ! an initial ice patch spanning fewer cells than 'tol_ice_ncell' ! would ask for a change larger than itself. So even its complete ! disappearance could not meet the stopping criterion. ! ! The last two terms are deactivated by default ('tol_ice_ncell = 0.' ! and 'tol_ice_planetfrac = 0.'). !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: total_surface ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: tol ! Relative tolerance asked by the user [-] real(dp), intent(in) :: cell_area_ini ! Largest cell contribution to the reference area [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] real(dp), intent(in), optional :: tol_planetfrac ! Accepted change as a fraction of the total planet surface [-] integer(di), intent(out), optional :: win_tol ! Tolerance setting the accepted change ! RESULT ! ------ real(dp) :: tol_area ! Accepted change [m2] ! LOCAL VARIABLES ! --------------- real(dp) :: area_user ! Change worth the relative tolerance asked by the user [m2] real(dp) :: area_planet ! Change worth the accepted fraction of the total planet surface [m2] real(dp) :: area_grid ! Change worth 'tol_ice_ncell' cells, bounded from above [m2] ! CODE ! ---- area_user = tol*coverage_ref area_planet = 0._dp if (present(tol_planetfrac)) area_planet = tol_planetfrac*total_surface area_grid = min(tol_ice_ncell*cell_area_ini,tol_coverage_max*coverage_ref) tol_area = max(area_user,area_planet,area_grid) ! Identification of the term which sets the accepted change if (present(win_tol)) then if (tol_area <= area_user) then win_tol = TOL_FROM_USER else if (tol_area <= area_planet) then win_tol = TOL_FROM_PLANET else if (tol_ice_ncell*cell_area_ini > tol_coverage_max*coverage_ref) then win_tol = TOL_FROM_CAP else win_tol = TOL_FROM_GRID end if end if END FUNCTION eff_tol_area !======================================================================= !======================================================================= SUBROUTINE report_tol_source(win_tol,tol,name_msg) !----------------------------------------------------------------------- ! NAME ! report_tol_source ! ! DESCRIPTION ! Tell which term sets the change of an ice-covered area that is ! accepted before stopping the PEM. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid use display, only: print_msg, LVL_NFO use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer(di), intent(in) :: win_tol ! Tolerance setting the accepted change real(dp), intent(in) :: tol ! Relative tolerance asked by the user [-] character(*), intent(in) :: name_msg ! Name of the accepted change for the messages ! CODE ! ---- select case (win_tol) case (TOL_FROM_PLANET) call print_msg("The accepted "//name_msg//" is set by 'tol_ice_planetfrac' = "//real2str(tol_ice_planetfrac*100._dp)// & " % of the planet surface instead of the tolerance of "//real2str(tol*100._dp)//" %.",LVL_NFO) case (TOL_FROM_GRID) call print_msg("The accepted "//name_msg//" is set by the grid resolution through 'tol_ice_ncell' = "//real2str(tol_ice_ncell)// & " instead of the tolerance of "//real2str(tol*100._dp)//" %.",LVL_NFO) case (TOL_FROM_CAP) if (ngrid == 1) then ! 1D call print_msg("Surface ice lies on a single column: the criterion reduces to the presence of ice in that column.",LVL_NFO) else ! 3D call print_msg("The initial surface ice spans fewer cells than 'tol_ice_ncell': the accepted "//name_msg// & " is capped instead of following the grid resolution.",LVL_NFO) end if case default ! Nothing to say when this is the relative tolerance asked by the user end select END SUBROUTINE report_tol_source !======================================================================= !======================================================================= SUBROUTINE check_area_criterion(name_msg,tol,area_ini,area_now,coverage_ref,cell_area_ref,flag) !----------------------------------------------------------------------- ! NAME ! check_area_criterion ! ! DESCRIPTION ! Compare the change of a monitored area to what is accepted and raise ! the stopping flag if it is exceeded. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! The monitored area and the reference area are two different things: ! a criterion can watch a sub-patch (the initially sublimating ice, ! the frozen part of the lakes, etc) while measuring its change ! against the whole patch of the reservoirs holding the species. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: total_surface use display, only: print_msg, LVL_NFO use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- character(*), intent(in) :: name_msg ! Name of the monitored area for the messages real(dp), intent(in) :: tol ! Relative tolerance asked by the user [-] real(dp), intent(in) :: area_ini ! Monitored area at initialization [m2] real(dp), intent(in) :: area_now ! Monitored area now [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] real(dp), intent(in) :: cell_area_ref ! Largest cell contribution to the reference area [m2] logical(k4), intent(inout) :: flag ! Stopping flag raised if the criterion is met ! LOCAL VARIABLES ! --------------- real(dp) :: tol_area ! Accepted change [m2] integer(di) :: win_tol ! Tolerance setting the accepted change ! CODE ! ---- if (coverage_ref <= 0._dp) return tol_area = eff_tol_area(tol,cell_area_ref,coverage_ref,tol_ice_planetfrac,win_tol) ! Check of the criterion if (abs(area_now - area_ini) > tol_area) then flag = .true. call print_msg("Monitored "//name_msg//" (ini|now) [m2] = "//real2str(area_ini)//' | '//real2str(area_now),LVL_NFO) call print_msg("Reference area of the criterion [m2] = "//real2str(coverage_ref),LVL_NFO) call print_msg("Change (accepted|current) [m2] = +/- "//real2str(tol_area)//' | '//real2str(area_now - area_ini),LVL_NFO) call print_msg("Change (accepted|current) [% of the reference area] = +/- "//real2str(100._dp*tol_area/coverage_ref)//' | '//real2str(100._dp*(area_now - area_ini)/coverage_ref),LVL_NFO) call print_msg("Change (accepted|current) [% of the planet surface] = +/- "//real2str(100._dp*tol_area/total_surface)//' | '//real2str(100._dp*(area_now - area_ini)/total_surface),LVL_NFO) call report_tol_source(win_tol,tol,'change') end if END SUBROUTINE check_area_criterion !======================================================================= !======================================================================= SUBROUTINE report_area_criterion(name_msg,tol,area_ini,coverage_ref,cell_area_ref) !----------------------------------------------------------------------- ! NAME ! report_area_criterion ! ! DESCRIPTION ! Report the monitored area at initialization, its reference area and ! the change that will be accepted. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: total_surface use display, only: print_msg, LVL_NFO, LVL_WRN use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- character(*), intent(in) :: name_msg ! Name of the monitored area for the messages real(dp), intent(in) :: tol ! Relative tolerance asked by the user [-] real(dp), intent(in) :: area_ini ! Monitored area at initialization [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] real(dp), intent(in) :: cell_area_ref ! Largest cell contribution to the reference area [m2] ! LOCAL VARIABLES ! --------------- real(dp) :: tol_area ! Accepted change [m2] integer(di) :: win_tol ! Tolerance setting the accepted change ! CODE ! ---- call print_msg("# Stopping criterion for "//name_msg//":",LVL_NFO) call print_msg("### Initial monitored area [m2] = "//real2str(area_ini),LVL_NFO) if (coverage_ref <= 0._dp) then call print_msg("The reference area is null so the stopping criterion on the "//name_msg//" will not be checked!",LVL_WRN) return end if tol_area = eff_tol_area(tol,cell_area_ref,coverage_ref,tol_ice_planetfrac,win_tol) call print_msg("### Reference area [m2] = "//real2str(coverage_ref),LVL_NFO) call print_msg("### Accepted change [m2] = +/- "//real2str(tol_area),LVL_NFO) call print_msg("### Accepted change [% of the reference area | % of planet surface] = +/- "//real2str(100._dp*tol_area/coverage_ref)//' | '//real2str(100._dp*tol_area/total_surface),LVL_NFO) if (area_ini > 0._dp .and. tol_area >= area_ini) call print_msg("The accepted change exceeds that surface. Even its complete disappearance will not meet the stopping criterion!",LVL_WRN) call report_tol_source(win_tol,tol,'change') END SUBROUTINE report_area_criterion !======================================================================= !======================================================================= SUBROUTINE stopping_crit_h2oice(h2oice_sublim_coverage_ini,is_h2oice_sublim_ini,is_h2oice_ini,coverage_ref,h2o_ice,h2oice_resvr_min,stopcrit) !----------------------------------------------------------------------- ! NAME ! stopping_crit_h2oice ! ! DESCRIPTION ! Check if the H2O ice which was sublimating has disappeared too much ! or if new H2O ice reservoirs have appeared. ! ! AUTHORS & DATE ! R. Vandemeulebrouck, 07/2022 ! L. Lange, 09/2022 ! JB Clement, 09/2023 ! ! NOTES ! Unlike its CO2 counterpart, this criterion is also checked in 1D, ! where it degenerates into a test of the presence of H2O ice in the ! single column of the grid. ! It is kept because it is then the only monitoring left on the H2O ! reservoir in 1D. Without it, the PEM would keep extrapolating the ! tendencies while ice has already vanished. ! A new reservoir is counted only where the ice exceeds ! 'h2oice_resvr_min'. A thinner ice is handed to the PCM as frost and ! removed from the reservoirs. So it would be counted again as new ! at each PEM run. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use slopes, only: masked_area, subslope_area_max ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), dimension(:,:), intent(in) :: is_h2oice_sublim_ini ! Initial location of the sublimating H2O ice logical(k4), dimension(:,:), intent(in) :: is_h2oice_ini ! Initial location of the H2O ice reservoirs (present or announced condensing) real(dp), dimension(:,:), intent(in) :: h2o_ice ! Current H2O surface ice [kg/m2] real(dp), intent(in) :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] real(dp), intent(in) :: h2oice_resvr_min ! Minimum H2O ice to be a perennial reservoir [kg/m2] class(stopFlags), intent(inout) :: stopcrit ! CODE ! ---- ! Check to escape the subroutine (if the PEM should stop already) if (stopcrit%stop_code() > 0) return ! Disappearance of the patch which was sublimating at initialization call check_area_criterion('sublimating H2O ice surface',tol_h2oice_surf,h2oice_sublim_coverage_ini, & masked_area(is_h2oice_sublim_ini .and. h2o_ice > 0._dp), & coverage_ref,subslope_area_max,stopcrit%surface_h2oice_change) ! Appearance of H2O ice outside of what the PCM announced call check_area_criterion('surface of the new H2O ice reservoirs',tol_h2oice_surf,0._dp, & masked_area(.not. is_h2oice_ini .and. h2o_ice > h2oice_resvr_min), & coverage_ref,subslope_area_max,stopcrit%newresvr_h2oice_change) END SUBROUTINE stopping_crit_h2oice !======================================================================= !======================================================================= SUBROUTINE report_stopping_crit_h2oice(h2oice_sublim_coverage_ini,coverage_ref) !----------------------------------------------------------------------- ! NAME ! report_stopping_crit_h2oice ! ! DESCRIPTION ! Report the reference areas of the H2O ice criteria and the changes ! that will be accepted. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use slopes, only: subslope_area_max ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] ! CODE ! ---- call report_area_criterion('sublimating H2O ice surface',tol_h2oice_surf,h2oice_sublim_coverage_ini,coverage_ref,subslope_area_max) call report_area_criterion('surface of the new H2O ice reservoirs',tol_h2oice_surf,0._dp,coverage_ref,subslope_area_max) END SUBROUTINE report_stopping_crit_h2oice !======================================================================= !======================================================================= SUBROUTINE stopping_crit_co2ice(co2ice_sublim_coverage_ini,is_co2ice_sublim_ini,is_co2ice_ini,coverage_ref,co2_ice,stopcrit) !----------------------------------------------------------------------- ! NAME ! stopping_crit_co2ice ! ! DESCRIPTION ! Check if the CO2 ice which was sublimating has disappeared too much, or if new CO2 ice reservoirs have appeared. ! ! AUTHORS & DATE ! R. Vandemeulebrouck, 07/2022 ! L. Lange, 09/2022 ! JB Clement, 09/2023 ! ! NOTES ! This criterion is not checked in 1D: what matters is the mass of ! CO2 exchanged with the atmosphere which is measured by the criterion ! on the surface pressure. ! The set of the cells expected to hold CO2 ice is frozen at ! initialization while 'd_co2ice' are updated throughout the run. ! So a cell which starts condensing later on is reported as a new ! reservoir. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid use slopes, only: masked_area, subslope_area_max ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- logical(k4), dimension(:,:), intent(in) :: is_co2ice_sublim_ini ! Initial location of the sublimating CO2 ice logical(k4), dimension(:,:), intent(in) :: is_co2ice_ini ! Initial location of the CO2 ice reservoirs (present or announced condensing) real(dp), dimension(:,:), intent(in) :: co2_ice ! Current CO2 surface ice [kg/m2] real(dp), intent(in) :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] class(stopFlags), intent(inout) :: stopcrit ! CODE ! ---- ! Check to escape the subroutine (not relevant in 1D or if the PEM should stop already) if (ngrid == 1 .or. stopcrit%stop_code() > 0) return ! Disappearance of the patch which was sublimating at initialization call check_area_criterion('sublimating CO2 ice surface',tol_co2ice_surf,co2ice_sublim_coverage_ini, & masked_area(is_co2ice_sublim_ini .and. co2_ice > 0._dp), & coverage_ref,subslope_area_max,stopcrit%surface_co2ice_change) ! Appearance of CO2 ice outside of what the PCM announced call check_area_criterion('surface of the new CO2 ice reservoirs',tol_co2ice_surf,0._dp, & masked_area(.not. is_co2ice_ini .and. co2_ice > 0._dp), & coverage_ref,subslope_area_max,stopcrit%newresvr_co2ice_change) END SUBROUTINE stopping_crit_co2ice !======================================================================= !======================================================================= SUBROUTINE report_stopping_crit_co2ice(co2ice_sublim_coverage_ini,coverage_ref) !----------------------------------------------------------------------- ! NAME ! report_stopping_crit_co2ice ! ! DESCRIPTION ! Report the reference areas of the CO2 ice criteria and the changes ! that will be accepted. ! ! AUTHORS & DATE ! JB Clement, 08/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use slopes, only: subslope_area_max ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2] real(dp), intent(in) :: coverage_ref ! Reference area for the tolerance [m2] ! CODE ! ---- call report_area_criterion('sublimating CO2 ice surface',tol_co2ice_surf,co2ice_sublim_coverage_ini,coverage_ref,subslope_area_max) call report_area_criterion('surface of the new CO2 ice reservoirs',tol_co2ice_surf,0._dp,coverage_ref,subslope_area_max) END SUBROUTINE report_stopping_crit_co2ice !======================================================================= !======================================================================= SUBROUTINE stopping_crit_pressure(ps_avg_glob_ini,ps_avg_global,stopcrit) !----------------------------------------------------------------------- ! NAME ! stopping_crit_pressure ! ! DESCRIPTION ! Check if pressure has changed too much. ! ! AUTHORS & DATE ! R. Vandemeulebrouck, 07/2022 ! L. Lange, 09/2022 ! JB Clement, 09/2023 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use display, only: print_msg, LVL_NFO, LVL_WRN use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: ps_avg_glob_ini, ps_avg_global class(stopFlags), intent(inout) :: stopcrit ! CODE ! ---- ! Check to escape the subroutine (not relevant if the PEM should stop already) if (stopcrit%stop_code() > 0) return if (ps_avg_glob_ini <= 0._dp) then call print_msg("Initial global pressure is null or negative. The stopping criterion for pressure will not be checked!",LVL_WRN) return end if ! Check of the criterion if (abs(ps_avg_global - ps_avg_glob_ini) > tol_ps_avg*ps_avg_glob_ini) then stopcrit%pressure_change = .true. call print_msg("Global pressure (ini|now) [Pa] = "//real2str(ps_avg_glob_ini)//' | '//real2str(ps_avg_global),LVL_NFO) call print_msg("Change (accepted|current) [%] = +/- "//real2str(tol_ps_avg*100._dp)//' | '//real2str(100._dp*(ps_avg_global - ps_avg_glob_ini)/ps_avg_glob_ini),LVL_NFO) end if END SUBROUTINE stopping_crit_pressure !======================================================================= END MODULE stopping_crit_core