MODULE criterion_pem_stop_mod

implicit none

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Criterions to check if the PEM needs to call the PCM
!!! Author: RV & LL, 02/2023
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)

use time_evol_mod, only: water_ice_criterion
use comslope_mod,  only: subslope_dist, nslope

implicit none

!=======================================================================
!
!  Routine that checks if the water ice criterion to stop the PEM is reached
!
!=======================================================================

!   arguments:
!   ----------
!   INPUT
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) :: qsurf           ! Actual density of water ice
real,                          intent(in) :: ini_surf        ! Initial surface of h2o ice that was sublimating
real, dimension(ngrid,nslope), intent(in) :: initial_h2o_ice ! Grid point that initialy were covered by h2o_ice
!   OUTPUT
logical, intent(out) :: STOPPING ! Is the criterion reached?
!   local:
!   ------
integer :: i, islope    ! Loop
real    :: present_surf ! Initial/Actual surface of water ice

!=======================================================================
! Initialisation to false
STOPPING = .false.

! Computation of the present surface of water ice sublimating
present_surf = 0.
do i = 1,ngrid
    do islope = 1,nslope
        !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
        if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
    enddo
enddo

! Check of the criterion
if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then
    STOPPING = .true.
    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
    write(*,*) "Current surface of water ice sublimating =", present_surf
    write(*,*) "Initial surface of water ice sublimating =", ini_surf
    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion))
endif

if (ini_surf < 1.e-5 .and. ini_surf > -1.e-5) STOPPING = .false.

END SUBROUTINE criterion_waterice_stop

! ----------------------------------------------------------------------

SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)

use time_evol_mod, only: co2_ice_criterion, ps_criterion
use comslope_mod,  only: subslope_dist

implicit none

!=======================================================================
!
!  Routine that checks if the co2 and pressure criteria to stop the PEM are reached
!
!=======================================================================

!   arguments:
!   ----------
!   INPUT
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) :: qsurf                ! Actual density of water ice
real,                          intent(in) :: ini_surf             ! Initial surface of co2 ice that was sublimating
real, dimension(ngrid,nslope), intent(in) :: initial_co2_ice      ! Grid point that initialy were covered by co2_ice
real,                          intent(in) :: global_ave_press_GCM ! Planet average pressure from the PCM start files
real,                          intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations
!   OUTPUT
logical, intent(out) :: STOPPING_ice ! Is the criterion for co2 ice reached?
logical, intent(out) :: STOPPING_ps  ! Is the criterion for pressure reached?
!   local:
!   ------
integer :: i, islope    ! Loop
real    :: present_surf ! Initial/Actual surface of water ice

!=======================================================================
! Initialisation to false
STOPPING_ice = .false.
STOPPING_ps = .false.

! Computation of the actual surface
present_surf = 0.
do i = 1,ngrid
    do islope = 1,nslope
        if (initial_co2_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
    enddo
enddo

! Check of the criterion
if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then
    STOPPING_ice = .true.
    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
    write(*,*) "Current surface of co2 ice sublimating =", present_surf
    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion))
endif

if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.

if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
    STOPPING_ps = .true.
    write(*,*) "Reason of stopping: the global pressure reach the threshold"
    write(*,*) "Current global pressure =", global_ave_press_new
    write(*,*) "PCM global pressure =", global_ave_press_GCM
    write(*,*) "Percentage of change accepted =", ps_criterion*100.
    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion))
endif

END SUBROUTINE criterion_co2_stop

END MODULE