MODULE sea_ice
!-----------------------------------------------------------------------
! NAME
!     sea_ice
!
! DESCRIPTION
!     Module to manage perennial sea ice in the PEM.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! VARIABLES
! ---------
! Sea ice in the PCM at the beginning and seasonal deviation to restore for PCM handoff
real(dp), dimension(:,:), allocatable, protected :: seaice_PCM
real(dp), dimension(:,:), allocatable, protected :: seaice_frac_PCM
real(dp), dimension(:,:), allocatable            :: seaice_dev4PCM
real(dp), dimension(:,:), allocatable            :: seaice_frac_dev4PCM

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

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

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

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

! CODE
! ----
allocate(seaice_PCM(ngrid,nslope))
allocate(seaice_frac_PCM(ngrid,nslope))
allocate(seaice_dev4PCM(ngrid,nslope))
allocate(seaice_frac_dev4PCM(ngrid,nslope))
seaice_PCM(:,:) = 0._dp
seaice_frac_PCM(:,:) = 0._dp
seaice_dev4PCM(:,:) = 0._dp
seaice_frac_dev4PCM(:,:) = 0._dp

END SUBROUTINE ini_sea_ice
!=======================================================================

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

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

! CODE
! ----
if (allocated(seaice_PCM)) deallocate(seaice_PCM)
if (allocated(seaice_frac_PCM)) deallocate(seaice_frac_PCM)
if (allocated(seaice_dev4PCM)) deallocate(seaice_dev4PCM)
if (allocated(seaice_frac_dev4PCM)) deallocate(seaice_frac_dev4PCM)

END SUBROUTINE end_sea_ice
!=======================================================================

!=======================================================================
SUBROUTINE set_seaice_PCM(seaice_PCM_in, seaice_frac_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_seaice_PCM
!
! DESCRIPTION
!     Setter for PCM sea-ice state read from "startfi.nc".
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

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

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

! CODE
! ----
do islope = 1,nslope
    seaice_PCM(:,islope) = max(0._dp,seaice_PCM_in(:))
    seaice_frac_PCM(:,islope) = min(1._dp,max(0._dp,seaice_frac_PCM_in(:)))
end do
call clamp_seaice_pair(seaice_PCM,seaice_frac_PCM)

END SUBROUTINE set_seaice_PCM
!=======================================================================

!=======================================================================
SUBROUTINE compute_sea_ice_dev4PCM(min_yrPCM_seaice, min_yrPCM_seaice_frac)
!-----------------------------------------------------------------------
! NAME
!     compute_sea_ice_dev4PCM
!
! DESCRIPTION
!     Compute seasonal sea-ice deviations 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_seaice, min_yrPCM_seaice_frac

! CODE
! ----
seaice_dev4PCM(:,:) = 0._dp
where (seaice_PCM(:,:) > 0._dp) seaice_dev4PCM(:,:) = seaice_PCM(:,:) - min_yrPCM_seaice(:,:)

seaice_frac_dev4PCM(:,:) = 0._dp
where (seaice_frac_PCM(:,:) > 0._dp) seaice_frac_dev4PCM(:,:) = seaice_frac_PCM(:,:) - min_yrPCM_seaice_frac(:,:)

call print_msg('Yearly minimum of sea ice [kg/m2] (min|max): '//real2str(minval(min_yrPCM_seaice))//' | '// &
               real2str(maxval(min_yrPCM_seaice)),LVL_NFO)
call print_msg('Raw input of sea ice      [kg/m2] (min|max): '//real2str(minval(seaice_PCM))//' | '// &
               real2str(maxval(seaice_PCM)),LVL_NFO)
call print_msg('Yearly minimum of sea ice fraction (min|max): '//real2str(minval(min_yrPCM_seaice_frac))//' | '// &
               real2str(maxval(min_yrPCM_seaice_frac)),LVL_NFO)
call print_msg('Raw input of sea ice fraction      (min|max): '//real2str(minval(seaice_frac_PCM))//' | '// &
               real2str(maxval(seaice_frac_PCM)),LVL_NFO)

END SUBROUTINE compute_sea_ice_dev4PCM
!=======================================================================

!=======================================================================
SUBROUTINE compute_sea_ice_tendencies(min_yrPCM_seaice, min_yrPCM_seaice_frac, seaice, seaice_frac, d_seaice, d_seaice_frac)
!-----------------------------------------------------------------------
! NAME
!     compute_sea_ice_tendencies
!
! DESCRIPTION
!     Compute sea-ice and sea-ice-fraction tendencies from yearly minima.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:,:), intent(in)  :: min_yrPCM_seaice, min_yrPCM_seaice_frac
real(dp), dimension(:,:),   intent(in)  :: seaice, seaice_frac
real(dp), dimension(:,:),   intent(out) :: d_seaice, d_seaice_frac

! CODE
! ----
d_seaice(:,:) = min_yrPCM_seaice(:,:,2) - min_yrPCM_seaice(:,:,1)
where (abs(d_seaice) < eps) d_seaice = 0._dp
where (d_seaice(:,:) < 0._dp .and. abs(seaice(:,:)) < eps) d_seaice(:,:) = 0._dp

d_seaice_frac(:,:) = min_yrPCM_seaice_frac(:,:,2) - min_yrPCM_seaice_frac(:,:,1)
where (abs(d_seaice_frac) < eps) d_seaice_frac = 0._dp
where (d_seaice_frac(:,:) < 0._dp .and. abs(seaice_frac(:,:)) < eps) d_seaice_frac(:,:) = 0._dp

END SUBROUTINE compute_sea_ice_tendencies
!=======================================================================

!=======================================================================
SUBROUTINE evolve_sea_ice(seaice, seaice_frac, d_seaice, d_seaice_frac)
!-----------------------------------------------------------------------
! NAME
!     evolve_sea_ice
!
! DESCRIPTION
!     Evolve perennial sea ice over one PEM time step.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use evolution, only: get_dt_yr

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(inout) :: seaice, seaice_frac
real(dp), dimension(:,:), intent(in)    :: d_seaice, d_seaice_frac

! LOCAL VARIABLES
! ---------------
real(dp) :: dt_yr ! Time step in Planetary years

! CODE
! ----
dt_yr = get_dt_yr()
seaice(:,:) = seaice(:,:) + d_seaice(:,:)*dt_yr
seaice_frac(:,:) = seaice_frac(:,:) + d_seaice_frac(:,:)*dt_yr
call clamp_seaice_pair(seaice,seaice_frac)

END SUBROUTINE evolve_sea_ice
!=======================================================================

!=======================================================================
SUBROUTINE build4PCM_sea_ice(seaice, seaice_frac, seaice4PCM, seaice_frac4PCM)
!-----------------------------------------------------------------------
! NAME
!     build4PCM_sea_ice
!
! DESCRIPTION
!     Reconstruct sea 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)  :: seaice, seaice_frac
real(dp), dimension(:),   intent(out) :: seaice4PCM, seaice_frac4PCM

! LOCAL VARIABLES
! ---------------
real(dp), dimension(size(seaice,1),size(seaice,2)) :: seaice_tmp, seaice_frac_tmp

! CODE
! ----
call print_msg('> Building sea ice for the PCM',LVL_NFO)
seaice_tmp(:,:) = seaice(:,:) + seaice_dev4PCM(:,:)
seaice_frac_tmp(:,:) = seaice_frac(:,:) + seaice_frac_dev4PCM(:,:)
call clamp_seaice_pair(seaice_tmp,seaice_frac_tmp)

seaice4PCM(:) = seaice_tmp(:,1)
seaice_frac4PCM(:) = seaice_frac_tmp(:,1)

END SUBROUTINE build4PCM_sea_ice
!=======================================================================

!=======================================================================
SUBROUTINE clamp_seaice_pair(seaice, seaice_frac)
!-----------------------------------------------------------------------
! NAME
!     clamp_seaice_pair
!
! DESCRIPTION
!     Enforce physical bounds and keep sea ice and fraction coupled.
!
! AUTHORS & DATE
!     C. Metz, 06/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(inout) :: seaice, seaice_frac

! CODE
! ----
where (seaice(:,:) < 0._dp) seaice(:,:) = 0._dp
where (seaice_frac(:,:) < 0._dp) seaice_frac(:,:) = 0._dp
where (seaice_frac(:,:) > 1._dp) seaice_frac(:,:) = 1._dp
where (seaice_frac(:,:) <= eps) seaice(:,:) = 0._dp
where (seaice(:,:) <= eps) seaice_frac(:,:) = 0._dp

END SUBROUTINE clamp_seaice_pair
!=======================================================================

END MODULE sea_ice