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 :: albedodat_PCM  ! Albedo of bare ground
real(dp), dimension(:,:), allocatable, protected :: 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
! ----
if (.not. allocated(albedo_PCM)) allocate(albedo_PCM(ngrid,nslope))
if (.not. allocated(albedodat_PCM)) allocate(albedodat_PCM(ngrid))
if (.not. allocated(emissivity_PCM)) allocate(emissivity_PCM(ngrid,nslope))

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_ice,co2_ice,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
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, latitudes
use frost,    only: h2o_frost4PCM, co2_frost4PCM
use display,  only: print_msg

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: h2o_ice, co2_ice
real(dp), dimension(:,:), intent(out) :: albedo4PCM, emissivity4PCM

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

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

! Reconstruction Loop
call print_msg('> Building albedo and emmissivity for the PCM')
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

        ! H2O ice
        if (h2o_ice(i,islope) > 0._dp) then
            albedo4PCM(i,islope) = albedo_h2oice
            emissivity4PCM(i,islope) = 1._dp
        end if

        ! CO2 ice (dominant over H2O ice)
        if (co2_ice(i,islope) > 0._dp) then
            albedo4PCM(i,islope) = albedo_co2ice(icap)
            emissivity4PCM(i,islope) = albedo_co2ice(icap)
        end if

        ! H2O frost (subject to threshold)
        if (h2o_frost4PCM(i,islope) > h2ofrost_albedo_threshold) then
            albedo4PCM(i,islope) = albedo_h2ofrost
            emissivity4PCM(i,islope) = 1._dp
        end if

        ! CO2 frost (final priority)
        if (co2_frost4PCM(i,islope) > 0.) then
            albedo4PCM(i,islope) = albedo_co2frost(icap)
            emissivity4PCM(i,islope) = emissivity_co2ice(icap)
        end if
    end do
end do

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

!=======================================================================
SUBROUTINE read_albedo_emis(albedo_h2oice,albedo_h2ofrost,h2ofrost_albedo_threshold, &
                            albedo_co2ice,albedo_co2frost,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
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
#ifdef CPP_IOIPSL
use IOIPSL,          only: getin
#else
use ioipsl_getincom, only: getin
#endif
use config,          only: callphys_name
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_h2ofrost, h2ofrost_albedo_threshold, emissivity_baresoil
real(dp), dimension(2), intent(out) :: albedo_co2frost, emissivity_co2ice, albedo_co2ice

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

! CODE
! ----
inquire(file = callphys_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find "'//callphys_name//'"!', 1)

! Read albedos of H2O frost and perennial ice from "callphys.def"
albedo_h2oice = 0.35_dp ! Default
call getin('albedo_h2o_cap',albedo_h2oice)
if (albedo_h2oice < 0._dp .or. albedo_h2oice > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_h2o_cap'' in "'//callphys_name//'" is out of bounds [0,1]!',1)

albedo_h2ofrost = 0.35_dp ! Default
call getin('albedo_h2ofrost',albedo_h2ofrost)
if (albedo_h2ofrost < 0._dp .or. albedo_h2ofrost > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_h2ofrost'' in "'//callphys_name//'" is out of bounds [0,1]!',1)

h2ofrost_albedo_threshold = 0.005_dp ! Default [kg/m2]
call getin('frost_albedo_threshold',h2ofrost_albedo_threshold)

cst_h2oice_albedo = .false. ! Default
call getin('cst_cap_albedo',cst_h2oice_albedo)
if (cst_h2oice_albedo) albedo_h2ofrost = albedo_h2oice ! If true, then we don't account for water frost albedo effect

! Read albedos of CO2 perennial ice from "callphys.def"
albedo_co2ice(1) = 0.6_dp ! Default, north hemisphere
albedo_co2ice(2) = 0.85_dp ! Default, south hemisphere
call getin('albedo_co2ice_north',albedo_co2ice(1))
call getin('albedo_co2ice_south',albedo_co2ice(2))
if (albedo_co2ice(1) < 0._dp .or. albedo_co2ice(1) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2ice_north'' in "'//callphys_name//'" is out of bounds [0,1]!',1)
if (albedo_co2ice(2) < 0._dp .or. albedo_co2ice(2) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2ice_south'' in "'//callphys_name//'" is out of bounds [0,1]!',1)

! Read albedos of CO2 frost from "callphys.def"
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)
albedo_co2frost(1)   = controle(22)
albedo_co2frost(2)   = controle(23)
emissivity_co2ice(1) = controle(24)
emissivity_co2ice(2) = controle(25)
emissivity_baresoil  = controle(26)
deallocate(controle)
do i = 1,2
    if (albedo_co2frost(i) < 0._dp .or. albedo_co2frost(i) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2frost'' from ''controle(22:23)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)
    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