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_surfice4PCM,co2_surfice4PCM,peren_h2o_ice_frac4PCM,peren_h2o_liq_frac4PCM, &
                                   is_ocean4PCM,lakeice4PCM,seaice4PCM,seaice_frac4PCM,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 surf_rad_prop_core, only: albedodat_PCM
use display,            only: print_msg, LVL_NFO

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:), intent(in)  :: h2o_surfice4PCM, co2_surfice4PCM
real(dp),    dimension(:),   intent(in)  :: peren_h2o_liq_frac4PCM, lakeice4PCM, seaice4PCM, seaice_frac4PCM
logical(k4), dimension(:),   intent(in)  :: is_ocean4PCM
real(dp),    dimension(:),   intent(in)  :: peren_h2o_ice_frac4PCM
real(dp),    dimension(:,:), intent(out) :: albedo4PCM, emissivity4PCM

! PARAMETERS
! ----------
real(dp), parameter :: snowlayer = 33._dp    ! PCM critical snow amount for albedo [kg/m2] (3.3 cm of snow)
real(dp), parameter :: snow_min = 15._dp     ! PCM critical snow amount over sea/lake ice [kg/m2] (5 cm of snow)
real(dp), parameter :: alb_ice_min = 0.08_dp ! PCM minimum albedo for thin-ice asymptotic value [kg/m2]
real(dp), parameter :: h_alb_ice = 275.1_dp  ! PCM scale for ice-thickness albedo [kg/m2]

! LOCAL VARIABLES
! ---------------
integer(di)            :: i, islope, icap
real(dp)               :: emissivity_baresoil, albedo_h2oice, albedo_co2ice, albedo_ocean
real(dp)               :: frac_snow, albedo_ice, albedo_land, albedo_lake, frac_shelf, frac_sic
real(dp)               :: frac_land, frac_lake, frac_glacier
real(dp), dimension(2) :: emissivity_co2ice

! CODE
! ----
! Fetch parameters from "callphys.def" and "startfi.nc"
call read_albedo_emis(albedo_h2oice,albedo_co2ice,albedo_ocean,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
        ! PCM sequence step 1 (condense_co2): initialize albedo/emissivity from CO2 cover or bare ground
        if (co2_surfice4PCM(i,islope) > 1._dp) then ! CO2 ice is taken into account if its thickness is > 1 mm
            albedo4PCM(i,islope) = albedo_co2ice
            emissivity4PCM(i,islope) = emissivity_co2ice(icap)
        else ! Bare ground
            albedo4PCM(i,islope) = albedodat_PCM(i)
            emissivity4PCM(i,islope) = emissivity_baresoil
        end if

        ! PCM sequence step 2 (hydrol): rebuild albedo from hydrology state (this can overwrite step 1 albedo)
        ! H2O ice on the land
        albedo_land = albedodat_PCM(i)
        if (h2o_surfice4PCM(i,islope) > 0._dp) then
            if (h2o_surfice4PCM(i,islope) >= snowlayer) then
                albedo_land = albedo_h2oice
            else
                ! In the Generic PCM, if H2O ice is present but less than 3.3 cm,
                ! the albedo is a weighted average between the bare ground and H2O ice albedo.
                albedo_land = albedo_land + (albedo_h2oice - albedo_land)*h2o_surfice4PCM(i,islope)/snowlayer
            end if
        end if

        ! Ocean grid cell (dominant over land)
        if (is_ocean4PCM(i)) then
            frac_shelf = min(peren_h2o_ice_frac4PCM(i),1._dp) ! The ice shelf covers a fraction of the cell
            albedo4PCM(i,islope) = frac_shelf*albedo_h2oice + (1._dp - frac_shelf)*albedo_ocean
            ! Sea ice (dominant over ocean), on the open water left by the ice shelf
            if (seaice_frac4PCM(i) > 0._dp .and. seaice4PCM(i) > 0._dp) then
                frac_sic = min(seaice_frac4PCM(i),1._dp - frac_shelf)
                frac_snow = max(0._dp,min(1._dp,h2o_surfice4PCM(i,islope)/snow_min))
                albedo_ice = albedo_h2oice - (albedo_h2oice - alb_ice_min)*exp(-seaice4PCM(i)/h_alb_ice)
                albedo4PCM(i,islope) = frac_sic*(albedo_h2oice*frac_snow + albedo_ice*(1._dp - frac_snow)) + &
                                       frac_shelf*albedo_h2oice + (1._dp - frac_sic - frac_shelf)*albedo_ocean
            end if
        else ! Land grid cell with lakes
            ! Match the PCM convention: lakes and glaciers overlap only when their summed areas exceed the cell.
            frac_land = 1._dp - min(peren_h2o_liq_frac4PCM(i) + peren_h2o_ice_frac4PCM(i),1._dp)
            frac_lake = max(0._dp,1._dp - frac_land - peren_h2o_ice_frac4PCM(i))
            frac_glacier = peren_h2o_ice_frac4PCM(i)
            frac_snow = max(0._dp,min(1._dp,h2o_surfice4PCM(i,islope)/snow_min))
            albedo_lake = albedo_ocean + frac_snow*(albedo_h2oice - albedo_ocean)
            ! Lake ice (dominant over lake)
            if (frac_lake > 0._dp .and. lakeice4PCM(i) > 0._dp) then
                albedo_ice = albedo_h2oice - (albedo_h2oice - alb_ice_min)*exp(-lakeice4PCM(i)/h_alb_ice)
                albedo_lake = albedo_h2oice*frac_snow + albedo_ice*(1._dp - frac_snow)
            end if
            albedo4PCM(i,islope) = frac_land*albedo_land + frac_lake*albedo_lake + frac_glacier*albedo_h2oice
        end if

        ! PCM sequence step 3 (hydrol end): re-add CO2 albedo override after hydrology update
        if (co2_surfice4PCM(i,islope) > 1._dp) albedo4PCM(i,islope) = albedo_co2ice

        ! Guard against unphysical values
        albedo4PCM(i,islope) = max(0._dp,min(1._dp,albedo4PCM(i,islope)))
        emissivity4PCM(i,islope) = max(0._dp,min(1._dp,emissivity4PCM(i,islope)))
    end do
end do

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

!=======================================================================
SUBROUTINE read_albedo_emis(albedo_h2oice,albedo_co2ice,albedo_ocean,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
!     Be careful: in the PEM, albedo is read from "callphys.def" as a
!     scalar value, but in the Generic PCM, it is calculated spectrally.
!-----------------------------------------------------------------------

! 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, albedo_ocean, 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"
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 ocean albedo from "callphys.def"
albedo_ocean = 0.07_dp ! Default
call getin('alb_ocean',albedo_ocean)
if (albedo_ocean < 0._dp .or. albedo_ocean > 1._dp) call stop_clean(__FILE__,__LINE__,'''alb_ocean'' 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 surf_rad_prop