MODULE surf_rad_prop
!-----------------------------------------------------------------------
! NAME
!     surf_rad_prop
!
! DESCRIPTION
!     Surface radiative properties (albedo and emissivity) for the PCM.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! VISIBILITY
! ----------
private :: read_albedo_emis

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

!=======================================================================
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
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,           only: ngrid, nslope, latitudes
use surf_rad_prop_core, only: albedodat_PCM
use frost,              only: h2o_frost4PCM, co2_frost4PCM
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, 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 emissivity 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

        ! H2O ice
        if (h2o_ice4PCM(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_ice4PCM(i,islope) > 0._dp) then
            albedo4PCM(i,islope) = albedo_co2ice(icap)
            emissivity4PCM(i,islope) = emissivity_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._dp) 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 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)                         :: cst_h2oice_albedo
integer(di)                         :: i, nindex
real(dp), dimension(:), allocatable :: controle

! CODE
! ----
! 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'' 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'' 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'' 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'' is out of bounds [0,1]!',1)

! Read albedos of CO2 frost 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)
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 surf_rad_prop