MODULE surf_temp

implicit none

!=======================================================================
contains
!=======================================================================

SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)

use grid_conversion, only: vect2lonlat, lonlat2vect

implicit none

! Inputs:
integer,                          intent(in) :: nlon, nlat, nslope, ngrid
real,    dimension(ngrid,nslope), intent(in) :: co2_ice
real,    dimension(ngrid),        intent(in) :: latitude
logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_ini
! Outputs:
real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg
logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared
! 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

! 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)
! Priority given to equator-ward slope which are most likely to hold no ice

implicit none

! Inputs:
real,    intent(in) :: lat
integer, intent(in) :: nslope, islope
! Outputs:
integer, dimension(nslope - 1), intent(out) :: priority
! Locals:
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