source: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90 @ 3578

MODULE compute_soiltemp_mod

implicit none
!-----------------------------------------------------------------------
!  Author: LL
!  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation.
!
!  Note: depths of layers and mid-layers, soil thermal inertia and
!        heat capacity are commons in comsoil_PEM.h
!-----------------------------------------------------------------------
contains
!=======================================================================

SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,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
!
!  Note: depths of layers and mid-layers, soil thermal inertia and
!        heat capacity are commons in comsoil_PEM.h
!-----------------------------------------------------------------------

#include "dimensions.h"

! Inputs:
! -------
integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
integer,                      intent(in) :: nsoil     ! number of soil layers
logical,                      intent(in) :: firstcall ! identifier for initialization call
real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
real,                         intent(in) :: timestep  ! time step [s]
real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
! Outputs:
!---------
real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
! Local:
!-------
integer :: ig, ik

! 0. Initialisations and preprocessing step
 if (firstcall) then
! 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_PEM
    ! mu_PEM
    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))

    ! q_{1/2}
    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
    ! q_{k+1/2}
    do ik = 1,nsoil - 1
        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
    enddo

    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
endif ! of if (firstcall)

if (.not. firstcall) THEN
! 2. Compute soil temperatures
! 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
endif

!  2. Compute beta_PEM coefficients (preprocessing for next time step)
! 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

END SUBROUTINE compute_tsoil_pem

!=======================================================================

SUBROUTINE ini_tsoil_pem(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: Initialize the soil with the solution of the stationnary problem of Heat Conduction. Boundarry conditions: Tsurf averaged from the PCM; Geothermal flux at the bottom layer
!
!  Note: depths of layers and mid-layers, soil thermal inertia and
!        heat capacity are commons in comsoil_PEM.h
!-----------------------------------------------------------------------

#include "dimensions.h"

! 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:
!-------
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_PEM
! 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
    do ig = 1,ngrid
        ! First layer:
        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 ini_tsoil_pem

!=======================================================================

SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil)

use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
use math_mod,      only: solve_steady_heat

implicit none

!-----------------------------------------------------------------------
!  Author: JBC
!  Purpose: Shifting the soil temperature profile to follow the surface evolution due to ice condensation/sublimation
!-----------------------------------------------------------------------
! Inputs:
! -------
integer,                       intent(in) :: ngrid       ! number of (horizontal) grid-points
integer,                       intent(in) :: nsoil       ! number of soil layers
integer,                       intent(in) :: nslope      ! number of sub-slopes
real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m]
real, dimension(ngrid,nslope), intent(in) :: zlag        ! newly built lag thickness [m]
real, dimension(ngrid,nslope), intent(in) :: tsurf       ! surface temperature [K]
! Outputs:
! --------
real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
! Local:
! ------
integer                             :: ig, isoil, islope, ishift, ilag, iz
real                                :: a, z, zshift_surfloc
real, dimension(ngrid,nsoil,nslope) :: tsoil_old

write(*,*)"Shifting soil temperature to surface"
tsoil_old = tsoil

do ig = 1,ngrid
    do islope = 1,nslope
        zshift_surfloc = zshift_surf(ig,islope)

        if (zshift_surfloc >= 0.) then ! In case of the surface is higher than initially
            if (zshift_surfloc < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
                ! Nothing to do; we keep the soil temperature profile
            else if (zshift_surfloc >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
                tsoil(ig,:,islope) = tsurf(ig,islope)
            else
                ishift = 0
                do while (mlayer_PEM(ishift) <= zshift_surfloc)
                    ishift = ishift + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
                enddo
                ! The "new soil" temperature is set to tsurf
                tsoil(ig,:ishift,islope) = tsurf(ig,islope)

                do isoil = ishift + 1,nsoil
                    ! Position in the old discretization of the depth
                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
                    iz = 0
                    do while (mlayer_PEM(iz) <= z)
                        iz = iz + 1
                    enddo
                    ! Interpolation of the temperature profile from the old discretization
                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
                enddo
            endif
        else ! In case of the surface is lower than initially
            if (abs(zshift_surfloc) < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
                ! Nothing to do; we keep the soil temperature profile
            else if (abs(zshift_surfloc) >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
                call solve_steady_heat(nsoil,mlayer_PEM,layer_PEM,mthermdiff_PEM(ig,:),thermdiff_PEM(ig,:),tsurf(ig,islope),fluxgeo,tsoil(ig,:,islope))
            else
                if (zlag(ig,islope) < mlayer_PEM(0)) then ! The lag is too thin to be taken into account
                    ilag = 0
                else
                    ilag = 0
                    do while (mlayer_PEM(ilag) <= zlag(ig,islope))
                        ilag = ilag + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
                    enddo
                    ! Position of the lag bottom in the old discretization of the depth
                    z = zlag(ig,islope) - zshift_surfloc
                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
                    iz = 0
                    do while (mlayer_PEM(iz) <= z)
                        iz = iz + 1
                    enddo
                    ! The "new lag" temperature is set to the ice temperature just below
                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
                    tsoil(ig,:ilag,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
                endif

                ishift = nsoil - 1
                z = mlayer_PEM(nsoil - 1) + zshift_surfloc
                do while (mlayer_PEM(ishift) >= z)
                    ishift = ishift - 1
                enddo
                ishift = ishift + 1 ! Adding 1 is needed to match the good index for tsoil!
                do isoil = ilag + 1,ishift
                    ! Position in the old discretization of the depth
                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
                    iz = 0
                    do while (mlayer_PEM(iz) <= z)
                        iz = iz + 1
                    enddo
                    ! Interpolation of the temperature profile from the old discretization
                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
                enddo

                ! The "new deepest layers" temperature is set by solving the steady heat equation
                call solve_steady_heat(nsoil - ishift + 1,mlayer_PEM(ishift - 1:),layer_PEM(ishift:),mthermdiff_PEM(ig,ishift - 1:),thermdiff_PEM(ig,ishift:),tsoil(ig,ishift,islope),fluxgeo,tsoil(ig,ishift:,islope))
            endif
        endif
    enddo
enddo

END SUBROUTINE shift_tsoil2surf

END MODULE compute_soiltemp_mod