MODULE ice_sheet
!-----------------------------------------------------------------------
! NAME
!     ice_sheet
!
! DESCRIPTION
!     Compute the flow of the H2O ice caps/sheets between the grid cells.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     The ice sitting on a cell can flow towards its eight neighbouring
!     cells taking into account:
!       - the true spherical geometry of the planet;
!       - the topography of the ice surface, where the bed elevation is
!         derived from the surface geopotential;
!       - the temperature of the ice, described by a steady conduction
!         profile;
!       - the shallow-ice approximation (SIA) of Glen's flow law for the
!         internal deformation, with a two-branch Arrhenius law for
!         the ice softness. The ice flows whatever the stress but the rate
!         is set by the ice thickness, the surface slope and the ice
!         temperature;
!       - a Weertman sliding law for basal sliding of a temperate bed, with
!         a reference speed 'u_slide_ref';
!       - a sub-grid ice margin to represent the ice edge, that is smoothed
!         by the grid resolution. A plastic ice margin with a yield stress
!         'tau_d' is used to set the steepness of the ice edge. It is graded
!         by the thickness ratio so that the resolved slope is used in the
!         interior of the ice sheet.
!
!     The levelling bound and the caps on the outflow of each cell make
!     the scheme unconditionally stable.
!     Mass is conserved exactly.
!
!     The flow can bring ice onto a previously ice-free cell, which may
!     cause the PEM run to stop.
!
!     Some known limitations:
!       - basal melting is not represented: where the bed is temperate,
!         the heat only switches the sliding on;
!       - the sliding law has no dependence on the effective pressure
!         nor on any subglacial hydrology;
!       - the ice margin is a sub-grid parameterization, not a resolved
!         slope;
!       - the 8-direction cubic splitting is not a consistent
!         discretization of the SIA equation because the slope is cubed
!         in each direction separately rather than once for the whole
!         gradient;
!       - there is no isostatic adjustment of the bed under the ice load;
!       - there is no thermal advection of the ice temperature with the
!         flow, so the ice temperature is always a conduction profile;
!       - there is no calving of ice into the ocean or lakes, so ice
!         simply grounds there;
!       - SIA approximations hold: no longitudinal stresses, i.e. no ice
!         streams, no ice shelves, no ice tongues, no ice rises, no ice
!         rumples;
!       - 'gamma' and 'hice' are taken at the cell centre rather than
!         reconstructed at the interface.
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
! Rheology and flow parameters
integer(di), parameter, private :: n_glen     = 3_di         ! Stress exponent of Glen's flow law, Glen 1955
real(dp),    parameter, private :: A_cold     = 3.985e-13_dp ! Arrhenius pre-exponential factor below 'T_arrh' [Pa-3.s-1], Greve & Blatter 2009
real(dp),    parameter, private :: A_warm     = 1.916e3_dp   ! Arrhenius pre-exponential factor above 'T_arrh' [Pa-3.s-1], Greve & Blatter 2009
real(dp),    parameter, private :: Q_cold     = 6.0e4_dp     ! Activation energy below 'T_arrh' [J/mol], Greve & Blatter 2009
real(dp),    parameter, private :: Q_warm     = 1.39e5_dp    ! Activation energy above 'T_arrh' [J/mol], Greve & Blatter 2009
real(dp),    parameter, private :: T_arrh     = 263.15_dp    ! Temperature of the Arrhenius branch change [K], Greve & Blatter 2009
real(dp),    parameter, private :: T_melt     = 273.15_dp    ! Upper bound for the ice temperature [K]
real(dp),    parameter, private :: T_floor    = 50._dp       ! Lower bound for the ice temperature [K]
real(dp),    parameter, private :: h_flow_min = 1._dp        ! Minimum ice thickness to allow flow [m]

! Ice margin and basal sliding parameters
real(dp), parameter, private :: tau_d       = 1.e5_dp ! Yield stress setting the ice margin slope [Pa], Cuffey & Paterson 2010
real(dp), parameter, private :: u_slide_ref = 100._dp ! Sliding speed of a temperate ice margin [m/plnt-yr], chosen value

! Thermal parameters of the ice column
real(dp),    parameter, private :: A_cond   = 651._dp    ! Conductivity coefficient [W/m] so that thermal conductivity k(T) = A_cond/T, Petrenko & Whitworth 1999
real(dp),    parameter, private :: T_triple = 273.16_dp  ! Triple-point temperature of H2O [K]
real(dp),    parameter, private :: p_triple = 611.657_dp ! Triple-point pressure of H2O [Pa]
real(dp),    parameter, private :: T_Ih_min = 251.165_dp ! Melting temperature of ice Ih at the triple point with ice III [K]
real(dp),    parameter, private :: a1_melt  = 0.626e6_dp ! Coefficients of the melting curve of ice Ih, Wagner et al. 1994
real(dp),    parameter, private :: a2_melt  = 0.197135e6_dp
real(dp),    parameter, private :: b2_melt  = 21.2_dp
integer(di), parameter, private :: n_quad   = 32_di      ! Number of intervals (even) across the depth for the integration of the softness

! VARIABLES
! ---------
logical(k4),                         protected :: h2oice_flow = .true. ! Flag to compute the flow of the H2O ice sheets
real(dp), dimension(:), allocatable, protected :: zbed                 ! Bed elevation [m]
real(dp), dimension(:), allocatable, protected :: d_h2oice_flow        ! H2O ice tendency due to the ice sheet flow [kg/m2/y]

! Connectivity of the grid cells with their neighbours
integer(di), dimension(:),   allocatable, protected, private :: n_nghb               ! Number of neighbours of each cell [-]
integer(di), dimension(:,:), allocatable, protected, private :: i_nghb               ! Indices of the neighbours [-]
real(dp),    dimension(:,:), allocatable, protected, private :: d_nghb               ! Great-circle distance between the cell centres [m]
real(dp),    dimension(:,:), allocatable, protected, private :: w_nghb               ! Effective width of the interface [m]
logical(k4),                              protected, private :: nghb_built = .false. ! Flag to true when the connectivity is built

! VISIBILITY
! ----------
private :: build_cell_neighbours, add_neighbour, great_circle
private :: arrhenius, t_melt_ice, column_softness

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

!=======================================================================
SUBROUTINE ini_ice_sheet()
!-----------------------------------------------------------------------
! NAME
!     ini_ice_sheet
!
! DESCRIPTION
!     Allocate 'ice_sheet' module arrays.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nlon

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

! LOCAL VARIABLES
! ---------------
integer(di) :: n_nghb_max ! Maximum number of neighbours of a cell

! CODE
! ----
! The polar cells have 'nlon' neighbours, the regular ones have 8
n_nghb_max = max(8_di,nlon)

allocate(zbed(ngrid))
allocate(d_h2oice_flow(ngrid))
allocate(n_nghb(ngrid))
allocate(i_nghb(ngrid,n_nghb_max))
allocate(d_nghb(ngrid,n_nghb_max))
allocate(w_nghb(ngrid,n_nghb_max))

zbed(:) = 0._dp
d_h2oice_flow(:) = 0._dp
n_nghb(:) = 0_di
i_nghb(:,:) = 0_di
d_nghb(:,:) = 0._dp
w_nghb(:,:) = 0._dp
nghb_built = .false.

END SUBROUTINE ini_ice_sheet
!=======================================================================

!=======================================================================
SUBROUTINE end_ice_sheet()
!-----------------------------------------------------------------------
! NAME
!     end_ice_sheet
!
! DESCRIPTION
!     Deallocate 'ice_sheet' module arrays.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
if (allocated(zbed)) deallocate(zbed)
if (allocated(d_h2oice_flow)) deallocate(d_h2oice_flow)
if (allocated(n_nghb)) deallocate(n_nghb)
if (allocated(i_nghb)) deallocate(i_nghb)
if (allocated(d_nghb)) deallocate(d_nghb)
if (allocated(w_nghb)) deallocate(w_nghb)
nghb_built = .false.

END SUBROUTINE end_ice_sheet
!=======================================================================

!=======================================================================
SUBROUTINE set_ice_sheet_config(h2oice_flow_in)
!-----------------------------------------------------------------------
! NAME
!     set_ice_sheet_config
!
! DESCRIPTION
!     Setter for "ice_sheet" configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility,  only: bool2str
use geometry, only: ngrid
use display,  only: print_msg, LVL_NFO, LVL_WRN

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

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

! CODE
! ----
h2oice_flow = h2oice_flow_in
if (h2oice_flow .and. ngrid == 1) then
    call print_msg('''h2oice_flow = .true.'' but ngrid = 1. So there will be no flow!',LVL_WRN)
    h2oice_flow = .false.
end if
call print_msg('h2oice_flow = '//bool2str(h2oice_flow),LVL_NFO)

END SUBROUTINE set_ice_sheet_config
!=======================================================================

!=======================================================================
SUBROUTINE set_zbed(phisfi)
!-----------------------------------------------------------------------
! NAME
!     set_zbed
!
! DESCRIPTION
!     Setter for 'zbed' from the surface geopotential, then build the
!     connectivity of the grid cells with their neighbours.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     Must be called once "startfi.nc" has been read, so that the planet
!     radius and the gravity are known.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid
use physics,  only: g
use utility,  only: real2str
use display,  only: print_msg, LVL_NFO
use stoppage, only: stop_clean

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

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in) :: phisfi ! Surface geopotential [m2/s2]

! CODE
! ----
if (size(phisfi) /= ngrid) call stop_clean(__FILE__,__LINE__,'problem of dimensions!',1)
if (g <= 0._dp) call stop_clean(__FILE__,__LINE__,'gravity is not initialized. Call ''init_physics'' first.',1)

zbed(:) = phisfi(:)/g
call print_msg('Bed elevation [m] within ['//real2str(minval(zbed))//','//real2str(maxval(zbed))//']',LVL_NFO)

call build_cell_neighbours()

END SUBROUTINE set_zbed
!=======================================================================

!=======================================================================
SUBROUTINE build_cell_neighbours()
!-----------------------------------------------------------------------
! NAME
!     build_cell_neighbours
!
! DESCRIPTION
!     Build the connectivity of each grid cell with its neighbours.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     The effective widths generalize the multiple-flow-direction contour
!     lengths of Quinn et al. (1991) (0.5*dx for a cardinal direction and
!     dx/(2*sqrt(2)) for a diagonal one) to a lon-lat grid on a sphere:
!       - cardinal north/south: half of the shared latitudinal edge;
!       - cardinal east/west  : half of the shared longitudinal edge;
!       - diagonal            : (L_lat + L_lon)/(4*sqrt(2)).
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nlon, nlat, latitudes, longitudes
use physics,  only: rad
use maths,    only: pi
use utility,  only: int2str
use display,  only: print_msg, LVL_NFO
use stoppage, only: stop_clean

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

! LOCAL VARIABLES
! ---------------
integer(di)                         :: ig, jg, kg, ilon, ilat, jlon, jlat, d_ilon, d_ilat, k, kk
real(dp), dimension(:), allocatable :: lat_row  ! Latitude of the centre of each row [rad]
real(dp), dimension(:), allocatable :: lat_edge ! Latitude of the interface between two rows [rad]
real(dp), dimension(:), allocatable :: len_lat  ! Length of the latitudinal edges [m]
real(dp), dimension(:), allocatable :: len_lon  ! Length of the longitudinal edges [m]
real(dp)                            :: dlon     ! Longitude increment [rad]
real(dp)                            :: w        ! Width of an interface [m]

! CODE
! ----
if (ngrid == 1) then ! Nothing to connect in 1D
    nghb_built = .true.
    return
end if
if (rad <= 0._dp) call stop_clean(__FILE__,__LINE__,'planet radius is not initialized!',1)
if (nlat < 3) call stop_clean(__FILE__,__LINE__,'the grid must have at least 3 latitudes to build the connectivity!',1)

call print_msg('> Building the connectivity of the grid cells for the H2O ice sheet flow',LVL_NFO)

! Latitude of the centre of each row
allocate(lat_row(nlat),lat_edge(nlat - 1),len_lat(nlat - 1),len_lon(nlat))
lat_row(1) = latitudes(1)
lat_row(nlat) = latitudes(ngrid)
do ilat = 2,nlat - 1
    lat_row(ilat) = latitudes(2 + (ilat - 2)*nlon)
end do

! Latitude of the interface between two rows
do ilat = 1,nlat - 1
    lat_edge(ilat) = 0.5_dp*(lat_row(ilat) + lat_row(ilat + 1))
end do

! Length of the edges: a cell on the sphere is a trapezoid
dlon = 2._dp*pi/real(nlon,dp)
do ilat = 1,nlat - 1
    len_lat(ilat) = rad*cos(lat_edge(ilat))*dlon ! The parallels shrink towards the poles
end do
len_lon(:) = 0._dp
do ilat = 2,nlat - 1
    len_lon(ilat) = rad*(lat_edge(ilat - 1) - lat_edge(ilat))
end do

! Neighbours of the regular cells
do ilat = 2,nlat - 1
    do ilon = 1,nlon
        ig = 1 + ilon + (ilat - 2)*nlon
        do d_ilat = -1,1
            do d_ilon = -1,1
                if (d_ilat == 0 .and. d_ilon == 0) cycle

                ! Index of the neighbour
                jlat = ilat + d_ilat
                jlon = ilon + d_ilon
                if (jlon < 1) jlon = jlon + nlon ! Longitude periodicity
                if (jlon > nlon) jlon = jlon - nlon ! Longitude periodicity
                if (jlat == 1) then ! North pole
                    jg = 1
                else if (jlat == nlat) then ! South pole
                    jg = ngrid
                else
                    jg = 1 + jlon + (jlat - 2)*nlon
                end if
                if (jg == ig) cycle ! Safeguard

                ! Width of the interface
                if (d_ilat == 0) then ! Cardinal east/west
                    w = 0.5_dp*len_lon(ilat)
                else if (d_ilat == -1) then ! Northward
                    if (d_ilon == 0) then ! Cardinal north
                        w = 0.5_dp*len_lat(ilat - 1)
                    else ! Diagonal north
                        w = (len_lat(ilat - 1) + len_lon(ilat))/(4._dp*sqrt(2._dp))
                    end if
                else ! Southward
                    if (d_ilon == 0) then ! Cardinal south
                        w = 0.5_dp*len_lat(ilat)
                    else ! Diagonal south
                        w = (len_lat(ilat) + len_lon(ilat))/(4._dp*sqrt(2._dp))
                    end if
                end if

                call add_neighbour(ig,jg,w)
            end do
        end do
    end do
end do

! Neighbours of the polar cells
do ilon = 1,nlon
    call add_neighbour(1_di,1 + ilon,len_lat(1)) ! North pole
    call add_neighbour(ngrid,1 + ilon + (nlat - 3)*nlon,len_lat(nlat - 1)) ! South pole
end do

! Great-circle distance between the cell centres
do ig = 1,ngrid
    do k = 1,n_nghb(ig)
        jg = i_nghb(ig,k)
        d_nghb(ig,k) = great_circle(latitudes(ig),longitudes(ig),latitudes(jg),longitudes(jg))
        if (d_nghb(ig,k) <= 0._dp) call stop_clean(__FILE__,__LINE__,'null distance between two neighbouring cells!',1)
    end do
end do

! Symmetrize the widths so that the flux does not depend on the direction
do ig = 1,ngrid
    do k = 1,n_nghb(ig)
        jg = i_nghb(ig,k)
        if (jg <= ig) cycle ! Already done
        kk = 0_di
        do kg = 1,n_nghb(jg)
            if (i_nghb(jg,kg) == ig) then ! Found the reverse neighbour
                kk = kg
                exit
            end if
        end do
        if (kk == 0_di) call stop_clean(__FILE__,__LINE__,'the connectivity of the grid cells is not symmetric!',1)
        w = 0.5_dp*(w_nghb(ig,k) + w_nghb(jg,kk))
        w_nghb(ig,k) = w
        w_nghb(jg,kk) = w
    end do
end do

deallocate(lat_row,lat_edge,len_lat,len_lon)
nghb_built = .true.

END SUBROUTINE build_cell_neighbours
!=======================================================================

!=======================================================================
SUBROUTINE add_neighbour(ig,jg,w)
!-----------------------------------------------------------------------
! NAME
!     add_neighbour
!
! DESCRIPTION
!     Add the cell 'jg' to the list of the neighbours of the cell 'ig'
!     with the interface width 'w'.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     If 'jg' is already a neighbour of 'ig', the width is accumulated
!     instead of appending a new entry.
!     The accumulation allows to handle the case of cells next to a pole
!     whose three northward/southward neighbours are actually the same
!     polar cell.
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
integer(di), intent(in) :: ig, jg ! Indices of the cell and its neighbour
real(dp),    intent(in) :: w      ! Width of the interface [m]

! LOCAL VARIABLES
! ---------------
integer(di) :: k

! CODE
! ----
do k = 1,n_nghb(ig)
    if (i_nghb(ig,k) == jg) then
        w_nghb(ig,k) = w_nghb(ig,k) + w
        return
    end if
end do

if (n_nghb(ig) >= size(i_nghb,2)) call stop_clean(__FILE__,__LINE__,'too many neighbours for a grid cell!',1)
n_nghb(ig) = n_nghb(ig) + 1
i_nghb(ig,n_nghb(ig)) = jg
w_nghb(ig,n_nghb(ig)) = w

END SUBROUTINE add_neighbour
!=======================================================================

!=======================================================================
FUNCTION great_circle(lat1,lon1,lat2,lon2) RESULT(d)
!-----------------------------------------------------------------------
! NAME
!     great_circle
!
! DESCRIPTION
!     Great-circle distance between two points of the sphere surface.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     Use of the Haversine formula.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use physics, only: rad

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: lat1, lon1, lat2, lon2 ! Coordinates of the two points [rad]

! RESULT
! ------
real(dp) :: d ! Great-circle distance [m]

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

! CODE
! ----
a = sin(0.5_dp*(lat2 - lat1))**2 + cos(lat1)*cos(lat2)*sin(0.5_dp*(lon2 - lon1))**2
d = 2._dp*rad*asin(min(sqrt(max(0._dp,a)),1._dp))

END FUNCTION great_circle
!=======================================================================

!=======================================================================
FUNCTION arrhenius(Tice) RESULT(A)
!-----------------------------------------------------------------------
! NAME
!     arrhenius
!
! DESCRIPTION
!     Compute the softness of the H2O ice.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     Two-branch Arrhenius law for the rate factor (softness) in Glen's
!     flow law (Paterson & Budd 1982).
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use physics, only: Rgas

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: Tice ! Ice temperature [K]

! RESULT
! ------
real(dp) :: A ! Softness of the ice [Pa-3.s-1]

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

! CODE
! ----
T = min(max(T_floor,Tice),T_melt)
if (T <= T_arrh) then
    A = A_cold*exp(-Q_cold/(Rgas*T))
else
    A = A_warm*exp(-Q_warm/(Rgas*T))
end if

END FUNCTION arrhenius
!=======================================================================

!=======================================================================
FUNCTION t_melt_ice(p) RESULT(T_melt_point)
!-----------------------------------------------------------------------
! NAME
!     t_melt_ice
!
! DESCRIPTION
!     Melting temperature of ice Ih under a given pressure.
!
! AUTHORS & DATE
!     JB Clement, 09/2026
!
! NOTES
!     Newton method to invert the melting curve of Wagner et al. 1994,
!     bounded by the triple points with the vapour and with ice III:
!       p_melt(T) = p_triple*(1 − a1*(1 − 1/theta^3) + a2*(1 − theta^b2)),
!           where theta = T/T_triple.
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: p ! Pressure [Pa]

! RESULT
! ------
real(dp) :: T_melt_point ! Melting temperature [K]

! LOCAL VARIABLES
! ---------------
integer(di) :: it
real(dp)    :: theta, dpdT, dT_melt

! CODE
! ----
T_melt_point = T_triple
if (p <= p_triple) return

! Initial guess from the slope of the curve at the triple point
T_melt_point = max(T_triple - (p - p_triple)*T_triple/(p_triple*(3._dp*a1_melt + b2_melt*a2_melt)),T_Ih_min)

! Newton method
do it = 1,20
    theta = T_melt_point/T_triple
    dpdT = -p_triple/T_triple*(3._dp*a1_melt/theta**4 + b2_melt*a2_melt*theta**(b2_melt - 1._dp))
    dT_melt = -(p_triple*(1._dp - a1_melt*(1._dp - 1._dp/theta**3) + a2_melt*(1._dp - theta**b2_melt)) - p)/dpdT
    dT_melt = min(max(T_Ih_min,T_melt_point + dT_melt),T_triple) - T_melt_point ! Keep the solution within the bounds of the melting curve
    T_melt_point = T_melt_point + dT_melt
    if (abs(dT_melt) < 1.e-9_dp) exit
end do

END FUNCTION t_melt_ice
!=======================================================================

!=======================================================================
FUNCTION column_softness(tsurf,hice,psurf) RESULT(A_col)
!-----------------------------------------------------------------------
! NAME
!     column_softness
!
! DESCRIPTION
!     Softness of Glen's flow law averaged over the ice column with the
!     weight of the shallow-ice flux.
!
! AUTHORS & DATE
!     JB Clement, 09/2026
!
! NOTES
!     Integration of Glen's law over the column is given in Hutter 1983
!     and Greve & Blatter 2009. The weight x^(n_glen + 1) favors the base
!     of the ice column where the flux is the largest.
!     The composite Simpson rule is used to integrate the softness over
!     the depth. It is normalized so that an isothermal column gives back
!     its own softness exactly.
!
!     Conduction profile with a uniform upward flux (Abbot & Switzer 2011):
!       T(z) = Tsurf*exp(flux_geo_deep*z/A_cond)
!
!     The Arrhenius law is fitted on ice near atmospheric pressure, where
!     the melting point is about 273 K. The softness depends on the
!     temperature relative to the local melting point rather than on the
!     absolute temperature (Greve & Blatter 2009, Cuffey & Paterson 2010).
!     Hence the rescaling.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use surf_ice, only: rho_h2oice
use physics,  only: g
use soil,     only: flux_geo_deep

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: tsurf ! Surface temperature [K]
real(dp), intent(in) :: hice  ! Ice thickness [m]
real(dp), intent(in) :: psurf ! Surface pressure [Pa]

! RESULT
! ------
real(dp) :: A_col ! Softness of the ice column [Pa-3.s-1]

! LOCAL VARIABLES
! ---------------
integer(di) :: k
real(dp)    :: x, w, wsum, T, T_melt_point

! CODE
! ----
A_col = 0._dp
wsum = 0._dp
do k = 0,n_quad
    x = real(k,dp)/real(n_quad,dp)
    if (k == 0 .or. k == n_quad) then
        w = 1._dp
    else if (mod(k,2_di) == 1) then
        w = 4._dp
    else
        w = 2._dp
    end if
    w = w*x**(n_glen + 1) ! Weight of the shallow-ice flux to each depth

    ! Conduction profile, rising with depth, capped at the local melting point
    T_melt_point = t_melt_ice(psurf + rho_h2oice*g*x*hice)
    T = min(max(tsurf,T_floor)*exp(flux_geo_deep*x*hice/A_cond),T_melt_point)
    A_col = A_col + w*arrhenius(T + T_triple - T_melt_point) ! Softness relative to the melting point
    wsum = wsum + w
end do
A_col = A_col/wsum

END FUNCTION column_softness
!=======================================================================

!=======================================================================
SUBROUTINE flow_h2oice(h2o_ice,tsurf,psurf,zshift_surf)
!-----------------------------------------------------------------------
! NAME
!     flow_h2oice
!
! DESCRIPTION
!     Flow the H2O ice from a cell towards its eight neighbouring cells.
!
! AUTHORS & DATE
!     JB Clement, 10/2026
!
! NOTES
!     All the fluxes are computed from the state at the beginning of the
!     step, so the result does not depend on the order of the cells.
!
!     Weertman's sliding law says the sliding speed of ice over its bed
!     is a power of the basal shear stress (Weertman 1957):
!       u_b = A_s*tau_b^n
!             with A_s the sliding parameter [m/s/Pa^n],
!                  tau_b the basal shear stress [Pa],
!                  n = 3 the stress exponent of Glen's flow law [-].
!
!     The shallow-ice approximation (SIA) of the Glen flow law says the
!     deformation flux of ice is a power of the surface slope (Hutter 1983,
!     Morland & Johnson 1980):
!       q = -gamma*h^(n+2)*|grad(s)|^(n - 1)*grad(s)
!             with gamma = 2*A*(rho*g)^n/(n + 2) the prefactor of the
!                          Glen/SIA flux [m-3.s-1],
!                  h the ice thickness [m],
!                  n = 3 the stress exponent of Glen's flow law [-],
!                  s the ice surface elevation [m], slope(s) = |grad(s)|,
!                  A the softness of the ice [Pa-3.s-1],
!                  rho the ice density [kg/m3],
!                  g the gravity [m/s2].
!     1. Ice flows down the slope of its surface, not of its bed.
!     2. The law is stiff. For n = 3, the flux grows as h^5*s^3.
!     3. The flux scales as (rho*g)^3, so it is very sensitive to the
!        planet gravity and the ice density.
!     4. Ice creeps at any stress, slowly when thin, flat or cold (no threshold).
!     5. The approximation degrades near margins and over steep or rough
!        beds, where longitudinal stresses are no longer negligible.
!
!     Conduction profile with a uniform upward flux (Abbot & Switzer 2011):
!       T(z) = Tsurf*exp(flux_geo_deep*z/A_cond)
!
!     The slope margin is a sub-grid correction to the slope of the ice
!     surface. If ice is a perfectly plastic material, it does not deform
!     below a yield stress tau_d (Nye 1951, 1952):
!       S_margin = tau_d/(rho*g*h)
!     This law gives a parabolic profile of the ice surface near the margin
!     of an ice cap.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,  only: ngrid, nslope, cell_area
use slopes,    only: slope_weight
use surf_ice,  only: rho_h2oice
use physics,   only: g
use soil,      only: flux_geo_deep
use orbit,     only: yr_len_sol, sol_len_s
use evolution, only: dt, get_dt_yr
use display,   only: print_msg, LVL_NFO
use stoppage,  only: stop_clean

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(inout) :: h2o_ice     ! H2O ice [kg/m2]
real(dp), dimension(:,:), intent(in)    :: tsurf       ! Surface temperature [K]
real(dp), dimension(:),   intent(in)    :: psurf       ! Surface pressure [Pa]
real(dp), dimension(:,:), intent(inout) :: zshift_surf ! Elevation shift for the surface [m]

! LOCAL VARIABLES
! ---------------
integer(di)                         :: ig, jg, islope, k
real(dp), dimension(ngrid)          :: mass_ice, mass_ice_new ! Column mass of H2O ice [kg/m2]
real(dp), dimension(ngrid)          :: hice         ! Ice thickness [m]
real(dp), dimension(ngrid)          :: zice         ! Ice surface elevation [m]
real(dp), dimension(ngrid)          :: Tice         ! Ice temperature [K]
real(dp), dimension(ngrid)          :: dV           ! Volume of ice gained by each cell [m3]
real(dp), dimension(ngrid)          :: wslope_tot   ! Total weight of the sub-slopes [-]
real(dp), dimension(size(i_nghb,2)) :: V_nghb       ! Volume of ice given to each neighbour [m3]
real(dp)                            :: gamma        ! Prefactor of the Glen/SIA flux [m-3.s-1]
real(dp)                            :: slope_nghb   ! Slope of the ice surface towards a neighbour [-]
real(dp)                            :: edge_ice     ! Amount of ice edge at an interface [-]
real(dp)                            :: slope_eff    ! Slope of the ice surface corrected for the sub-grid margin [-]
real(dp)                            :: tau_b        ! Basal shear stress [Pa]
real(dp)                            :: A_s          ! Sliding parameter of the Weertman sliding law [m/s/Pa3]
real(dp)                            :: T_bed        ! Temperature at the base of the ice [K]
real(dp)                            :: T_melt_point ! Pressure-melting point at the base of the ice [K]
logical(k4)                         :: temperate    ! Flag to true if the bed is at the melting point
real(dp)                            :: V_glen       ! Volume of ice moved by the Glen/SIA flux [m3]
real(dp)                            :: V_sliding    ! Volume of ice moved by the basal sliding [m3]
real(dp)                            :: V_lvl        ! Volume of ice levelling the two ice surfaces [m3]
real(dp)                            :: V_sum        ! Total volume of ice leaving a cell [m3]
real(dp)                            :: Vmax         ! Upper bound for the volume of ice leaving a cell [m3]
real(dp)                            :: zice_max     ! Highest ice surface among the downhill neighbours [m]
real(dp)                            :: dt_yr, coef

! CODE
! ----
call print_msg('> Flow of the H2O ice sheets',LVL_NFO)
if (.not. nghb_built) call stop_clean(__FILE__,__LINE__,'the connectivity of the grid cells is not built!',1)

d_h2oice_flow(:) = 0._dp
if (ngrid == 1) return ! Nothing to flow in 1D

! Aggregate the sub-slopes into a cell mean
do ig = 1,ngrid
    mass_ice(ig) = 0._dp
    Tice(ig) = 0._dp
    wslope_tot(ig) = 0._dp
    do islope = 1,nslope
        mass_ice(ig) = mass_ice(ig) + h2o_ice(ig,islope)*slope_weight(ig,islope)
        Tice(ig) = Tice(ig) + tsurf(ig,islope)*slope_weight(ig,islope)
        wslope_tot(ig) = wslope_tot(ig) + slope_weight(ig,islope)
    end do
    if (wslope_tot(ig) <= 0._dp) call stop_clean(__FILE__,__LINE__,'null total weight of the sub-slopes of a grid cell!',1)
    Tice(ig) = Tice(ig)/wslope_tot(ig)
    hice(ig) = mass_ice(ig)/rho_h2oice
    zice(ig) = zbed(ig) + hice(ig)
end do

! Sliding parameter of the Weertman sliding law
! Set so that a temperate ice margin (where the basal shear stress equals 'tau_d' by construction) slides at 'u_slide_ref' whatever the ice thickness
A_s = u_slide_ref/(yr_len_sol*sol_len_s)/tau_d**n_glen

! Compute the volumes of ice exchanged between the cells
dV(:) = 0._dp
do ig = 1,ngrid
    if (hice(ig) < h_flow_min) cycle ! No ice to flow

    ! Prefactor of the Glen/SIA flux, averaged over the ice column
    gamma = 2._dp*column_softness(Tice(ig),hice(ig),psurf(ig))*(rho_h2oice*g)**n_glen/real(n_glen + 2,dp)

    ! Is the bed at the pressure-melting point? Cold-based ice does not slide
    T_bed = max(Tice(ig),T_floor)*exp(flux_geo_deep*hice(ig)/A_cond) ! Conduction profile
    T_melt_point = t_melt_ice(psurf(ig) + rho_h2oice*g*hice(ig))
    temperate = (T_bed >= T_melt_point)

    ! Candidate volume given to each downhill neighbour
    V_nghb(:) = 0._dp
    V_sum = 0._dp
    zice_max = -huge(1._dp)
    do k = 1,n_nghb(ig)
        jg = i_nghb(ig,k)
        slope_nghb = (zice(ig) - zice(jg))/d_nghb(ig,k)
        if (slope_nghb <= 0._dp) cycle ! Uphill neighbour so no outflow

        ! Sub-grid margin slope: the grid cannot resolve an ice edge whose steepness is set by the yield stress
        ! Graded by the thickness ratio from 0 (neighbour as thick as the donor: interior, no correction) to 1 (bare neighbour: full margin slope)
        edge_ice = min(max(0._dp,(hice(ig) - hice(jg))/hice(ig)),1._dp)
        slope_eff = max(slope_nghb,edge_ice*tau_d/(rho_h2oice*g*hice(ig)))

        ! Internal deformation (Glen/SIA) and basal sliding (Weertman) over the time step
        V_glen = gamma*hice(ig)**(n_glen + 2)*slope_eff**n_glen*w_nghb(ig,k)*dt
        V_sliding = 0._dp
        if (temperate) then
            tau_b = rho_h2oice*g*hice(ig)*slope_eff ! Basal shear stress
            V_sliding = A_s*tau_b**n_glen*hice(ig)*w_nghb(ig,k)*dt
        end if

        ! Levelling volume that would make the two ice surfaces equal
        V_lvl = (zice(ig) - zice(jg))/(1._dp/cell_area(ig) + 1._dp/cell_area(jg))
        V_nghb(k) = min(V_glen + V_sliding,V_lvl) ! Enforce that the two ice surfaces do not invert
        V_sum = V_sum + V_nghb(k)
        zice_max = max(zice_max,zice(jg))
    end do
    if (V_sum <= 0._dp) cycle

    ! Never move more ice than the cell holds, nor drop its surface below the highest receiver
    Vmax = min(hice(ig),zice(ig) - zice_max)*cell_area(ig)
    if (Vmax <= 0._dp) cycle
    if (V_sum > Vmax) V_nghb(:) = V_nghb(:)*Vmax/V_sum

    ! Exchange the volumes
    do k = 1,n_nghb(ig)
        if (V_nghb(k) <= 0._dp) cycle
        jg = i_nghb(ig,k)
        dV(ig) = dV(ig) - V_nghb(k)
        dV(jg) = dV(jg) + V_nghb(k)
    end do
end do

! Back to the column mass
do ig = 1,ngrid
    mass_ice_new(ig) = mass_ice(ig) + dV(ig)*rho_h2oice/cell_area(ig)
    if (mass_ice_new(ig) < 0._dp) then
        if (mass_ice_new(ig) < -eps*max(1._dp,mass_ice(ig))) call stop_clean(__FILE__,__LINE__,'negative H2O ice produced by the ice sheet flow!',1)
        mass_ice_new(ig) = 0._dp ! Rounding error
    end if
end do

! Put the change back on the sub-slopes and shift the surface elevation
dt_yr = get_dt_yr()
do ig = 1,ngrid
    if (abs(mass_ice_new(ig) - mass_ice(ig)) < eps) cycle
    if (mass_ice_new(ig) < mass_ice(ig)) then
        ! Net loss: scale down the ice available, keeping the sub-slope repartition
        if (mass_ice(ig) <= 0._dp) cycle ! Safeguard
        coef = mass_ice_new(ig)/mass_ice(ig)
        do islope = 1,nslope
            zshift_surf(ig,islope) = zshift_surf(ig,islope) + (coef - 1._dp)*h2o_ice(ig,islope)/rho_h2oice
            h2o_ice(ig,islope) = coef*h2o_ice(ig,islope)
        end do
    else
        ! Net gain: spread uniformly over the true area of the cell
        coef = (mass_ice_new(ig) - mass_ice(ig))/wslope_tot(ig)
        do islope = 1,nslope
            zshift_surf(ig,islope) = zshift_surf(ig,islope) + coef/rho_h2oice
            h2o_ice(ig,islope) = h2o_ice(ig,islope) + coef
        end do
    end if
    d_h2oice_flow(ig) = (mass_ice_new(ig) - mass_ice(ig))/dt_yr
end do

END SUBROUTINE flow_h2oice
!=======================================================================

END MODULE ice_sheet