MODULE lake_ice
!-----------------------------------------------------------------------
! NAME
!     lake_ice
!
! DESCRIPTION
!     Module to manage lake ice.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
real(dp), protected, private :: k_snow = 0.31_dp ! Snow conductivity [W/m/K], PCM default from phygeneric/hydrol.F90

! Lake ice in the PCM at the beginning [kg/m2] and seasonal deviation to restore for PCM handoff
real(dp), dimension(:,:), allocatable, protected :: lakeice_PCM
real(dp), dimension(:,:), allocatable            :: lake_frost4PCM

! Fixed Generic PCM constants mirrored locally for PEM lake-ice evolution.
real(dp), parameter, private :: T_h2o_ice_liq = 273.16_dp  ! Temperature of water/ice interface [K], from phygeneric/comcstfi_mod.F90
real(dp), parameter, private :: RLFTT         = 3.334e5_dp ! Latent heat of freezing [J/kg], from phygeneric/surfdat_h.F90
real(dp), parameter, private :: rhowaterice   = 920._dp    ! Density of water ice [kg/m3], from phygeneric/hydrol.F90
real(dp), parameter, private :: rho_snow      = 300._dp    ! Density of snow [kg/m3], from phygeneric/hydrol.F90
real(dp), parameter, private :: k_ice         = 2.17_dp    ! Thermal conductivity of ice [W/m/K], from phygeneric/hydrol.F90

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

!=======================================================================
SUBROUTINE ini_lake_ice()
!-----------------------------------------------------------------------
! NAME
!     ini_lake_ice
!
! DESCRIPTION
!     Allocate 'lake_ice' module arrays.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
allocate(lakeice_PCM(ngrid,nslope))
allocate(lake_frost4PCM(ngrid,nslope))
lakeice_PCM(:,:) = 0._dp
lake_frost4PCM(:,:) = 0._dp

END SUBROUTINE ini_lake_ice
!=======================================================================

!=======================================================================
SUBROUTINE end_lake_ice()
!-----------------------------------------------------------------------
! NAME
!     end_lake_ice
!
! DESCRIPTION
!     Deallocate 'lake_ice' arrays.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
if (allocated(lakeice_PCM)) deallocate(lakeice_PCM)
if (allocated(lake_frost4PCM)) deallocate(lake_frost4PCM)

END SUBROUTINE end_lake_ice
!=======================================================================

!=======================================================================
SUBROUTINE set_k_snow(k_snow_in)
!-----------------------------------------------------------------------
! NAME
!     set_k_snow
!
! DESCRIPTION
!     Setter for 'k_snow'.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility,  only: real2str
use display,  only: print_msg, LVL_NFO
use stoppage, only: stop_clean

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

! ARGUMENTS
! ---------
real(dp), intent(in) :: k_snow_in

! CODE
! ----
k_snow = k_snow_in
call print_msg('k_snow = '//real2str(k_snow),LVL_NFO)
if (k_snow <= 0._dp) call stop_clean(__FILE__,__LINE__,'''k_snow'' must be positive!',1)

END SUBROUTINE set_k_snow
!=======================================================================

!=======================================================================
SUBROUTINE set_lakeice_PCM(lakeice_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_lakeice_PCM
!
! DESCRIPTION
!     Setter for 'lakeice_PCM'.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in) :: lakeice_PCM_in

! CODE
! ----
lakeice_PCM(:,:) = lakeice_PCM_in(:,:)

END SUBROUTINE set_lakeice_PCM
!=======================================================================

!=======================================================================
SUBROUTINE compute_lake_frost4PCM(min_yrPCM_lakeice)
!-----------------------------------------------------------------------
! NAME
!     compute_lake_frost4PCM
!
! DESCRIPTION
!     Compute the seasonal lake-ice deviation to give back to the PCM.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use display, only: print_msg, LVL_NFO
use utility, only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in) :: min_yrPCM_lakeice

! CODE
! ----
where (lakeice_PCM(:,:) > 0._dp) lake_frost4PCM(:,:) = lakeice_PCM(:,:) - min_yrPCM_lakeice(:,:)
call print_msg('Yearly minimum of lake ice [kg/m2] (min|max): '//real2str(minval(min_yrPCM_lakeice))//' | '//real2str(maxval(min_yrPCM_lakeice)),LVL_NFO)
call print_msg('Raw input of lake ice      [kg/m2] (min|max): '//real2str(minval(lakeice_PCM))//' | '//real2str(maxval(lakeice_PCM)),LVL_NFO)
call print_msg('Residual lake frost        [kg/m2] (min|max): '//real2str(minval(lake_frost4PCM))//' | '//real2str(maxval(lake_frost4PCM)),LVL_NFO)

END SUBROUTINE compute_lake_frost4PCM
!=======================================================================

!=======================================================================
SUBROUTINE build4PCM_lake_ice(lakeice,lakeice4PCM)
!-----------------------------------------------------------------------
! NAME
!     build4PCM_lake_ice
!
! DESCRIPTION
!     Reconstructs lake ice plus seasonal deviation for the PCM.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use display, only: print_msg, LVL_NFO

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: lakeice
real(dp), dimension(:,:), intent(out) :: lakeice4PCM

! CODE
! ----
call print_msg('> Building lake ice for the PCM',LVL_NFO)
lakeice4PCM(:,:) = lakeice(:,:) + lake_frost4PCM(:,:)
where (lakeice4PCM(:,:) < 0._dp) lakeice4PCM(:,:) = 0._dp

END SUBROUTINE build4PCM_lake_ice
!=======================================================================

!=======================================================================
SUBROUTINE evolve_lake_ice(lakeice,d_lakeice,h2o_ice,tsurf_avg,lake_weighted_heatcap_surf)
!-----------------------------------------------------------------------
! NAME
!     evolve_lake_ice
!
! DESCRIPTION
!     Evolve lake ice over one PEM time step with the lake-ice growth
!     formula used for perennial liquid-water cells.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,  only: ngrid, nslope
use evolution, only: dt, get_dt_yr
use hydrology, only: lake_frac, compute_lake_water_mass_per_area
use stoppage,  only: stop_clean
use display,   only: print_msg, LVL_WRN
use utility,   only: int2str, real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)    :: h2o_ice, tsurf_avg, lake_weighted_heatcap_surf
real(dp), dimension(:,:), intent(inout) :: lakeice
real(dp), dimension(:,:), intent(out)   :: d_lakeice

! LOCAL VARIABLES
! ---------------
integer(di)                :: ig, islope
real(dp)                   :: dt_yr ! Time step in Planetary years
real(dp)                   :: h0, h_snow, F, R, radicand, h_new, lakeice_old
real(dp), dimension(ngrid) :: lake_water_mass_per_area

! CODE
! ----
call compute_lake_water_mass_per_area(lake_water_mass_per_area)
d_lakeice(:,:) = 0._dp
dt_yr = get_dt_yr()

do islope = 1,nslope
    do ig = 1,ngrid
        lakeice_old = max(0._dp,lakeice(ig,islope))
        if (lake_frac(ig) <= 0._dp) then
            lakeice(ig,islope) = 0._dp
            d_lakeice(ig,islope) = (lakeice(ig,islope) - lakeice_old)/dt_yr
            cycle
        end if

        if (lake_weighted_heatcap_surf(ig,islope) <= 0._dp) then
            call stop_clean(__FILE__,__LINE__,'''lake_weighted_heatcap_surf'' must be positive in lake-ice evolution!',1)
        end if

        h0 = lakeice_old/rhowaterice
        h_snow = max(0._dp,h2o_ice(ig,islope))/rho_snow
        F = (T_h2o_ice_liq - tsurf_avg(ig,islope))/(rhowaterice*RLFTT)
        R = lake_frac(ig)*dt/lake_weighted_heatcap_surf(ig,islope) + h_snow/k_snow
        radicand = (h0 + k_ice*R)**2 + 2._dp*k_ice*F*dt
        if (radicand <= 0._dp) then
            h_new = 0._dp
        else
            h_new = max(0._dp,sqrt(radicand) - k_ice*R)
        end if
        lakeice(ig,islope) = h_new*rhowaterice

        if (lakeice(ig,islope) > lake_water_mass_per_area(ig)) then
            call print_msg('Lake ice exceeds available lake water at grid '//int2str(ig)//', slope '//int2str(islope)// &
                           ': lakeice='//real2str(lakeice(ig,islope))//' kg/m2, available='//                           &
                           real2str(lake_water_mass_per_area(ig))//' kg/m2.',LVL_WRN)
        end if
        d_lakeice(ig,islope) = (lakeice(ig,islope) - lakeice_old)/dt_yr
    end do
end do

END SUBROUTINE evolve_lake_ice
!=======================================================================

END MODULE lake_ice