MODULE surf_temp

implicit none

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

SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,iim_input,jjm_input,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)

implicit none

! Inputs:
integer,                          intent(in) :: iim_input, jjm_input, nslope, ngrid
real,    dimension(ngrid,nslope), intent(in) :: co2_ice
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(iim_input + 1,jjm_input + 1,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
real, dimension(ngrid)                              :: tmp

write(*,*) "> Updating surface temperature where ice disappeared"

! Convert from reduced grid to lon-lat grid
#ifndef CPP_1D
do islope = 1,nslope
    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,co2_ice(:,islope),co2ice_ll(:,:,islope))
    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
enddo
#else
    tsurf_ll(1,1,:) = tsurf_avg(1,:)
    co2ice_ll(1,1,:) = co2_ice(1,:)
    mask_co2ice_ini(1,1,:) = merge(1.,0.,is_co2ice_ini(1,:))
    co2ice_disappeared_ll(1,1,:) = merge(1.,0.,co2ice_disappeared(1,:))
#endif

! For each point where ice disappeared
rmax = max(iim_input + 1,jjm_input + 1)
do j = 1,jjm_input + 1
    do i = 1,iim_input + 1
        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.
                do k = 1,nslope
                    if (k /= islope .and. mask_co2ice_ini(i,j,k) < 0.5) then
                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,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
                                if (ii >= 1 .and. ii <= iim_input + 1 .and. jj >= 1 .and. jj <= jjm_input + 1) then
                                    do k = 1,nslope
                                        if (mask_co2ice_ini(ii,jj,k) < 0.5) then
                                            tsurf_ll(i,j,islope) = tsurf_ll(i,j,k)
                                            found = .true.
                                            exit
                                        endif
                                    enddo
                                endif
                            endif
                            if (found) exit
                        enddo
                        if (found) exit
                    enddo
                    radius = radius + 1
                enddo
            endif

        enddo
    enddo
enddo

! Convert back from lon-lat grid to reduced grid
#ifndef CPP_1D
do islope = 1,nslope
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
    where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
enddo
#else
    tsurf_avg(1,:) = tsurf_ll(1,1,:)
    where(co2ice_disappeared_ll(1,1,:) > 0.5) co2ice_disappeared(1,:) = .true.
#endif

END SUBROUTINE update_tsurf_nearest_baresoil

END MODULE surf_temp