source: trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 @ 3996

MODULE stopping_crit
!-----------------------------------------------------------------------
! NAME
!     stopping_crit
!
! DESCRIPTION
!     Stopping criteria for PEM.
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! MODULE VARIABLES
! ----------------
real :: h2oice_crit ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
real :: co2ice_crit ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
real :: ps_crit     ! Percentage of change of averaged surface pressure before stopping the PEM

type :: stopFlags
    logical :: surface_h2oice_change = .false.
    logical :: zero_dh2oice          = .false.
    logical :: surface_co2ice_change = .false.
    logical :: pressure_change       = .false.
    logical :: nyears_max_reached    = .false.
    logical :: nyears_maxorb_reached = .false.
    logical :: nyears_tot_reached    = .false.
    logical :: time_limit_reached    = .false.
contains
    procedure :: is_any_set
    procedure :: stop_code
    procedure :: stop_message
end type

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
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, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
class(stopFlags), intent(in) :: this    ! Stopping flags object
logical                      :: stopPEM ! True if any flag is set

! CODE
! ----
stopPEM = this%surface_h2oice_change .or. &
          this%zero_dh2oice          .or. &
          this%surface_co2ice_change .or. &
          this%pressure_change       .or. &
          this%nyears_max_reached    .or. &
          this%nyears_maxorb_reached .or. &
          this%nyears_tot_reached    .or. &
          this%time_limit_reached

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, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
class(StopFlags), intent(in) :: this ! Stopping flags object
integer                      :: code ! Code corresponding to active flag (0 if none)

! CODE
! ----
code = 0
if     (this%surface_h2oice_change) then; code = 1
elseif (this%zero_dh2oice)          then; code = 2
elseif (this%surface_co2ice_change) then; code = 3
elseif (this%pressure_change)       then; code = 4
elseif (this%nyears_max_reached)    then; code = 5
elseif (this%nyears_maxorb_reached) then; code = 6
elseif (this%nyears_tot_reached)    then; code = 7
elseif (this%time_limit_reached)    then; code = 8
endif

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, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
class(StopFlags), intent(in) :: this ! Stopping flags object
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)."
elseif (this%zero_dh2oice)          then; msg = " **** STOPPING because tendencies on h2o ice = 0 (see message above)."
elseif (this%surface_co2ice_change) then; msg = " **** STOPPING because surface of co2 ice has changed too much (see message above)."
elseif (this%pressure_change)       then; msg = " **** STOPPING because surface global pressure has changed too much (see message above)."
elseif (this%nyears_max_reached)    then; msg = " **** STOPPING because maximum number of years is reached."
elseif (this%nyears_maxorb_reached) then; msg = " **** STOPPING because maximum number of years due to orbital parameters is reached."
elseif (this%nyears_tot_reached)    then; msg = " **** STOPPING because the number of years to be simulated is reached."
elseif (this%time_limit_reached)    then; msg = " **** STOPPING because the job time limit is going to be reached."
endif

END FUNCTION stop_message
!=======================================================================

!=======================================================================
SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid)

!-----------------------------------------------------------------------
! NAME
!     stopping_crit_h2o_ice
!
! DESCRIPTION
!     Check if H2O ice surface change criterion to stop PEM is reached.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use comslope_mod,  only: subslope_dist, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer,                          intent(in)    :: ngrid                ! # of physical grid points
real,    dimension(ngrid),        intent(in)    :: cell_area            ! Area of the cells
real,    dimension(ngrid,nslope), intent(in)    :: h2o_ice              ! Actual density of h2o ice
real,                             intent(in)    :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
logical, dimension(ngrid,nslope), intent(in)    :: is_h2oice_sublim_ini ! Grid points where h2o ice was initially sublimating
type(stopFlags),                  intent(inout) :: stopCrit             ! Stopping criterion

! LOCAL VARIABLES
! ---------------
integer :: i, islope       ! Loop
real    :: h2oice_now_surf ! Current surface of h2o ice

! CODE
! ----
if (stopCrit%stop_code() > 0) return

! Computation of the present surface of h2o ice still sublimating
h2oice_now_surf = 0.
do i = 1,ngrid
    do islope = 1,nslope
        if (is_h2oice_sublim_ini(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)
    enddo
enddo

! Check of the criterion
if (h2oice_now_surf < h2oice_ini_surf*(1. - h2oice_crit)) then
    stopCrit%surface_h2oice_change = .true.
    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
    write(*,*) "h2oice_now_surf < h2oice_ini_surf*(1. - h2oice_crit)", h2oice_now_surf < h2oice_ini_surf*(1. - h2oice_crit)
    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
    write(*,*) "Percentage of change accepted =", h2oice_crit*100
else if (h2oice_now_surf > h2oice_ini_surf*(1. + h2oice_crit)) then
    stopCrit%surface_h2oice_change = .true.
    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
    write(*,*) "h2oice_now_surf > h2oice_ini_surf*(1. + h2oice_crit)", h2oice_now_surf > h2oice_ini_surf*(1. + h2oice_crit)
    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
    write(*,*) "Percentage of change accepted =", h2oice_crit*100
endif

END SUBROUTINE stopping_crit_h2o_ice
!=======================================================================

!=======================================================================
SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopCrit,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
!-----------------------------------------------------------------------
! NAME
!     stopping_crit_co2
!
! DESCRIPTION
!     Check if CO2 ice surface and pressure change criteria are reached.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use comslope_mod,  only: subslope_dist

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer,                          intent(in)    :: ngrid, nslope          ! # of grid physical grid points
real,    dimension(ngrid),        intent(in)    :: cell_area              ! Area of the cells
real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                ! Actual density of co2 ice
real,                             intent(in)    :: co2ice_sublim_surf_ini ! Initial surface of sublimatingco2 ice
logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_sublim_ini   ! Grid points where co2 ice was initially sublimating
real,                             intent(in)    :: ps_avg_global_ini      ! Planet average pressure from the PCM start files
real,                             intent(in)    :: ps_avg_global          ! Planet average pressure from the PEM computations
type(stopFlags),                  intent(inout) :: stopCrit               ! Stopping criterion

! LOCAL VARIABLES
! ---------------
integer :: i, islope       ! Loop
real    :: co2ice_now_surf ! Current surface of co2 ice

! CODE
! ----
if (stopCrit%stop_code() > 0) return

! Computation of the present surface of co2 ice still sublimating
co2ice_now_surf = 0.
do i = 1,ngrid
    do islope = 1,nslope
        if (is_co2ice_sublim_ini(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)
    enddo
enddo

! Check of the criterion
if (co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2ice_crit)) then
    stopCrit%surface_co2ice_change = .true.
    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    write(*,*) "co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2ice_crit)", co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2ice_crit)
    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
    write(*,*) "Percentage of change accepted =", co2ice_crit*100.
else if (co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2ice_crit)) then
    stopCrit%surface_co2ice_change = .true.
    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    write(*,*) "co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2ice_crit)", co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2ice_crit)
    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
    write(*,*) "Percentage of change accepted =", co2ice_crit*100.
endif

if (ps_avg_global < ps_avg_global_ini*(1. - ps_crit)) then
    stopCrit%pressure_change = .true.
    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
    write(*,*) "ps_avg_global < ps_avg_global_ini*(1. - ps_crit)", ps_avg_global < ps_avg_global_ini*(1. - ps_crit)
    write(*,*) "Initial global pressure =", ps_avg_global_ini
    write(*,*) "Current global pressure =", ps_avg_global
    write(*,*) "Percentage of change accepted =", ps_crit*100.
else if (ps_avg_global > ps_avg_global_ini*(1. + ps_crit)) then
    stopCrit%pressure_change = .true.
    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
    write(*,*) "ps_avg_global > ps_avg_global_ini*(1. + ps_crit)", ps_avg_global > ps_avg_global_ini*(1. + ps_crit)
    write(*,*) "Initial global pressure =", ps_avg_global_ini
    write(*,*) "Current global pressure =", ps_avg_global
    write(*,*) "Percentage of change accepted =", ps_crit*100.
endif

END SUBROUTINE stopping_crit_co2
!=======================================================================

!=======================================================================
SUBROUTINE stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
!-----------------------------------------------------------------------
! NAME
!     stopping_crit_h2o
!
! DESCRIPTION
!     Check if PEM oscillates infinitely with H2O (no net exchange).
!
! AUTHORS & DATE
!     L. Lange, 2024
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use evolution,      only: dt
use comslope_mod,   only: subslope_dist, def_slope_mean
use phys_constants, only: pi

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer,                       intent(in)    :: ngrid, nslope            ! # of grid points, # of subslopes
real, dimension(ngrid),        intent(in)    :: cell_area                ! Area of each mesh grid (m^2)
real, dimension(ngrid),        intent(in)    :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
real, dimension(ngrid),        intent(in)    :: delta_h2o_icetablesublim ! Mass of H2O that condensed/sublimated at the ice table
real, dimension(ngrid,nslope), intent(in)    :: h2o_ice                  ! H2O ice (kg/m^2)
real, dimension(ngrid,nslope), intent(in)    :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
real,                          intent(out)   :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
type(stopFlags),               intent(inout) :: stopCrit                 ! Stopping criterion

! LOCAL VARIABLES
! ---------------
integer :: i, islope ! Loop

! CODE
! ----
! Checks to escape the subroutine (adjusment not relevant in 1D or if the PEM should stop already)
if (ngrid == 1 .or. stopCrit%stop_code() > 0) then
    S_atm_2_h2o = 1.
    S_h2o_2_atm = S_atm_2_h2o
    S_atm_2_h2oice = 1.
    S_h2oice_2_atm = S_atm_2_h2oice
    return
endif

! We compute the amount of water going out from the atmosphere (S_atm_2_h2o) and going into the atmophere (S_h2o_2_atm)
S_atm_2_h2o = 0.
S_h2o_2_atm = 0.
S_atm_2_h2oice = 0.
S_h2oice_2_atm = 0.
do i = 1,ngrid
    if (delta_h2o_adsorbed(i) > 0.) then
        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_adsorbed(i)*cell_area(i)
    else
        S_h2o_2_atm = S_h2o_2_atm + delta_h2o_adsorbed(i)*cell_area(i)
    endif
    if (delta_h2o_icetablesublim(i) > 0.) then
        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_icetablesublim(i)*cell_area(i)
    else
        S_h2o_2_atm = S_h2o_2_atm + delta_h2o_icetablesublim(i)*cell_area(i)
    endif
    do islope = 1,nslope
        if (d_h2oice(i,islope) > 0.) then
            S_atm_2_h2o = S_atm_2_h2o + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
            S_atm_2_h2oice = S_atm_2_h2oice + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
        else if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
            S_h2o_2_atm = S_h2o_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
            S_h2oice_2_atm = S_h2oice_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
        endif
    enddo
enddo

! Since relative atmospheric water is kept constant, we need to equate condensing reservoirs to the sublimating ones.
! It is not possible if one of them is missing.
if ((abs(S_atm_2_h2oice) < 1.e-10 .or. abs(S_h2oice_2_atm) < 1.e-10)) then
    write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
    write(*,*) "This can be due to the absence of h2o ice in the PCM run."
    write(*,*) "Amount of condensing ice  =", S_atm_2_h2oice
    write(*,*) "Amount of sublimating ice =", S_h2oice_2_atm
    stopCrit%zero_dh2oice = .true.
endif

END SUBROUTINE stopping_crit_h2o
!=======================================================================

END MODULE stopping_crit