source: LMDZ4/branches/V3_test/libf/phylmd/soil.F @ 3651

!
! $Header$
!
      SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil,
     s          pcapcal, pfluxgrd)
      use dimphy
      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   -------------

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

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

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

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

      INTEGER ig,jk
c@$$      REAL zdz2(nsoilmx),z1(klon)
      REAL zdz2(nsoilmx),z1(klon,nbsrf)
      REAL,SAVE :: min_period,dalph_soil
      REAL ztherm_i(klon)

c   local saved variables:
c   ----------------------
      REAL dz1(nsoilmx),dz2(nsoilmx)
c@$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
cym      REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)
      REAL,ALLOCATABLE,SAVE ::  zc(:,:,:),zd(:,:,:)
c$OMP THREADPRIVATE(zc,zd)
      REAL lambda
cym      SAVE dz1,dz2,zc,zd,lambda
      SAVE dz1,dz2,lambda
c$OMP THREADPRIVATE(dz1,dz2,lambda)
      LOGICAL firstcall, firstsurf(nbsrf)
      SAVE firstcall, firstsurf
c$OMP THREADPRIVATE(firstcall, firstsurf)
      REAL isol,isno,iice
      SAVE isol,isno,iice
c$OMP THREADPRIVATE(isol,isno,iice)
      DATA firstcall/.true./
      DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./

      DATA isol,isno,iice/2000.,2000.,2000./
      LOGICAL,SAVE :: First=.true.
c$OMP THREADPRIVATE(First)
c-----------------------------------------------------------------------
c   Depthts:
c   --------

      REAL fz,rk,fz1,rk1,rk2
      fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
      pfluxgrd(:) = 0.
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 (first) THEN
        allocate(zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf))
        first=.false.
      ENDIF
      
      IF (indice.EQ.is_sic) THEN
         DO ig = 1, knon
            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, knon
            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, knon
            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, knon
            ztherm_i(ig)   = iice
         ENDDO
      ELSE
         PRINT*, "valeur d indice non prevue", indice
         CALL abort
      ENDIF


c@$$      IF (firstcall) THEN
      IF (firstsurf(indice)) 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.
c$OMP MASTER
         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$OMP END MASTER
c$OMP BARRIER

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
C PB
         firstsurf(indice) = .FALSE. 
c@$$         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,knon
            ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
     s      (lambda*(1.-zd(ig,1,indice))+1.)
         ENDDO

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

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

c@$$  PB ajout pour cas glace de mer
      IF (indice .EQ. is_sic) THEN
          DO ig = 1 , knon
            ptsoil(ig,nsoilmx) = RTT - 1.8
          END DO 
      ENDIF 

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

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

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

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

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

      RETURN
      END