MODULE evol_ice_mod

implicit none

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

SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)

use time_evol_mod, only: dt
use glaciers_mod,  only: rho_co2ice

implicit none

!=======================================================================
!
! Routine to compute the evolution of CO2 ice
!
!=======================================================================
! arguments:
! ----------
! INPUT
integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes

! OUTPUT
real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! Previous and actual density of CO2 ice
real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of perennial ice over one year
real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]

! local:
! ------
real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
!=======================================================================
! Evolution of CO2 ice for each physical point
write(*,*) '> Evolution of CO2 ice'

co2_ice_old = co2_ice
co2_ice = co2_ice + d_co2ice*dt
where (co2_ice < 0.)
    co2_ice = 0.
    d_co2ice = -co2_ice_old/dt
end where

zshift_surf = d_co2ice*dt/rho_co2ice

END SUBROUTINE evol_co2_ice

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

SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)

use time_evol_mod,     only: dt
use glaciers_mod,      only: rho_h2oice
use stopping_crit_mod, only: stopping_crit_h2o

implicit none

!=======================================================================
!
! Routine to compute the evolution of h2o ice
!
!=======================================================================
! arguments:
! ----------
! INPUT
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 have condensed/sublimated at the ice table (kg/m^2)

! OUTPUT
real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! H2O ice (kg/m^2)
real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
integer,                       intent(inout) :: stopPEM  ! Stopping criterion code
real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]

! local:
! ------
integer                       :: i, islope                                                ! Loop variables
real                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
real, dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
!=======================================================================
write(*,*) '> Evolution of H2O ice'

if (ngrid == 1) then ! In 1D
    h2o_ice = h2o_ice + d_h2oice*dt
    zshift_surf = d_h2oice*dt/rho_h2oice
else ! In 3D
    call 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,stopPEM)
    if (stopPEM > 0) return

    call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
    h2o_ice = h2o_ice + d_h2oice_new*dt
    zshift_surf = d_h2oice_new*dt/rho_h2oice
endif

END SUBROUTINE evol_h2o_ice

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

SUBROUTINE balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)

use time_evol_mod, only: dt

implicit none

!=======================================================================
!
! Routine to balance the H2O flux from and into the atmosphere accross reservoirs
!
!=======================================================================
! arguments:
! ----------
! INPUT
integer,                       intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
real,                          intent(in) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! ! Variables to conserve H2O
real, dimension(ngrid,nslope), intent(in) :: h2o_ice ! H2O ice (kg/m^2)

! OUTPUT
real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
real, dimension(ngrid,nslope), intent(out) :: d_h2oice_new ! Adjusted tendencies to keep balance between donor and recipient reservoirs

! local:
! ------
integer :: i, islope
real    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target
!=======================================================================
S_target = (S_atm_2_h2o + S_h2o_2_atm)/2.
S_target_cond_h2oice = S_atm_2_h2oice + S_target - S_atm_2_h2o
S_target_subl_h2oice = S_h2oice_2_atm + S_target - S_h2o_2_atm

d_h2oice_new = 0.
S_ghostice = 0.
do i = 1,ngrid
    do islope = 1,nslope
        if (d_h2oice(i,islope) > 0.) then ! Condensing
            d_h2oice_new(i,islope) = d_h2oice_new(i,islope)*S_target_cond_h2oice/S_atm_2_h2oice
        else if (d_h2oice(i,islope) < 0.) then ! Sublimating
            d_target = d_h2oice(i,islope)*S_target_subl_h2oice/S_h2oice_2_atm
            if (abs(d_target*dt) < h2o_ice(i,islope)) then ! Enough ice to sublimate everything
                d_h2oice_new(i,islope) = d_target
            else ! Not enough ice to sublimate everything
                ! We sublimate what we can
                d_h2oice_new(i,islope) = h2o_ice(i,islope)/dt
                ! It means the tendency is zero next time
                d_h2oice(i,islope) = 0.
                ! We compute the amount of H2O ice that we could not make sublimate
                S_ghostice = S_ghostice + abs(d_target*dt) - h2o_ice(i,islope)
            endif
        endif
    enddo
enddo

! We need to remove this ice unable to sublimate from places where ice condenseds in order to keep balance
where (d_h2oice_new > 0.) d_h2oice_new = d_h2oice_new*(S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice

END SUBROUTINE balance_h2oice_reservoirs

END MODULE evol_ice_mod