MODULE soil_TIfeedback_PEM_mod

implicit none

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

SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,newtherm_i)

use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM

implicit none

!=======================================================================
!   Description :
!       Surface water ice / Thermal inertia feedback.
!
!   When surface water-ice is thick enough, this routine creates a new
!   soil thermal inertia with three different layers :
!   - One layer of surface water ice (the thickness is given
!     by the variable icecover (in kg of ice per m2) and the thermal
!     inertia is prescribed by inert_h2o_ice (see surfdat_h));
!   - A transitional layer of mixed thermal inertia;
!   - A last layer of regolith below the ice cover whose thermal inertia
!     is equal to inertiedat.
!
!  To use the model :
!      SET THE tifeedback LOGICAL TO ".true." in callphys.def.
!
!  Author: Adapted from J.-B. Madeleine Mars 2008 ( Updated November 2012) by LL, 2022
!=======================================================================

!Inputs
!------
integer,                intent(in) :: ngrid    ! Number of horizontal grid points
integer,                intent(in) :: nsoil    ! Number of soil layers
real, dimension(ngrid), intent(in) :: icecover ! tracer on the surface (kg.m-2)
!Outputs
!-------
real,intent(inout) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
!Local variables
!---------------
integer :: ig                       ! Grid point (ngrid)
integer :: ik                       ! Grid point (nsoil)
integer :: iref                     ! Ice/Regolith boundary index
real    :: icedepth                 ! Ice cover thickness (m)
real    :: inert_h2o_ice = 800.     ! surface water ice thermal inertia [SI]
real    :: rho_ice = 920.           ! density of water ice [kg/m^3]
real    :: prev_thermi(ngrid,nsoil) ! previous thermal inertia [SI]

prev_thermi(:,:) = newtherm_i(:,:)

!Creating the new soil thermal inertia table
!-------------------------------------------
do ig = 1,ngrid
    ! Calculating the ice cover thickness
    icedepth = icecover(ig)/rho_ice

    ! If the ice cover is too thick or watercaptag = true, the entire column is changed:
    if (icedepth >= layer_PEM(nsoil)) then
        do ik = 1,nsoil
            newtherm_i(ig,ik) = max(inert_h2o_ice,prev_thermi(ig,ik))
        enddo
    ! We neglect the effect of a very thin ice cover:
    else if (icedepth < layer_PEM(1)) then
        do ik = 1,nsoil
            newtherm_i(ig,ik) = inertiedat_PEM(ig,ik)
        enddo
    else
        ! Ice/regolith boundary index:
        iref = 1
        ! Otherwise, we find the ice/regolith boundary:
        do ik=1,nsoil-1
            if ((icedepth >= layer_PEM(ik)) .and. (icedepth < layer_PEM(ik + 1))) then
                iref = ik + 1
                exit
            endif
        enddo
        ! And we change the thermal inertia:
        do ik = 1,iref - 1
            newtherm_i(ig,ik) = max(inert_h2o_ice,prev_thermi(ig,ik))
        enddo
        ! Transition (based on the equations of thermal conduction):
        newtherm_i(ig,iref) = sqrt((layer_PEM(iref) - layer_PEM(iref-1))/(((icedepth - layer_PEM(iref - 1))/newtherm_i(ig,iref - 1)**2) &
                              + ((layer_PEM(iref) - icedepth)/inertiedat_PEM(ig,ik)**2) ) )
        ! Underlying regolith:
        do ik = iref + 1,nsoil
            newtherm_i(ig,ik) = inertiedat_PEM(ig,ik)
        enddo
    endif ! icedepth
enddo ! ig

return

END SUBROUTINE soil_TIfeedback_PEM

END MODULE soil_TIfeedback_PEM_mod