MODULE stopping_crit_mod

implicit none

!=======================================================================
contains
!=======================================================================

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

SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)

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

implicit none

!=======================================================================
!
! Routine to check if the h2o 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) :: h2o_ice   ! Actual density of h2o ice
real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
!   OUTPUT
integer, intent(inout) :: stopPEM ! Stopping criterion code
!   local:
!   ------
integer :: i, islope ! Loop
real    :: surf_now  ! Current surface of h2o ice

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

! Check of the criterion
if (surf_now < surf_ini*(1. - h2o_ice_crit)) then
    stopPEM = 1
    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
    write(*,*) "surf_now < surf_ini*(1. - h2o_ice_crit)", surf_now < surf_ini*(1. - h2o_ice_crit)
    write(*,*) "Current surface of h2o ice sublimating =", surf_now
    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
else if (surf_now > surf_ini*(1. + h2o_ice_crit)) then
    stopPEM = 1
    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
    write(*,*) "surf_now > surf_ini*(1. + h2o_ice_crit)", surf_now > surf_ini*(1. + h2o_ice_crit)
    write(*,*) "Current surface of h2o ice sublimating =", surf_now
    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
endif

if (abs(surf_ini) < 1.e-5) stopPEM = 0

END SUBROUTINE stopping_crit_h2o_ice

!=======================================================================

SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)

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

implicit none

!=======================================================================
!
! Routine to check 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) :: co2_ice              ! Actual density of h2o ice
real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
real,                          intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
!   OUTPUT
integer, intent(inout) :: stopPEM ! Stopping criterion code

!   local:
!   ------
integer :: i, islope ! Loop
real    :: surf_now  ! Current surface of co2 ice

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

! Check of the criterion
if (surf_now < surf_ini*(1. - co2_ice_crit)) then
    stopPEM = 3
    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    write(*,*) "surf_now < surf_ini*(1. - co2_ice_crit)", surf_now < surf_ini*(1. - co2_ice_crit)
    write(*,*) "Current surface of co2 ice sublimating =", surf_now
    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
else if (surf_now > surf_ini*(1. + co2_ice_crit)) then
    stopPEM = 3
    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
    write(*,*) "surf_now > surf_ini*(1. + co2_ice_crit)", surf_now > surf_ini*(1. + co2_ice_crit)
    write(*,*) "Current surface of co2 ice sublimating =", surf_now
    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
endif

if (abs(surf_ini) < 1.e-5) stopPEM = 0

if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
    stopPEM = 4
    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
    write(*,*) "global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)", global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)
    write(*,*) "Current global pressure =", global_avg_press_new
    write(*,*) "PCM global pressure =", global_avg_press_PCM
    write(*,*) "Percentage of change accepted =", ps_criterion*100.
else if (global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)) then
    stopPEM = 4
    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
    write(*,*) "global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)", global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)
    write(*,*) "Current global pressure =", global_avg_press_new
    write(*,*) "PCM global pressure =", global_avg_press_PCM
    write(*,*) "Percentage of change accepted =", ps_criterion*100.
endif

END SUBROUTINE stopping_crit_co2

END MODULE