MODULE surface
!-----------------------------------------------------------------------
! NAME
!     surface
!
! DESCRIPTION
!     Contains global parameters used for the surface.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
real(dp), dimension(:),   allocatable, protected, private :: albedodat_PCM  ! Albedo of bare ground
real(dp), dimension(:,:), allocatable, protected, private :: albedo_PCM     ! Surface albedo_PCM
real(dp), dimension(:,:), allocatable, protected          :: emissivity_PCM ! Thermal IR surface emissivity_PCM

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

!=======================================================================
SUBROUTINE ini_surface()
!-----------------------------------------------------------------------
! NAME
!     ini_surface
!
! DESCRIPTION
!     Initialize the parameters of module 'surface'.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
allocate(albedo_PCM(ngrid,nslope))
allocate(albedodat_PCM(ngrid))
allocate(emissivity_PCM(ngrid,nslope))
albedo_PCM(:,:) = 0._dp
albedodat_PCM(:) = 0._dp
emissivity_PCM(:,:) = 1._dp

END SUBROUTINE ini_surface
!=======================================================================

!=======================================================================
SUBROUTINE end_surface()
!-----------------------------------------------------------------------
! NAME
!     end_surface
!
! DESCRIPTION
!     Deallocate surface arrays.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
if (allocated(albedo_PCM)) deallocate(albedo_PCM)
if (allocated(albedodat_PCM)) deallocate(albedodat_PCM)
if (allocated(emissivity_PCM)) deallocate(emissivity_PCM)

END SUBROUTINE end_surface
!=======================================================================

!=======================================================================
SUBROUTINE set_albedodat_PCM(albedodat_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_albedodat_PCM
!
! DESCRIPTION
!     Setter for 'albedodat_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
albedodat_PCM(:) = albedodat_PCM_in(:)

END SUBROUTINE set_albedodat_PCM
!=======================================================================

!=======================================================================
SUBROUTINE set_albedo_PCM(albedo_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_albedo_PCM
!
! DESCRIPTION
!     Setter for 'albedo_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
albedo_PCM(:,:) = albedo_PCM_in(:,:)

END SUBROUTINE set_albedo_PCM
!=======================================================================

!=======================================================================
SUBROUTINE set_emissivity_PCM(emissivity_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_emissivity_PCM
!
! DESCRIPTION
!     Setter for 'emissivity_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
emissivity_PCM(:,:) = emissivity_PCM_in(:,:)

END SUBROUTINE set_emissivity_PCM
!=======================================================================

!=======================================================================
SUBROUTINE build4PCM_surf_rad_prop(h2o_ice4PCM,co2_ice4PCM,albedo4PCM,emissivity4PCM)
!-----------------------------------------------------------------------
! NAME
!     build4PCM_surf_rad_prop
!
! DESCRIPTION
!     Reconstructs albedo and emissivity for the PCM.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!     C. Metz, 01/2026 & 04/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: h2o_ice4PCM, co2_ice4PCM
real(dp), dimension(:,:), intent(out) :: albedo4PCM, emissivity4PCM

! LOCAL VARIABLES
! ---------------
integer(di)            :: i, islope, icap
real(dp)               :: albedo_h2oice, emissivity_baresoil, albedo_co2ice
real(dp), dimension(2) :: emissivity_co2ice

! CODE
! ----
! Fetch parameters from "callphys.def" and "startfi.nc"
call read_albedo_emis(albedo_h2oice,albedo_co2ice,emissivity_co2ice,emissivity_baresoil)

! Reconstruction Loop
call print_msg('> Building albedo and emmissivity for the PCM',LVL_NFO)
do i = 1,ngrid
    ! Determine hemisphere: 1 = Northern, 2 = Southern
    if (latitudes(i) < 0._dp) then
        icap = 2
    else
        icap = 1
    end if

    do islope = 1,nslope
        ! Bare ground (default)
        albedo4PCM(i,islope) = albedodat_PCM(i)
        emissivity4PCM(i,islope) = emissivity_baresoil

        ! CO2 frost
        ! cmetz co2 cond called before hydrol in physiq_mod, so water albedo layers over co2
        ! It's gt 1 or some small number in the gPCM
        ! cmetz I think it's gt 1 or some small number in the gPCM
        if (co2_ice4PCM(i,islope) > 1._dp) then
            albedo4PCM(i,islope) = albedo_co2ice
            emissivity4PCM(i,islope) = emissivity_co2ice(icap)
        end if

        ! H2O frost (dominant over CO2 ice)
        ! cmetz in the PCM we use a bucket model kind of thing
        ! 33 comes from parameter "snowlayer" in PCM, which is the size of the albedo bucket, and is hardcoded in hydrol.F90
        if (h2o_ice4PCM(i,islope) > 0._dp) then
            albedo4PCM(i,islope) = albedo4PCM(i,islope) + (albedo_h2oice - albedo4PCM(i,islope))*h2o_ice4PCM(i,islope)/33.0
            emissivity4PCM(i,islope) = 1._dp
        end if
    end do
end do

END SUBROUTINE build4PCM_surf_rad_prop
!=======================================================================

!=======================================================================
SUBROUTINE read_albedo_emis(albedo_h2oice,albedo_co2ice,emissivity_co2ice,emissivity_baresoil)
!-----------------------------------------------------------------------
! NAME
!     read_albedo_emis
!
! DESCRIPTION
!     Reads albedo/emissivity parameters from "callphys.def" and
!     "startfi.nc" ('controle').
!
! AUTHORS & DATE
!     C. Metz, 01/2026 & 04/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
#ifdef CPP_IOIPSL
use IOIPSL,          only: getin
#else
use ioipsl_getincom, only: getin
#endif
use io_netcdf,       only: open_nc, close_nc, startfi_name, get_dim_nc, get_var_nc
use stoppage,        only: stop_clean

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

! ARGUMENTS
! ---------
real(dp),               intent(out) :: albedo_h2oice, albedo_co2ice, emissivity_baresoil
real(dp), dimension(2), intent(out) :: emissivity_co2ice

! LOCAL VARIABLES
! ---------------
integer(di)                         :: i, nindex
real(dp), dimension(:), allocatable :: controle

! CODE
! ----
! Read albedos of H2O ice from "callphys.def"
! cmetz right now we get the scalar values which are not used in the generic PCM, they are in fact calculated spectrally...
albedo_h2oice = 0.35_dp ! Default
call getin('albedosnow',albedo_h2oice)
if (albedo_h2oice < 0._dp .or. albedo_h2oice > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedosnow'' is out of bounds [0,1]!',1)

! Read albedos of CO2 ice from "callphys.def"
albedo_co2ice = 0.6_dp ! Default
call getin('albedoco2ice',albedo_co2ice)
if (albedo_co2ice < 0._dp .or. albedo_co2ice > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedoco2ice'' is out of bounds [0,1]!',1)

! Read emissivity of CO2 ice from "startfi.nc" ('controle' variable)
call open_nc(startfi_name,'read')
call get_dim_nc('index',nindex)
allocate(controle(nindex))
call get_var_nc('controle',controle)
call close_nc(startfi_name)
emissivity_co2ice(1) = controle(24)
emissivity_co2ice(2) = controle(25)
emissivity_baresoil  = controle(26)
deallocate(controle)
do i = 1,2
    if (emissivity_co2ice(i) < 0._dp .or. emissivity_co2ice(i) > 1._dp) call stop_clean(__FILE__,__LINE__,'''emissivity_co2ice'' from ''controle(24:25)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)
end do
if (emissivity_baresoil < 0._dp .or. emissivity_baresoil > 1._dp) call stop_clean(__FILE__,__LINE__,'''emissivity_baresoil'' from ''controle(26)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)

END SUBROUTINE read_albedo_emis
!=======================================================================

END MODULE surface