MODULE surf_temp
!-----------------------------------------------------------------------
! NAME
!     surf_temp
!
! DESCRIPTION
!     Surface temperature management.
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMATERS
! ----------
real(dp), dimension(:,:), allocatable, protected :: tsurf_PCM ! Surface temperature in the PCM at the beginning [K]

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

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

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

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

! CODE
! ----
if (.not. allocated(tsurf_PCM)) allocate(tsurf_PCM(ngrid,nslope))

END SUBROUTINE ini_surf_temp
!=======================================================================

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

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

! CODE
! ----
if (allocated(tsurf_PCM)) deallocate(tsurf_PCM)

END SUBROUTINE end_surf_temp
!=======================================================================

!=======================================================================
SUBROUTINE set_tsurf_PCM(tsurf_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_tsurf_PCM
!
! DESCRIPTION
!     Setter for 'tsurf_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! CODE
! ----
tsurf_PCM(:,:) = tsurf_PCM_in(:,:)

END SUBROUTINE set_tsurf_PCM
!=======================================================================

!=======================================================================
SUBROUTINE adapt_tsurf2disappearedice(surfice,is_ice_ini,ice_disappeared,tsurf_avg)
!-----------------------------------------------------------------------
! NAME
!     adapt_tsurf2disappearedice
!
! DESCRIPTION
!     Update surface temperature where CO2 ice disappeared using nearby
!     bare soil temperature.
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nlon, nlat, latitudes, vect2lonlat, lonlat2vect
use display,  only: print_msg
use utility,  only: int2str

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

! ARGUMENTS
! ---------
real(dp),    dimension(:,:), intent(in)    :: surfice
logical(k4), dimension(:,:), intent(in)    :: is_ice_ini
real(dp),    dimension(:,:), intent(inout) :: tsurf_avg
logical(k4), dimension(:,:), intent(inout) :: ice_disappeared

! LOCAL VARIABLES
! ---------------
integer(di)                              :: islope, i, j, k, radius, rmax, d_i, d_j, ii, jj
logical(k4)                              :: found
real(dp),    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, mask_co2ice_disappeared
real(dp),    dimension(nlon,nlat)        :: latitude_ll
real(dp),    dimension(ngrid)            :: tmp
integer(di), dimension(nslope - 1)       :: priority

! CODE
! ----
! Check to escape the subroutine (not relevant in 1D)
if (ngrid == 1) return

call print_msg("> Adapting surface temperature where ice disappeared")
! Convert from reduced grid to lon-lat grid
call vect2lonlat(nlon,nlat,ngrid,latitudes,latitude_ll)
do islope = 1,nslope
    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,surfice(:,islope),co2ice_ll(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,is_ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,ice_disappeared(:,islope)),mask_co2ice_disappeared(:,:,islope))
end do

! For each point where ice disappeared
rmax = max(nlon,nlat)
do j = 1,nlat
    do i = 1,nlon
        do islope = 1,nslope
            if (mask_co2ice_ini(i,j,islope) > 0.5_dp .and. co2ice_ll(i,j,islope) < minieps .and. mask_co2ice_disappeared(i,j,islope) < 0.5_dp) then
                found = .false.
                mask_co2ice_disappeared(i,j,islope) = 1._dp
                call get_slope_priority(latitude_ll(i,j),islope,priority)
                do k = 1,nslope - 1
                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5_dp) then
                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
                        found = .true.
                        exit
                    end if
                end do

                radius = 1
                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
                    do d_j = -radius,radius
                        do d_i = -radius,radius
                            if (abs(d_i) + abs(d_j) == radius) then
                                ii = i + d_i
                                jj = j + d_j
                                ! Longitudinal periodicity
                                if (ii < 1) then
                                    ii = ii + nlon
                                else if (ii > nlon) then
                                    ii = ii - nlon
                                end if
                                ! Latitude boundaries
                                if (jj >= 1 .and. jj <= nlat) then
                                    call get_slope_priority(latitude_ll(ii,jj),islope,priority)
                                    do k = 1,nslope - 1
                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5_dp) then
                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
                                            found = .true.
                                            exit
                                        end if
                                    end do
                                end if
                            end if
                            if (found) exit
                        end do
                        if (found) exit
                    end do
                    radius = radius + 1
                end do
                if (.not. found) call print_msg("Warning: no bare soil found for ice disappeared at point ("//int2str(i)//','//int2str(j)//','//int2str(islope)//'!')
            end if
        end do
    end do
end do

! Convert back from lon-lat grid to reduced grid
do islope = 1,nslope
    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
    call lonlat2vect(nlon,nlat,ngrid,mask_co2ice_disappeared(:,:,islope),tmp)
    where (tmp > 0.5_dp) ice_disappeared(:,islope) = .true.
end do

END SUBROUTINE adapt_tsurf2disappearedice
!=======================================================================

!=======================================================================
SUBROUTINE get_slope_priority(lat,islope,priority)
!-----------------------------------------------------------------------
! NAME
!     get_slope_priority
!
! DESCRIPTION
!     Determine slope priority based on latitudes (equator-ward favored).
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!     Equator-ward slopes are most likely to hold no ice.
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
real(dp),                  intent(in)  :: lat      ! Latitude [degrees]
integer(di),               intent(in)  :: islope
integer(di), dimension(:), intent(out) :: priority ! Priority ordering of slopes

! LOCAL VARIABLES
! ---------------
integer(di) :: i, k

! CODE
! ----
k = 1

! Northern hemisphere
if (lat > 0.) then
    ! Equator-ward slopes
    do i = islope - 1,1,-1
        priority(k) = i
        k = k + 1
    end do
    ! Pole-ward slopes
    do i = islope + 1,nslope
        priority(k) = i
        k = k + 1
    end do
else ! Southern hemisphere
    ! Equator-ward slopes
    do i = islope + 1,nslope
        priority(k) = i
        k = k + 1
    end do
    ! Pole-ward slopes
    do i = islope - 1,1,-1
        priority(k) = i
        k = k + 1
    end do
end if

END SUBROUTINE get_slope_priority
!=======================================================================

!=======================================================================
SUBROUTINE build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)
!-----------------------------------------------------------------------
! NAME
!     build4PCM_tsurf
!
! DESCRIPTION
!     Reconstructs surface temperature for the PCM.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use display, only: print_msg

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: tsurf_avg, tsurf_dev
real(dp), dimension(:,:), intent(out) :: tsurf4PCM

! CODE
! ----
call print_msg('> Building surface temperature for the PCM')
tsurf4PCM(:,:) = tsurf_avg(:,:) + tsurf_dev(:,:)

END SUBROUTINE build4PCM_tsurf
!=======================================================================

END MODULE surf_temp