source: LMDZ.3.3/trunk/libf/phylmd/soil.F @ 5415

      SUBROUTINE soil(ptimestep, indice, snow, ptsrf, ptsoil,
     s          pcapcal, pfluxgrd)
      IMPLICIT NONE

c=======================================================================
c
c   Auteur:  Frederic Hourdin     30/01/92
c   -------
c
c   objet:  computation of : the soil temperature evolution
c   ------                   the surfacic heat capacity "Capcal"
c                            the surface conduction flux pcapcal
c
c
c   Method: implicit time integration
c   -------
c   Consecutive ground temperatures are related by:
c           T(k+1) = C(k) + D(k)*T(k)  (1)
c   the coefficients C and D are computed at the t-dt time-step.
c   Routine structure:
c   1)new temperatures are computed  using (1)
c   2)C and D coefficients are computed from the new temperature
c     profile for the t+dt time-step
c   3)the coefficients A and B are computed where the diffusive
c     fluxes at the t+dt time-step is given by
c            Fdiff = A + B Ts(t+dt)
c     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
c            with F0 = A + B (Ts(t))
c                 Capcal = B*dt
c           
c   Interface:
c   ----------
c
c   Arguments:
c   ----------
c   ptimestep            physical timestep (s)
c   indice               sub-surface index
c   snow(klon,nbsrf)     snow
c   ptsrf(klon)          surface temperature at time-step t (K)
c   ptsoil(klon,nsoilmx) temperature inside the ground (K)
c   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
c   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
c   
c=======================================================================
c   declarations:
c   -------------

#include "dimensions.h"
#include "dimphy.h"
#include "dimsoil.h"
#include "indicesol.h"

c-----------------------------------------------------------------------
c  arguments
c  ---------

      REAL ptimestep
      INTEGER indice
      REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
      REAL pcapcal(klon),pfluxgrd(klon)

c-----------------------------------------------------------------------
c  local arrays
c  ------------

      INTEGER ig,jk
      REAL zdz2(nsoilmx),z1(klon)
      REAL min_period,dalph_soil
      REAL ztherm_i(klon)

c   local saved variables:
c   ----------------------
      REAL dz1(nsoilmx),dz2(nsoilmx)
      REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
      REAL lambda
      SAVE dz1,dz2,zc,zd,lambda
      LOGICAL firstcall
      SAVE firstcall
      REAL isol,isno,iice
      SAVE isol,isno,iice

      DATA firstcall/.true./

      DATA isol,isno,iice/2000.,2000.,2000./

c-----------------------------------------------------------------------
c   Depthts:
c   --------

      REAL fz,rk,fz1,rk1,rk2
      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)

c   calcul de l'inertie thermique a partir de la variable rnat.
c   on initialise a iice meme au-dessus d'un point de mer au cas 
c   ou le point de mer devienne point de glace au pas suivant
c   on corrige si on a un point de terre avec ou sans glace
c
      IF (indice.EQ.is_sic) THEN
         DO ig = 1, klon
            ztherm_i(ig)   = iice
            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
         ENDDO
      ELSE IF (indice.EQ.is_lic) THEN
         DO ig = 1, klon
            ztherm_i(ig)   = iice
            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
         ENDDO
      ELSE IF (indice.EQ.is_ter) THEN
         DO ig = 1, klon
            ztherm_i(ig)   = isol
            IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
         ENDDO
      ELSE IF (indice.EQ.is_oce) THEN
         DO ig = 1, klon
            ztherm_i(ig)   = iice
         ENDDO
      ELSE
         PRINT*, "valeur d indice non prevue", indice
         CALL abort
      ENDIF


      IF (firstcall) THEN

c-----------------------------------------------------------------------
c   ground levels 
c   grnd=z/l where l is the skin depth of the diurnal cycle:
c   --------------------------------------------------------

         min_period=1800. ! en secondes
         dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.

         OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
         READ(99,*) min_period
         READ(99,*) dalph_soil
         PRINT*,'Discretization for the soil model'
         PRINT*,'First level e-folding depth',min_period,
     s   '   dalph',dalph_soil
         CLOSE(99)
9999     CONTINUE

c   la premiere couche represente un dixieme de cycle diurne
         fz1=sqrt(min_period/3.14)

         DO jk=1,nsoilmx
            rk1=jk
            rk2=jk-1
            dz2(jk)=fz(rk1)-fz(rk2)
         ENDDO
         DO jk=1,nsoilmx-1
            rk1=jk+.5
            rk2=jk-.5
            dz1(jk)=1./(fz(rk1)-fz(rk2))
         ENDDO
         lambda=fz(.5)*dz1(1)
         PRINT*,'full layers, intermediate layers (seconds)'
         DO jk=1,nsoilmx
            rk=jk
            rk1=jk+.5
            rk2=jk-.5
            PRINT *,'fz=',
     .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
         ENDDO

         firstcall =.false.

c   Initialisations:
c   ----------------

      ELSE   !--not firstcall
c-----------------------------------------------------------------------
c   Computation of the soil temperatures using the Cgrd and Dgrd
c  coefficient computed at the previous time-step:
c  -----------------------------------------------

c    surface temperature
         DO ig=1,klon
            ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
     s      (lambda*(1.-zd(ig,1))+1.)
         ENDDO

c   other temperatures
         DO jk=1,nsoilmx-1
            DO ig=1,klon
               ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
            ENDDO
         ENDDO

      ENDIF !--not firstcall
c-----------------------------------------------------------------------
c   Computation of the Cgrd and Dgrd coefficient for the next step:
c   ---------------------------------------------------------------

      DO jk=1,nsoilmx
         zdz2(jk)=dz2(jk)/ptimestep
      ENDDO

      DO ig=1,klon
         z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1)
         zc(ig,nsoilmx-1)=zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig)
         zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig)
      ENDDO

      DO jk=nsoilmx-1,2,-1
         DO ig=1,klon
            z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
            zc(ig,jk-1)=
     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
            zd(ig,jk-1)=dz1(jk-1)*z1(ig)
         ENDDO
      ENDDO

c-----------------------------------------------------------------------
c   computation of the surface diffusive flux from ground and
c   calorific capacity of the ground:
c   ---------------------------------

      DO ig=1,klon
         pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*
     s   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
         pcapcal(ig)=ztherm_i(ig)*
     s   (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))
         z1(ig)=lambda*(1.-zd(ig,1))+1.
         pcapcal(ig)=pcapcal(ig)/z1(ig)
         pfluxgrd(ig)=pfluxgrd(ig)
     s   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
     s   /ptimestep
      ENDDO

      RETURN
      END