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

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

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

!=======================================================================
SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
!-----------------------------------------------------------------------
! NAME
!     update_tsurf_nearest_baresoil
!
! DESCRIPTION
!     Update surface temperature where CO2 ice disappeared using nearby
!     bare soil temperature.
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use grid_conversion, only: vect2lonlat, lonlat2vect

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

! ARGUMENTS
! ---------
integer,                          intent(in)    :: nlon, nlat, nslope, ngrid ! Grid dimensions
real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                   ! CO2 ice density
real,    dimension(ngrid),        intent(in)    :: latitude                  ! Latitude
logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_ini             ! Initial CO2 ice flag
real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg                 ! Average surface temperature
logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared        ! Ice disappeared flag

! LOCAL VARIABLES
! ---------------
real, parameter                      :: eps = 1.e-10
integer                              :: islope, i, j, k, radius, rmax, di, dj, ii, jj
logical                              :: found
real,    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
real,    dimension(nlon,nlat)        :: latitude_ll
real,    dimension(ngrid)            :: tmp
integer, dimension(nslope - 1)       :: priority

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

write(*,*) "> Updating surface temperature where ice disappeared"
! Convert from reduced grid to lon-lat grid
call vect2lonlat(nlon,nlat,ngrid,latitude,latitude_ll)
do islope = 1,nslope
    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,co2_ice(:,islope),co2ice_ll(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
enddo

! 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 .and. co2ice_ll(i,j,islope) < eps .and. co2ice_disappeared_ll(i,j,islope) < 0.5) then
                found = .false.
                co2ice_disappeared_ll(i,j,islope) = 1.
                call get_slope_priority(latitude_ll(i,j),nslope,islope,priority)
                do k = 1,nslope - 1
                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5) then
                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
                        found = .true.
                        exit
                    endif
                enddo

                radius = 1
                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
                    do dj = -radius,radius
                        do di = -radius,radius
                            if (abs(di) + abs(dj) == radius) then
                                ii = i + di
                                jj = j + dj
                                ! Longitudinal periodicity
                                if (ii < 1) then
                                    ii = ii + nlon
                                else if (ii > nlon) then
                                    ii = ii - nlon
                                endif
                                ! Latitude boundaries
                                if (jj >= 1 .and. jj <= nlat) then
                                    call get_slope_priority(latitude_ll(ii,jj),nslope,islope,priority)
                                    do k = 1,nslope - 1
                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5) then
                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
                                            found = .true.
                                            exit
                                        endif
                                    enddo
                                endif
                            endif
                            if (found) exit
                        enddo
                        if (found) exit
                    enddo
                    radius = radius + 1
                enddo
                if (.not. found) write(*,*) "WARNING: no bare soil found for ice disappeared on  point:",i,j,islope
            endif
        enddo
    enddo
enddo

! 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,co2ice_disappeared_ll(:,:,islope),tmp)
    where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
enddo

END SUBROUTINE update_tsurf_nearest_baresoil
!=======================================================================

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

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

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

! LOCAL VARIABLES
! ---------------
integer :: 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
    enddo
    ! Pole-ward slopes
    do i = islope + 1,nslope
        priority(k) = i
        k = k + 1
    enddo
else ! Southern hemisphere
    ! Equator-ward slopes
    do i = islope + 1,nslope
        priority(k) = i
        k = k + 1
    enddo
    ! Pole-ward slopes
    do i = islope - 1,1,-1
        priority(k) = i
        k = k + 1
    enddo
endif

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

END MODULE surf_temp