MODULE ice_table
!-----------------------------------------------------------------------
! NAME
!     ice_table
!
! DESCRIPTION
!     Ice table variables and methods to compute it (dynamic and static).
!
! AUTHORS & DATE
!     L. Lange, 02/2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
logical(k4), protected :: icetable_equilibrium ! Flag to compute the icetable depth according to equilibrium
logical(k4), protected :: icetable_dynamic     ! Flag to compute the icetable depth according to the dynamic method

! VARIABLES
! ---------
real(dp), allocatable, dimension(:,:)   :: icetable_depth     ! Depth of the ice table [m]
real(dp), allocatable, dimension(:,:)   :: icetable_thickness ! Thickness of the ice table [m]
real(dp), allocatable, dimension(:,:,:) :: ice_porefilling    ! Amount of porefilling in each layer in each grid [m^3/m^3]

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

!=======================================================================
SUBROUTINE ini_ice_table()
!-----------------------------------------------------------------------
! NAME
!     ini_ice_table
!
! DESCRIPTION
!     Allocate module arrays for ice table fields.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil

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

! CODE
! ----
allocate(icetable_depth(ngrid,nslope))
allocate(icetable_thickness(ngrid,nslope))
allocate(ice_porefilling(ngrid,nsoil,nslope))

END SUBROUTINE ini_ice_table
!=======================================================================

!=======================================================================
SUBROUTINE end_ice_table
!-----------------------------------------------------------------------
! NAME
!     end_ice_table
!
! DESCRIPTION
!     Deallocate ice table arrays.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
if (allocated(icetable_depth)) deallocate(icetable_depth)
if (allocated(icetable_thickness)) deallocate(icetable_thickness)
if (allocated(ice_porefilling)) deallocate(ice_porefilling)

END SUBROUTINE end_ice_table
!=======================================================================

!=======================================================================
SUBROUTINE set_ice_table_config(icetable_equilibrium_in,icetable_dynamic_in)
!-----------------------------------------------------------------------
! NAME
!     set_ice_table_config
!
! DESCRIPTION
!     Setter for 'ice_table' configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: bool2str
use display, only: print_msg

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

! ARGUMENTS
! ---------
logical(k4), intent(in) :: icetable_equilibrium_in, icetable_dynamic_in

! CODE
! ----
icetable_equilibrium = icetable_equilibrium_in
icetable_dynamic = icetable_dynamic_in
if (icetable_equilibrium .and. icetable_dynamic) then
    call print_msg('Ice table is asked to be computed both by the equilibrium and dynamic method. Then, the dynamic method is chosen.')
    icetable_equilibrium = .false.
end if
call print_msg('icetable_equilibrium = '//bool2str(icetable_equilibrium_in))
call print_msg('icetable_dynamic     = '//bool2str(icetable_dynamic_in))

END SUBROUTINE set_ice_table_config
!=======================================================================

!=======================================================================
SUBROUTINE computeice_table_equilibrium(watercaptag,rhowatersurf_avg,rhowatersoil_avg,regolith_inertia,ice_table_depth,ice_table_thickness)
!-----------------------------------------------------------------------
! NAME
!     computeice_table_equilibrium
!
! DESCRIPTION
!     Compute the ice table depth knowing the yearly average water
!     density at the surface and at depth. Computations are made following
!     the methods in Schorghofer et al., 2005.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 12/12/2025
!
! NOTES
!     This subroutine only gives the ice table at equilibrium and does
!     not consider exchange with the atmosphere.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,           only: ngrid, nslope, nsoil
use maths,              only: findroot
use soil,               only: mlayer, layer ! Depth of the vertical grid
use soil_therm_inertia, only: get_ice_TI

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

! ARGUMENTS
! ---------
logical(k4), dimension(:),     intent(in)  :: watercaptag         ! Boolean to check the presence of a perennial glacier
real(dp),    dimension(:,:),   intent(in)  :: rhowatersurf_avg    ! Water density at the surface, yearly averaged [kg/m^3]
real(dp),    dimension(:,:,:), intent(in)  :: rhowatersoil_avg    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
real(dp),    dimension(:,:),   intent(in)  :: regolith_inertia    ! Thermal inertia of the regolith layer [SI]
real(dp),    dimension(:,:),   intent(out) :: ice_table_depth     ! Ice table depth [m]
real(dp),    dimension(:,:),   intent(out) :: ice_table_thickness ! Ice table thickness [m]

! LOCAL VARIABLES
! ---------------
integer(di)                       :: ig, islope, isoil, isoilend
real(dp), dimension(nsoil)        :: diff_rho                ! Difference of water vapor density between the surface and at depth [kg/m^3]
real(dp)                          :: ice_table_end           ! Depth of the end of the ice table  [m]
real(dp), dimension(ngrid,nslope) :: previous_icetable_depth ! Ice table computed at previous ice depth [m]
real(dp)                          :: stretch                 ! Stretch factor to improve the convergence of the ice table
real(dp)                          :: wice_inertia            ! Water Ice thermal Inertia [USI]

! CODE
! ----
previous_icetable_depth = ice_table_depth
do ig = 1,ngrid
    if (watercaptag(ig)) then
        ice_table_depth(ig,:) = 0._dp
        ice_table_thickness(ig,:) = layer(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    else
        do islope = 1,nslope
            ice_table_depth(ig,islope) = -1._dp
            ice_table_thickness(ig,islope) = 0._dp
            do isoil = 1,nsoil
                diff_rho(isoil) = rhowatersurf_avg(ig,islope) - rhowatersoil_avg(ig,isoil,islope)
            end do
            if (diff_rho(1) > 0.) then ! ice is at the surface
                ice_table_depth(ig,islope) = 0._dp
                do isoilend = 2,nsoil - 1
                    if (diff_rho(isoilend) > 0._dp .and. diff_rho(isoilend + 1) < 0._dp) then
                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer(isoilend - 1),mlayer(isoilend),ice_table_end)
                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                        exit
                    end if
                end do
            else
                do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
                    if (diff_rho(isoil) < 0._dp .and. diff_rho(isoil + 1) > 0._dp) then
                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer(isoil),mlayer(isoil + 1),ice_table_depth(ig,islope))
                        ! Now let's find the end of the ice table
                        ice_table_thickness(ig,islope) = layer(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
                        do isoilend = isoil + 1,nsoil - 1
                            if (diff_rho(isoilend) > 0._dp .and. diff_rho(isoilend + 1) < 0._dp) then
                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer(isoilend - 1),mlayer(isoilend),ice_table_end)
                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                                exit
                            end if
                        end do
                    end if ! diff_rho(z) <0 & diff_rho(z+1) > 0
                end do
            end if ! diff_rho(1) > 0
        end do
    end if ! watercaptag
end do

! Small trick to speed up the convergence, Oded's idea.
do islope = 1,nslope
    do ig = 1,ngrid
        if (ice_table_depth(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0._dp) then
            call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),wice_inertia)
            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_depth(ig,islope) - &
                                             previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
            ice_table_depth(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch
        end if
    end do
end do

END SUBROUTINE computeice_table_equilibrium
!=======================================================================

!=======================================================================
SUBROUTINE compute_deltam_icetable(icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_icetable)
!-----------------------------------------------------------------------
! NAME
!     compute_deltam_icetable
!
! DESCRIPTION
!     Compute the mass of H2O that has sublimated from the ice table /
!     condensed.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil
use soil,     only: mlayer
use slopes,   only: subslope_dist, def_slope_mean
use planet,   only: porosity
use maths,    only: pi

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:),   intent(in)  :: icetable_thickness_old   ! Ice table thickness at the previous iteration [m]
real(dp),    dimension(:,:,:), intent(in)  :: ice_porefilling_old      ! Ice pore filling at the previous iteration [m]
real(dp),    dimension(:,:),   intent(in)  :: tsurf                    ! Surface temperature [K]
real(dp),    dimension(:,:,:), intent(in)  :: tsoil                    ! Soil temperature [K]
real(dp),    dimension(:),     intent(out) :: delta_icetable ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, islope, isoil
real(dp)    :: Tice ! Ice temperature [k]

! CODE
! ----
! Let's compute the amount of ice that has sublimated in each subslope
if (icetable_equilibrium) then
    delta_icetable = 0._dp
    do ig = 1,ngrid
        do islope = 1,nslope
            call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
            delta_icetable(ig) = delta_icetable(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
        end do
    end do
else if (icetable_dynamic) then
    delta_icetable = 0._dp
    do ig = 1,ngrid
        do islope = 1,nslope
            do isoil = 1,nsoil
                call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),mlayer(isoil - 1),Tice)
                delta_icetable(ig) = delta_icetable(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
            end do
        end do
    end do
end if

END SUBROUTINE compute_deltam_icetable
!=======================================================================

!=======================================================================
SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
!-----------------------------------------------------------------------
! NAME
!     find_layering_icetable
!
! DESCRIPTION
!     Compute layering between dry soil, pore filling ice and ice
!     sheet based on Schorghofer, Icarus (2010). Adapted from NS MSIM.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,           only: mlayer
use geometry,       only: nsoil
use maths,          only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
use subsurface_ice, only: constriction

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

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in)    :: porefill         ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
real(dp), dimension(:), intent(in)    :: psat_soil        ! Soil water pressure at saturation, yearly averaged [Pa]
real(dp),               intent(in)    :: psat_surf        ! Surface water pressure at saturation, yearly averaged [Pa]
real(dp),               intent(in)    :: pwat_surf        ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
real(dp),               intent(in)    :: psat_bottom      ! Boundary conditions for soil vapor pressure [Pa]
real(dp),               intent(in)    :: B                ! Constant (Eq. 8 from  Schorgofer, Icarus (2010).)
integer(di),            intent(in)    :: index_IS         ! Index of the soil layer where the ice sheet begins [1]
real(dp),               intent(inout) :: depth_filling    ! Depth where pore filling begins [m]
integer(di),            intent(out)   :: index_filling    ! Index where the pore filling begins [1]
integer(di),            intent(out)   :: index_geothermal ! Index where the ice table stops [1]
real(dp),               intent(out)   :: depth_geothermal ! Depth where the ice table stops [m]
real(dp), dimension(:), intent(out)   :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase

! LOCAL VARIABLES
! ---------------
real(dp), dimension(nsoil) :: eta                          ! Constriction
real(dp), dimension(nsoil) :: dz_psat                      ! First derivative of the vapor pressure at saturation
real(dp), dimension(nsoil) :: dz_eta                       ! \partial z \eta
real(dp), dimension(nsoil) :: dzz_psat                     ! \partial \partial psat
integer(di)                :: ilay, index_tmp              ! Index for loop
real(dp)                   :: old_depth_filling            ! Depth_filling saved
real(dp)                   :: Jdry                         ! Flux trought the dry layer
real(dp)                   :: Jsat                         ! Flux trought the ice layer
real(dp)                   :: Jdry_prevlay, Jsat_prevlay   ! Same but for the previous ice layer
integer(di)                :: index_firstice               ! First index where ice appears (i.e., f > 0)
real(dp)                   :: massfillabove, massfillafter ! H2O mass above and after index_geothermal
real(dp), parameter        :: pvap2rho = 18.e-3/8.314

! CODE
! ----
! 0. Compute constriction over the layer
! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
if (index_IS < 0) then
    index_tmp = nsoil
    do ilay = 1,nsoil
        eta(ilay) = constriction(porefill(ilay))
    end do
else
    index_tmp = index_IS
    do ilay = 1,index_IS - 1
        eta(ilay) = constriction(porefill(ilay))
    end do
    do ilay = index_IS,nsoil
        eta(ilay) = 0._dp
    end do
end if

! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
old_depth_filling = depth_filling

call deriv1(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dz_psat)

do ilay = 1,index_tmp
    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
    if (Jdry - Jsat <= 0) then
        index_filling = ilay
        exit
    end if
end do

if (index_filling == 1) then
    depth_filling = mlayer(1)
else if (index_filling > 1) then
    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer(index_filling - 1)
    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer(index_filling),mlayer(index_filling - 1),depth_filling)
end if

! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
index_firstice = -1
do ilay = 1,nsoil
    if (porefill(ilay) <= 0._dp) then
        index_firstice = ilay  ! first point with ice
        exit
    end if
end do

! 2.1: now we can compute
call deriv1(mlayer,nsoil,eta,1.,eta(nsoil - 1),dz_eta)

if (index_firstice > 0 .and. index_firstice < nsoil - 2) call deriv1_onesided(index_firstice,mlayer,nsoil,eta,dz_eta(index_firstice))

call deriv2_simple(mlayer,nsoil,psat_soil,psat_surf,psat_bottom,dzz_psat)
dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)

! 3. Ice table boundary due to geothermal heating
if (index_IS > 0) index_geothermal = -1
if (index_geothermal < 0) depth_geothermal = -1._dp
if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
    index_geothermal = -1
    do ilay = 2,nsoil
        if (dz_psat(ilay) > 0._dp) then  ! first point with reversed flux
            index_geothermal = ilay
            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
            exit
        end if
    end do
else
    index_geothermal = -1
end if
if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
    call  colint(porefill/eta,mlayer,nsoil,index_geothermal - 1,nsoil,massfillabove)
    index_tmp = -1
    do ilay = index_geothermal,nsoil
        if (minval(eta(ilay:nsoil)) <= 0._dp) cycle ! eta=0 means completely full
        call colint(porefill/eta,mlayer,nsoil,ilay,nsoil,massfillafter)
        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
            if (ilay > index_geothermal) then
                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer(ilay - 1),mlayer(ilay),depth_geothermal)
!                if (depth_geothermal > mlayer(ilay) .or. depth_geothermal < mlayer(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer(ilay - 1),depth_geothermal,mlayer(ilay)
              index_tmp = ilay
           end if
           ! otherwise leave depth_geothermal unchanged
           exit
        end if
        massfillabove = massfillafter
    end do
    if (index_tmp > 0) index_geothermal = index_tmp
end if

END SUBROUTINE find_layering_icetable
!=======================================================================

!=======================================================================
SUBROUTINE compute_Tice(ptsoil,ptsurf,ice_depth,Tice)
!-----------------------------------------------------------------------
! NAME
!     compute_Tice
!
! DESCRIPTION
!     Compute subsurface ice temperature by interpolating the
!     temperatures between the two adjacent cells.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 12/12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: nsoil
use soil,     only: mlayer
use stoppage, only: stop_clean

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

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in)  :: ptsoil    ! Soil temperature (K)
real(dp),               intent(in)  :: ptsurf    ! Surface temperature (K)
real(dp),               intent(in)  :: ice_depth ! Ice depth (m)
real(dp),               intent(out) :: Tice      ! Ice temperatures (K)

! LOCAL VARIABLES
! ---------------
integer(di) :: ik       ! Loop variables
integer(di) :: indexice ! Index of the ice

! CODE
! ----
indexice = -1
if (ice_depth >= mlayer(nsoil - 1)) then
    Tice = ptsoil(nsoil)
else
    if(ice_depth < mlayer(0)) then
        indexice = 0
    else
        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
            if (mlayer(ik) <= ice_depth .and. mlayer(ik + 1) > ice_depth) then
                indexice = ik + 1
                exit
            end if
        end do
    end if
    if (indexice < 0) then
        call stop_clean(__FILE__,__LINE__,"subsurface ice is below the last soil layer!",1)
    else
        if(indexice >= 1) then ! Linear inteprolation between soil temperature
            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer(indexice - 1) - mlayer(indexice))*(ice_depth - mlayer(indexice)) + ptsoil(indexice + 1)
        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
            Tice = (ptsoil(1) - ptsurf)/mlayer(0)*(ice_depth - mlayer(0)) + ptsoil(1)
        end if ! index ice >= 1
    end if !indexice < 0
end if ! icedepth > mlayer(nsoil - 1)

END SUBROUTINE compute_Tice
!=======================================================================

!=======================================================================
FUNCTION rho_ice(T,ice) RESULT(rho)
!-----------------------------------------------------------------------
! NAME
!     rho_ice
!
! DESCRIPTION
!     Compute subsurface ice density.
!
! AUTHORS & DATE
!     JB Clement, 12/12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

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

! ARGUMENTS
! ---------
real(dp),     intent(in) :: T
character(3), intent(in) :: ice

! LOCAL VARIABLES
! ---------------
real(dp) :: rho

! CODE
! ----
select case (trim(adjustl(ice)))
    case('h2o')
        rho = -3.5353e-4_dp*T**2 + 0.0351_dp*T + 933.5030_dp ! Rottgers 2012
    case('co2')
        rho = (1.72391_dp - 2.53e-4_dp*T - 2.87_dp*1.e-7_dp*T**2)*1.e3_dp ! Mangan et al. 2017
    case default
        call stop_clean(__FILE__,__LINE__,"type of ice unknown!",1)
end select

END FUNCTION rho_ice
!=======================================================================

!=======================================================================
SUBROUTINE evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable, &
                            q_h2o_ts,ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old)
!-----------------------------------------------------------------------
! NAME
!     evolve_ice_table
!
! DESCRIPTION
!     Update the subsurface ice table using either the equilibrium or
!     dynamic method and calculate the resulting H2O mass exchange.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!     C. Metz, 02/2026
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,       only: ngrid, nsoil, nslope, nday
use slopes,         only: subslope_dist, def_slope_mean
use surf_ice,       only: threshold_h2oice_cap
use maths,          only: pi
use soil,           only: TI
use subsurface_ice, only: dyn_ss_ice_m
use display,        only: print_msg

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

! ARGUMENTS
! ---------
real(dp), dimension(:),     intent(in)    :: ps_avg
real(dp), dimension(:,:),   intent(in)    :: h2o_ice, h2o_surfdensity_avg, tsurf_avg, q_h2o_ts
real(dp), dimension(:,:,:), intent(in)    :: h2o_soildensity_avg, tsoil_avg
real(dp), dimension(:,:),   intent(out)   :: icetable_depth_old
real(dp), dimension(:),     intent(out)   :: delta_icetable
real(dp), dimension(:,:),   intent(inout) :: icetable_depth, icetable_thickness
real(dp), dimension(:,:,:), intent(inout) :: ice_porefilling

! LOCAL VARIABLES
! ---------------
integer(di)                                :: i, islope
logical(k4), dimension(ngrid)              :: is_h2o_perice
real(dp),    dimension(ngrid,nslope)       :: icetable_thickness_old
real(dp),    dimension(ngrid,nsoil,nslope) :: ice_porefilling_old

! CODE
! ----
do i = 1,ngrid
    if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
        ! There is enough ice to be considered as an infinite reservoir
        is_h2o_perice(i) = .true.
    else
        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
        is_h2o_perice(i) = .false.
    end if
end do
if (icetable_equilibrium) then
    call print_msg("> Evolution of ice table (equilibrium method)")
    icetable_thickness_old(:,:) = icetable_thickness(:,:)
    call computeice_table_equilibrium(is_h2o_perice,h2o_surfdensity_avg,h2o_soildensity_avg,TI(:,1,:),icetable_depth,icetable_thickness)
else if (icetable_dynamic) then
    call print_msg("> Evolution of ice table (dynamic method)")
    ice_porefilling_old(:,:,:) = ice_porefilling(:,:,:)
    icetable_depth_old(:,:) = icetable_depth(:,:)
    do i = 1,ngrid
        do islope = 1,nslope
            call dyn_ss_ice_m(icetable_depth(i,islope),tsurf_avg(i,islope), tsoil_avg(i,:,islope),nsoil,TI(i,1,islope),ps_avg, &
                              (/sum(q_h2o_ts(i,:))/real(nday,dp)/),ice_porefilling(i,:,islope),ice_porefilling(i,:,islope),icetable_depth(i,islope))
        end do
    end do
end if

! Mass of H2O exchange between the ssi and the atmosphere
call compute_deltam_icetable(icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_avg,delta_icetable)

END SUBROUTINE evolve_ice_table
!=======================================================================

END MODULE ice_table