module criterion_pem_stop_mod
  implicit none

  contains


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


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

  USE temps_mod_evol, 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 grid physical grid points 
  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
  REAL,    intent(in) :: qsurf(ngrid,nslope)           ! physical point field : Actual density of water ice
  REAL,    intent(in) :: ini_surf                      ! Initial surface of h2o ice that was sublimating
  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Grid point that initialy were covered by h2o_ice

!   OUTPUT
  LOGICAL, intent(out) :: STOPPING              ! Logical : 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).GT.0.5 .and. qsurf(i,islope).GT.0.) then
         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
      endif
    enddo
  enddo
  
!   check of the criterion
  if(present_surf.LT.ini_surf*(1-water_ice_criterion) .OR. &
     present_surf.GT.ini_surf*(1+water_ice_criterion)) then
    STOPPING=.TRUE. 
    print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
    print *, "Current surface of water ice sublimating=", present_surf
    print *, "Initial surface of water ice sublimating=", ini_surf
    print *, "Percentage of change accepted=", water_ice_criterion*100
    print *, "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion))
  endif

  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
    STOPPING=.FALSE.
  endif
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 temps_mod_evol, ONLY: co2_ice_criterion,ps_criterion
  use comslope_mod, ONLY: subslope_dist

      IMPLICIT NONE

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

!   arguments:
!   ----------

!   INPUT
  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
  REAL,    intent(in) :: ini_surf                      ! Initial surface of co2 ice that was sublimating
  REAL,    intent(in) :: initial_co2_ice(ngrid,nslope) ! Grid point that initialy were covered by co2_ice
  REAL,    intent(in) :: global_ave_press_GCM          ! Planet average pressure from the GCM start files 
  REAL,    intent(in) :: global_ave_press_new          ! Planet average pressure from the PEM computations

!   OUTPUT
  LOGICAL, intent(out) :: STOPPING_ice              ! Logical : is the criterion for ice reached?
  LOGICAL, intent(out) :: STOPPING_ps               ! Logical : 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).GT.0.5 .and. qsurf(i,islope).GT.0.) then
         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
      endif
   enddo
  enddo
  
!   check of the criterion
  if(present_surf.LT.ini_surf*(1-co2_ice_criterion) .OR. &
     present_surf.GT.ini_surf*(1+co2_ice_criterion)) then
    STOPPING_ice=.TRUE. 
    print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
    print *, "Current surface of co2 ice sublimating=", present_surf
    print *, "Initial surface of co2 ice sublimating=", ini_surf
    print *, "Percentage of change accepted=", co2_ice_criterion*100
    print *, "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion))
  endif

  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
       STOPPING_ice=.FALSE.
  endif

  if(global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion) .OR. &
     global_ave_press_new.GT.global_ave_press_GCM*(1+ps_criterion)) then
    STOPPING_ps=.TRUE. 
    print *, "Reason of stopping : The global pressure reach the threshold:"
    print *, "Current global pressure=", global_ave_press_new
    print *, "GCM global pressure=", global_ave_press_GCM
    print *, "Percentage of change accepted=", ps_criterion*100
    print *, "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion))
  endif

END SUBROUTINE criterion_co2_stop


END MODULE