MODULE soil_pem_ini_mod

implicit none

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

SUBROUTINE soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)

use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
use comsoil_h,     only: volcapa

implicit none

!-----------------------------------------------------------------------
!  Author: LL
!  Purpose: Compute soil temperature using an implict 1st order scheme, stationnary solution
!
!  Note: depths of layers and mid-layers, soil thermal inertia and
!        heat capacity are commons in comsoil_PEM.h
!-----------------------------------------------------------------------

#include "dimensions.h"

!-----------------------------------------------------------------------
!  arguments
!  ---------
!  inputs:
integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
integer,                      intent(in) :: nsoil   ! number of soil layers
real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]

! outputs:
real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
! local variables:
integer :: ig, ik, iloop

! 0. Initialisations and preprocessing step
! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
do ig = 1,ngrid
    do ik = 0,nsoil - 1
        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
    enddo
enddo

! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
do ig = 1,ngrid
    do ik = 1,nsoil - 1
        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
    enddo
enddo

! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
! mu_PEM
mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))

! q_{1/2}
coefq_PEM(:) = 0.
! q_{k+1/2}

do ig = 1,ngrid
    ! d_k
    do ik = 1,nsoil-1
        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
    enddo

    ! alph_PEM_{N-1}
    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
    ! alph_PEM_k
    do ik = nsoil - 2,1,-1
        alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
    enddo
enddo ! of do ig=1,ngrid

!  1. Compute beta_PEM coefficients
! Bottom layer, beta_PEM_{N-1}
do ig = 1,ngrid
        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
enddo
! Other layers
do ik = nsoil - 2,1,-1
    do ig = 1,ngrid
        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
    enddo
enddo

!  2. Compute soil temperatures
do iloop = 1,10 !just convergence
    ! First layer:
    do ig = 1,ngrid
        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
    ! Other layers:
    do ik = 1,nsoil - 1
        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
    enddo
enddo

enddo ! iloop

END SUBROUTINE soil_pem_ini

END MODULE soil_pem_ini_mod