MODULE glaciers
!-----------------------------------------------------------------------
! NAME
!     glaciers
!
! DESCRIPTION
!     Compute flow and transfer of CO2 and H2O ice glaciers on slopes
!     based on maximum ice thickness and basal drag conditions.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
logical(k4), protected :: h2oice_flow ! Flag to compute H2O ice flow
logical(k4), protected :: co2ice_flow ! Flag to compute CO2 ice flow

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

!=======================================================================
SUBROUTINE set_glaciers_config(h2oice_flow_in,co2ice_flow_in)
!-----------------------------------------------------------------------
! NAME
!     set_glaciers_config
!
! DESCRIPTION
!     Setter for 'glaciers' configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
logical(k4), intent(in) :: h2oice_flow_in, co2ice_flow_in

! CODE
! ----
h2oice_flow = h2oice_flow_in
co2ice_flow = co2ice_flow_in
if (h2oice_flow .and. nslope == 1) then
    call print_msg('Warning: h2oice_flow = .true. but nslope = 1. So there will be no flow!')
    h2oice_flow = .false.
end if
if (co2ice_flow .and. nslope == 1) then
    call print_msg('Warning: co2ice_flow = .true. but nslope = 1. So there will be no flow!')
    co2ice_flow = .false.
end if
call print_msg('h2oice_flow = '//bool2str(h2oice_flow))
call print_msg('co2ice_flow = '//bool2str(co2ice_flow))

END SUBROUTINE set_glaciers_config
!=======================================================================

!=======================================================================
SUBROUTINE flow_co2glaciers(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,co2ice,is_co2ice_flow)
!-----------------------------------------------------------------------
! NAME
!     flow_co2glaciers
!
! DESCRIPTION
!     Compute maximum thickness and transfer CO2 ice between subslopes.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope
use display,  only: print_msg

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:), intent(in)    :: vmr_co2_PEM     ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
real(dp),    dimension(:,:), intent(in)    :: ps_PCM          ! Grid points x Time field: surface pressure given by the PCM [Pa]
real(dp),                    intent(in)    :: ps_avg_glob_ini ! Global averaged surface pressure at the beginning [Pa]
real(dp),                    intent(in)    :: ps_avg_global   ! Global averaged surface pressure during the PEM iteration [Pa]
real(dp),    dimension(:,:), intent(inout) :: co2ice          ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
logical(k4), dimension(:,:), intent(out)   :: is_co2ice_flow  ! Flag to see if there is flow on the subgrid slopes

! LOCAL VARIABLES
! ---------------
real(dp), dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K]
real(dp), dimension(ngrid,nslope) :: hmax  ! Maximum thickness before flow

! CODE
! ----
call print_msg("> Flow of CO2 glaciers")
call computeTcondCO2(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,Tcond)
call compute_hmaxglaciers(Tcond,"CO2",hmax)
call transfer_ice_duringflow(hmax,Tcond,co2ice,is_co2ice_flow)

END SUBROUTINE flow_co2glaciers
!=======================================================================

!=======================================================================
SUBROUTINE flow_h2oglaciers(Tice,h2oice,is_h2oice_flow)
!-----------------------------------------------------------------------
! NAME
!     flow_h2oglaciers
!
! DESCRIPTION
!     Compute maximum thickness and transfer H2O ice between subslopes.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope
use display,  only: print_msg

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:), intent(in)    ::  Tice           ! Ice Temperature [K]
real(dp),    dimension(:,:), intent(inout) ::  h2oice         ! H2O ice on the subgrid slopes [kg/m^2]
logical(k4), dimension(:,:), intent(out)   ::  is_h2oice_flow ! Flow flag on subgrid slopes

! LOCAL VARIABLES
! ---------------
real(dp), dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow

! CODE
! ----
call print_msg("> Flow of H2O glaciers")
call compute_hmaxglaciers(Tice,"H2O",hmax)
call transfer_ice_duringflow(hmax,Tice,h2oice,is_h2oice_flow)

END SUBROUTINE flow_h2oglaciers
!=======================================================================

!=======================================================================
SUBROUTINE compute_hmaxglaciers(Tice,name_ice,hmax)
!-----------------------------------------------------------------------
! NAME
!     compute_hmaxglaciers
!
! DESCRIPTION
!     Compute the maximum thickness of CO2 and H2O glaciers given a
!     slope angle before initiating flow.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,  only: ngrid, nslope
use slopes,    only: iflat, def_slope_mean
use ice_table, only: rho_ice
use stoppage,  only: stop_clean
use maths,     only: pi
use physics,   only: g

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: Tice     ! Ice temperature [K]
character(3),             intent(in)  :: name_ice ! Nature of ice
real(dp), dimension(:,:), intent(out) :: hmax     ! Maximum thickness before flow [m]

! LOCAL VARIABLES
! ---------------
real(dp)    :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
integer(di) :: ig, islope
real(dp)    :: slo_angle

! CODE
! ----
select case (trim(adjustl(name_ice)))
    case('H2O')
        tau_d = 1.e5_dp
    case('CO2')
        tau_d = 5.e3_dp
    case default
        call stop_clean(__FILE__,__LINE__,"type of ice unknown!",1)
end select

do ig = 1,ngrid
    do islope = 1,nslope
        if (islope == iflat) then
            hmax(ig,islope) = 1.e8_dp
        else
            slo_angle = abs(def_slope_mean(islope)*pi/180._dp)
            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
        end if
    end do
end do

END SUBROUTINE compute_hmaxglaciers
!=======================================================================

!=======================================================================
SUBROUTINE transfer_ice_duringflow(hmax,Tice,qice,flag_flow)
!-----------------------------------------------------------------------
! NAME
!     transfer_ice_duringflow
!
! DESCRIPTION
!     Transfer the excess of ice from one subslope to another.
!     No transfer between mesh at the time.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,  only: ngrid, nslope
use slopes,    only: iflat, subslope_dist, def_slope_mean
use ice_table, only: rho_ice
use stoppage,  only: stop_clean
use maths,     only: pi

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:), intent(in)    :: hmax      ! Maximum height before initiating flow [m]
real(dp),    dimension(:,:), intent(in)    :: Tice      ! Ice temperature [K]
real(dp),    dimension(:,:), intent(inout) :: qice      ! Ice in the subslope [kg/m^2]
logical(k4), dimension(:,:), intent(out)   :: flag_flow ! Flow flag on subgrid slopes

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, islope
integer(di) :: iaval ! Index where ice is transferred

! CODE
! ----
flag_flow = .false.

do ig = 1,ngrid
    do islope = 1,nslope
        if (islope /= iflat) then ! ice can be infinite on flat ground
! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) then
! Second: determine the flatest slopes possible:
                if (islope > iflat) then
                    iaval = islope - 1
                else
                    iaval = islope + 1
                end if
                do while (iaval /= iflat .and. abs(subslope_dist(ig,iaval)) < minieps)
                    if (iaval > iflat) then
                        iaval = iaval - 1
                    else
                        iaval = iaval + 1
                    end if
                end do
                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) &
                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180._dp)/cos(pi*def_slope_mean(islope)/180._dp)

                qice(ig,islope) = rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)
                flag_flow(ig,islope) = .true.
            end if ! co2ice > hmax
        end if ! iflat
    end do !islope
end do !ig

END SUBROUTINE

!=======================================================================
SUBROUTINE computeTcondCO2(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,Tcond)
!-----------------------------------------------------------------------
! NAME
!     computeTcondCO2
!
! DESCRIPTION
!     Compute CO2 condensation temperature.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nday
use planet,   only: alpha_clap_co2, beta_clap_co2

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: vmr_co2_PEM       ! VMR of CO2 in the first layer [mol/mol]
real(dp), dimension(:,:), intent(in)  :: ps_PCM            ! Surface pressure in the PCM [Pa]
real(dp),                 intent(in)  :: ps_avg_glob_ini ! Global averaged surfacepressure in the PCM [Pa]
real(dp),                 intent(in)  :: ps_avg_global     ! Global averaged surface pressure computed during the PEM iteration
real(dp), dimension(:,:), intent(out) :: Tcond             ! Condensation temperature of CO2, yearly averaged

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, it

! CODE
! ----
do ig = 1,ngrid
    Tcond(ig,:) = 0._dp
    do it = 1,nday
        Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_glob_ini/ps_avg_global/100._dp))
    end do
    Tcond(ig,:) = Tcond(ig,:)/nday
end do

END SUBROUTINE computeTcondCO2
!=======================================================================

END MODULE glaciers