source: trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini.F90 @ 2909

      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,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia [SI]
      real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]

! outputs:
      real,intent(inout) :: tsoil(ngrid,nsoil) ! 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.0-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.0-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