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

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4

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

! PARAMETERS
! ----------
real(dp), protected :: h2oice_crit      ! Change of the surface of H2O ice sublimating before stopping the PEM
real(dp), protected :: co2ice_crit      ! Change of the surface of CO2 ice sublimating before stopping the PEM
real(dp), protected :: ps_crit          ! Change of averaged surface pressure before stopping the PEM
real(dp), protected :: lakesurf_crit    ! Change of the surface of lakes before stopping the PEM
real(dp), protected :: lakedistrib_crit ! Change of the distribution of lakes before stopping the PEM

! VARIABLES
! ---------
type :: stopFlags
    logical(k4) :: surface_lake_change       = .false.
    logical(k4) :: distrib_lake_change       = .false.
    logical(k4) :: surface_h2oice_change     = .false.
    logical(k4) :: h2o_flux_balance_unclosed = .false.
    logical(k4) :: surface_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.
contains
    procedure :: is_any_set
    procedure :: stop_code
    procedure :: stop_message
end type

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

!=======================================================================
SUBROUTINE set_stopping_crit_config(h2oice_crit_in,co2ice_crit_in,lakesurf_crit_in,lakedistrib_crit_in,ps_crit_in)
!-----------------------------------------------------------------------
! NAME
!     set_stopping_crit_config
!
! DESCRIPTION
!     Setter for "stopping_crit" 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) :: h2oice_crit_in, co2ice_crit_in, lakesurf_crit_in, lakedistrib_crit_in, ps_crit_in

! CODE
! ----
h2oice_crit = h2oice_crit_in
lakesurf_crit = lakesurf_crit_in
co2ice_crit = co2ice_crit_in
lakedistrib_crit = lakedistrib_crit_in
ps_crit = ps_crit_in
call print_msg('h2oice_crit      = '//real2str(h2oice_crit),LVL_NFO)
call print_msg('co2ice_crit      = '//real2str(co2ice_crit),LVL_NFO)
call print_msg('lakesurf_crit    = '//real2str(lakesurf_crit),LVL_NFO)
call print_msg('lakedistrib_crit = '//real2str(lakedistrib_crit),LVL_NFO)
call print_msg('ps_crit          = '//real2str(ps_crit),LVL_NFO)
if (h2oice_crit <= 0._dp .or. h2oice_crit > 1._dp) call stop_clean(__FILE__,__LINE__,'''h2oice_crit'' outside admissible range ]0.,1.]!',1)
if (co2ice_crit <= 0._dp .or. co2ice_crit > 1._dp) call stop_clean(__FILE__,__LINE__,'''co2ice_crit'' outside admissible range ]0.,1.]!',1)
if (lakesurf_crit <= 0._dp .or. lakesurf_crit > 1._dp) call stop_clean(__FILE__,__LINE__,'''lakesurf_crit'' outside admissible range ]0.,1.]!',1)
if (lakedistrib_crit <= 0._dp .or. lakedistrib_crit > 1._dp) call stop_clean(__FILE__,__LINE__,'''lakedistrib_crit'' outside admissible range ]0.,1.]!',1)
if (ps_crit <= 0._dp .or. ps_crit > 1._dp) call stop_clean(__FILE__,__LINE__,'''ps_crit'' 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, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
stopPEM = this%surface_h2oice_change     .or. &
          this%surface_co2ice_change     .or. &
          this%surface_lake_change       .or. &
          this%distrib_lake_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

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(di)                  :: code ! Code corresponding to active flag (0 if none)

! CODE
! ----
code = 0
if      (this%surface_h2oice_change)     then; code = 1
else if (this%h2o_flux_balance_unclosed) then; code = 2
else if (this%surface_co2ice_change)     then; code = 3
else if (this%pressure_change)           then; code = 4
else if (this%nmax_yr_run_reached)       then; code = 5
else if (this%nmax_yr_runorb_reached)    then; code = 6
else if (this%nmax_yr_sim_reached)       then; code = 7
else if (this%time_limit_reached)        then; code = 8
else if (this%surface_lake_change)       then; code = 9
else if (this%distrib_lake_change)       then; code = 10
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, 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)."
else if (this%surface_co2ice_change)     then; msg = "**** STOPPING because surface of CO2 ice has changed too much (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%h2o_flux_balance_unclosed) then; msg = "**** STOPPING because H2O flux imbalance cannot be closed under physical bounds (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%surface_lake_change)       then; msg = "**** STOPPING because surface of lakes has changed too much (see message above)."
else if (this%distrib_lake_change)       then; msg = "**** STOPPING because distribution of lakes has changed too much (see message above)."
else
    msg = "**** STOPPING for unknown reason (this should not happen!)."
end if

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

!=======================================================================
SUBROUTINE stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit)

!-----------------------------------------------------------------------
! NAME
!     stopping_crit_h2oice
!
! DESCRIPTION
!     Check if H2O ice surface area coverage has changed too much.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use slopes,   only: subslope_dist, def_slope_mean
use geometry, only: ngrid, nslope, cell_area
use maths,    only: pi
use display,  only: print_msg, LVL_NFO, LVL_WRN
use utility,  only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)    :: h2o_ice, d_h2oice
real(dp),                 intent(in)    :: h2oice_sublim_coverage_ini
type(stopFlags),          intent(inout) :: stopcrit

! LOCAL VARIABLES
! ---------------
integer(di) :: i, islope
real(dp)    :: h2oice_now_surf

! CODE
! ----
! Check to escape the subroutine (if the PEM should stop already)
if (stopcrit%stop_code() > 0) return
if (h2oice_sublim_coverage_ini <= 0._dp) then
    call print_msg("Initial surface of H2O ice sublimating is null or negative. The stopping criterion for H2O ice will not be checked!",LVL_WRN)
    return
end if

! Computation of the current surface of sublimating H2O ice
h2oice_now_surf = 0._dp
do i = 1,ngrid
    do islope = 1,nslope
        if (h2o_ice(i,islope) > 0._dp .and. d_h2oice(i,islope) < 0._dp) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
    end do
end do

! Check of the criterion
if (abs(h2oice_now_surf - h2oice_sublim_coverage_ini) > h2oice_crit*h2oice_sublim_coverage_ini) then
    stopcrit%surface_h2oice_change = .true.
    call print_msg("Surface of H2O ice sublimating (ini|now) [m2] = "//real2str(h2oice_sublim_coverage_ini)//' | '//real2str(h2oice_now_surf),LVL_NFO)
    call print_msg("Change (accepted|current) [%] = +/- "//real2str(h2oice_crit*100._dp)//' | '//real2str(100._dp*(h2oice_now_surf - h2oice_sublim_coverage_ini)/h2oice_sublim_coverage_ini),LVL_NFO)
end if

END SUBROUTINE stopping_crit_h2oice
!=======================================================================

!=======================================================================
SUBROUTINE stopping_crit_co2ice(co2ice_sublim_coverage_ini,co2_ice,d_co2ice,stopcrit)
!-----------------------------------------------------------------------
! NAME
!     stopping_crit_co2ice
!
! DESCRIPTION
!     Check if CO2 ice surface area coverage has changed too much.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use slopes,   only: subslope_dist, def_slope_mean
use geometry, only: ngrid, nslope, cell_area
use maths,    only: pi
use display,  only: print_msg, LVL_NFO, LVL_WRN
use utility,  only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)    :: co2_ice, d_co2ice
real(dp),                 intent(in)    :: co2ice_sublim_coverage_ini
type(stopFlags),          intent(inout) :: stopcrit

! LOCAL VARIABLES
! ---------------
integer(di) :: i, islope
real(dp)    :: co2ice_now_surf ! Current surface of CO2 ice

! 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
if (co2ice_sublim_coverage_ini <= 0._dp) then
    call print_msg("Initial surface of CO2 ice sublimating is null or negative. The stopping criterion for CO2 ice will not be checked!",LVL_WRN)
    return
end if

! Computation of the current surface of CO2 ice still sublimating
co2ice_now_surf = 0.
do i = 1,ngrid
    do islope = 1,nslope
        if (co2_ice(i,islope) > 0._dp .and. d_co2ice(i,islope) < 0._dp) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
    end do
end do

! Check of the criterion
if (abs(co2ice_now_surf - co2ice_sublim_coverage_ini) > co2ice_crit*co2ice_sublim_coverage_ini) then
    stopcrit%surface_co2ice_change = .true.
    call print_msg("Surface of CO2 ice sublimating (ini|now) [m2] = "//real2str(co2ice_sublim_coverage_ini)//' | '//real2str(co2ice_now_surf),LVL_NFO)
    call print_msg("Change (accepted|current) [%] = +/- "//real2str(co2ice_crit*100._dp)//' | '//real2str(100._dp*(co2ice_now_surf - co2ice_sublim_coverage_ini)/co2ice_sublim_coverage_ini),LVL_NFO)
end if

END SUBROUTINE 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
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! 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
type(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) > ps_crit*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(ps_crit*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
!=======================================================================

!=======================================================================
SUBROUTINE stopping_crit_lakes(lake_frac_ini,lake_frac,stopcrit)
!-----------------------------------------------------------------------
! NAME
!     stopping_crit_lakes
!
! DESCRIPTION
!     Check if lake surface area coverage has changed too much.
!
! AUTHORS & DATE
!     A. Gauvain, 06/2026
!     JB. Clement, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, cell_area
use display,  only: print_msg, LVL_NFO, LVL_WRN
use utility,  only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in)    :: lake_frac_ini, lake_frac
type(stopFlags),        intent(inout) :: stopcrit

! LOCAL VARIABLES
! ---------------
real(dp) :: lake_ini_surf, lake_now_surf ! Initial and current surface of lakes
real(dp) :: distrib_change               ! Change of the distribution of lakes

! 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

! Computation of the initial and current surface of lakes
lake_ini_surf = sum(lake_frac_ini(:)*cell_area(:))
lake_now_surf = sum(lake_frac(:)*cell_area(:))
if (lake_ini_surf <= 0._dp) then
    call print_msg("Initial surface of lakes is null or negative. The stopping criterion for lakes will not be checked!",LVL_WRN)
    return
end if

! Check of the criterion for surface change
if (abs(lake_now_surf - lake_ini_surf) > lakesurf_crit*lake_ini_surf) then
    stopcrit%surface_lake_change = .true.
    call print_msg("Surface of lakes (ini|now) [m2] = "//real2str(lake_ini_surf)//' | '//real2str(lake_now_surf),LVL_NFO)
    call print_msg("Change (accepted|current) [%] = +/- "//real2str(lakesurf_crit*100._dp)//' | '//real2str(100._dp*(lake_now_surf - lake_ini_surf)/lake_ini_surf),LVL_NFO)
end if

! Check of the criterion for distribution change
distrib_change = 0.5_dp*sum(abs(lake_frac(:) - lake_frac_ini(:)) * cell_area(:)) ! Divided by 2 to avoid double counting gain and loss of lake area
if (distrib_change > lakedistrib_crit*lake_ini_surf) then
    stopcrit%distrib_lake_change = .true.
    call print_msg("Lake redistribution [% of initial lake area] = "//real2str(100._dp*distrib_change/lake_ini_surf), LVL_NFO)
end if

END SUBROUTINE stopping_crit_lakes
!=======================================================================

END MODULE stopping_crit