source: LMDZ6/trunk/libf/phylmd/soil.f90 @ 5435

!
! $Header$
!
SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, &
     lon, lat, ptsoil, pcapcal, pfluxgrd)
  
  USE yomcst_mod_h
  USE dimphy
  USE mod_phys_lmdz_para
  USE indice_sol_mod
  USE print_control_mod, ONLY: lunout
  USE dimsoil_mod_h, ONLY: nsoilmx
  USE comsoil_mod_h
  IMPLICIT NONE

!=======================================================================
!
!   Auteur:  Frederic Hourdin     30/01/92
!   -------
!
!   Object:  Computation of : the soil temperature evolution
!   -------                   the surfacic heat capacity "Capcal"
!                            the surface conduction flux pcapcal
!
!   Update: 2021/07 : soil thermal inertia, formerly a constant value,
!   ------   can also be now a function of soil moisture (F Cheruy's idea)
!            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
!            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
!
!   Method: Implicit time integration
!   -------
!   Consecutive ground temperatures are related by:
!           T(k+1) = C(k) + D(k)*T(k)  (*)
!   The coefficients C and D are computed at the t-dt time-step.
!   Routine structure:
!   1) C and D coefficients are computed from the old temperature
!   2) new temperatures are computed using (*)
!   3) C and D coefficients are computed from the new temperature
!      profile for the t+dt time-step
!   4) the coefficients A and B are computed where the diffusive
!      fluxes at the t+dt time-step is given by
!             Fdiff = A + B Ts(t+dt)
!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
!             with F0 = A + B (Ts(t))
!                 Capcal = B*dt
!
!   Interface:
!   ----------
!
!   Arguments:
!   ----------
!   ptimestep            physical timestep (s)
!   indice               sub-surface index
!   snow(klon)           snow
!   ptsrf(klon)          surface temperature at time-step t (K)
!   qsol(klon)           soil moisture (kg/m2 or mm)
!   lon(klon)            longitude in radian
!   lat(klon)            latitude in radian
!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
!
!=======================================================================
! Arguments
! ---------
  REAL, INTENT(IN)                     :: ptimestep
  INTEGER, INTENT(IN)                  :: indice, knon !, knindex
  REAL, DIMENSION(klon), INTENT(IN)    :: snow
  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
  REAL, DIMENSION(klon), INTENT(IN)    :: qsol
  REAL, DIMENSION(klon), INTENT(IN)    :: lon
  REAL, DIMENSION(klon), INTENT(IN)    :: lat

  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd

!-----------------------------------------------------------------------
! Local variables
! ---------------
  INTEGER                             :: ig, jk, ierr
  REAL                                :: min_period,dalph_soil
  REAL, DIMENSION(nsoilmx)            :: zdz2
  REAL                                :: z1s
  REAL, DIMENSION(klon)               :: ztherm_i
  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef

! Local saved variables
! ---------------------
  REAL, SAVE                     :: lambda
!$OMP THREADPRIVATE(lambda)
  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
!$OMP THREADPRIVATE(dz1,dz2)
  LOGICAL, SAVE                  :: firstcall=.TRUE.
!$OMP THREADPRIVATE(firstcall)
    
!-----------------------------------------------------------------------
!   Depthts:
!   --------
  REAL fz,rk,fz1,rk1,rk2
  fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)


!-----------------------------------------------------------------------
! Calculation of some constants
! NB! These constants do not depend on the sub-surfaces
!-----------------------------------------------------------------------

  IF (firstcall) THEN
!-----------------------------------------------------------------------
!   ground levels 
!   grnd=z/l where l is the skin depth of the diurnal cycle:
!-----------------------------------------------------------------------

     min_period=1800. ! en secondes
     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
!$OMP MASTER
     IF (is_mpi_root) THEN
        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
        IF (ierr == 0) THEN ! Read file only if it exists
           READ(99,*) min_period
           READ(99,*) dalph_soil
           WRITE(lunout,*)'Discretization for the soil model'
           WRITE(lunout,*)'First level e-folding depth',min_period, &
                '   dalph',dalph_soil
           CLOSE(99)
        END IF
     ENDIF
!$OMP END MASTER
     CALL bcast(min_period)
     CALL bcast(dalph_soil)

!   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)
     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
     DO jk=1,nsoilmx
        rk=jk
        rk1=jk+.5
        rk2=jk-.5
        WRITE(lunout,*)'fz=', &
             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
     ENDDO

     firstcall =.FALSE.
  END IF


!-----------------------------------------------------------------------
!   Calcul de l'inertie thermique a partir de la variable rnat.
!   on initialise a inertie_sic meme au-dessus d'un point de mer au cas 
!   ou le point de mer devienne point de glace au pas suivant
!   on corrige si on a un point de terre avec ou sans glace
!
!   iophys can be used to write the ztherm_i variable in a phys.nc file
!   and check the results; to do so, add "CALL iophys_ini" in physiq_mod
!   and add knindex to the list of inputs in all the calls to soil.F90
!   (and to soil.F90 itself !)
!-----------------------------------------------------------------------

  IF (indice == is_sic) THEN
     DO ig = 1, knon
        ztherm_i(ig)   = inertie_sic
     ENDDO
     IF (iflag_sic == 0) THEN
       DO ig = 1, knon
         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
       ENDDO
!      Otherwise sea-ice keeps the same inertia, even when covered by snow
     ENDIF
!    CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', &
!      knon, knindex, ztherm_i)
  ELSE IF (indice == is_lic) THEN
     DO ig = 1, knon
        ztherm_i(ig)   = inertie_lic
        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
     ENDDO
!    CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', &
!      knon, knindex, ztherm_i)
  ELSE IF (indice == is_ter) THEN
     ! 
     ! La relation entre l'inertie thermique du sol et qsol change d'apres 
     !   iflag_inertie, defini dans physiq.def, et appele via comsoil.h
     ! 
     DO ig = 1, knon
        ! iflag_inertie=0 correspond au cas inertie=constant, comme avant
        IF (iflag_inertie==0) THEN         
           ztherm_i(ig)   = inertie_sol
        ELSE IF (iflag_inertie == 1) THEN
          ! I = a_qsol * qsol + b  modele lineaire deduit d'une
          ! regression lineaire I = a_mrsos * mrsos + b obtenue sur 
          ! sorties MO d'une simulation LMDZOR(CMIP6) sur l'annee 2000 
          ! sur tous les points avec frac_snow=0 
          ! Difference entre qsol et mrsos prise en compte par un
          ! facteur d'echelle sur le coefficient directeur de regression: 
          ! fact = 35./150. = mrsos_max/qsol_max
          ! et a_qsol = a_mrsos * fact (car a = dI/dHumidite)
            ztherm_i(ig) = 30.0 *35.0/150.0 *qsol(ig) +770.0
          ! AS : pour qsol entre 0 - 150, on a I entre 770 - 1820
        ELSE IF (iflag_inertie == 2) THEN
          ! deux regressions lineaires, sur les memes sorties,  
          ! distinguant le type de sol : sable ou autre (limons/argile)
          ! Implementation simple : regression type "sable" seulement pour 
          ! Sahara, defini par une "boite" lat/lon (NB : en radians !! )
          IF (lon(ig)>-0.35 .AND. lon(ig)<0.70 .AND. lat(ig)>0.17 .AND. lat(ig)<0.52) THEN 
              ! Valeurs theoriquement entre 728 et 2373 ; qsol valeurs basses
              ztherm_i(ig) = 47. *35.0/150.0 *qsol(ig) +728.  ! boite type "sable" pour Sahara
          ELSE 
              ! Valeurs theoriquement entre 550 et 1940 ; qsol valeurs moyennes et hautes
              ztherm_i(ig) = 41. *35.0/150.0 *qsol(ig) +505.
          ENDIF
        ELSE IF (iflag_inertie == 3) THEN
          ! AS : idee a tester : 
          ! si la relation doit etre une droite, 
          ! definissons-la en fonction des valeurs min et max de qsol (0:150),
          ! et de l'inertie (900 : 2000 ou 2400 ; choix ici: 2000)
          ! I = I_min + qsol * (I_max - I_min)/(qsol_max - qsol_min)
              ztherm_i(ig) = 900. + qsol(ig) * (2000. - 900.)/150.
        ELSE          
          WRITE (lunout,*) "Le choix iflag_inertie = ",iflag_inertie," n'est pas defini. Veuillez choisir un entier entre 0 et 3" 
        ENDIF
     !
     ! Fin de l'introduction de la relation entre l'inertie thermique du sol et qsol
     !-------------------------------------------
        !AS : donc le moindre flocon de neige sur un point de grid
        ! fait que l'inertie du point passe a la valeur pour neige ! 
        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
       
     ENDDO
!    CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', &
!      knon, knindex, ztherm_i)
  ELSE IF (indice == is_oce) THEN
     DO ig = 1, knon
!       This is just in case, but SST should be used by the model anyway
        ztherm_i(ig)   = inertie_sic
     ENDDO
!    CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', &
!      knon, knindex, ztherm_i)
  ELSE
     WRITE(lunout,*) "valeur d indice non prevue", indice
     call abort_physic("soil", "", 1)
  ENDIF


!-----------------------------------------------------------------------
! 1)
! Calculation of Cgrf and Dgrd coefficients using soil temperature from 
! previous time step.
!
! These variables are recalculated on the local compressed grid instead 
! of saved in restart file.
!-----------------------------------------------------------------------
  DO jk=1,nsoilmx
     zdz2(jk)=dz2(jk)/ptimestep
  ENDDO
  
  DO ig=1,knon
     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
     C_coef(ig,nsoilmx-1,indice)= &
          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
  ENDDO
  
  DO jk=nsoilmx-1,2,-1
     DO ig=1,knon
        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
             *(1.-D_coef(ig,jk,indice)))
        C_coef(ig,jk-1,indice)= &
             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
     ENDDO
  ENDDO

!-----------------------------------------------------------------------
! 2)
! Computation of the soil temperatures using the Cgrd and Dgrd
! coefficient computed above
!
!-----------------------------------------------------------------------

!    Surface temperature
  DO ig=1,knon
     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
          (lambda*(1.-D_coef(ig,1,indice))+1.)
  ENDDO
  
!   Other temperatures
  DO jk=1,nsoilmx-1
     DO ig=1,knon
        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
             *ptsoil(ig,jk)
     ENDDO
  ENDDO

  IF (indice == is_sic) THEN
     DO ig = 1 , knon
        ptsoil(ig,nsoilmx) = RTT - 1.8
     END DO
  ENDIF

!-----------------------------------------------------------------------
! 3)
! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil 
! temperature
!-----------------------------------------------------------------------
  DO ig=1,knon
     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
  ENDDO
  
  DO jk=nsoilmx-1,2,-1
     DO ig=1,knon
        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
             *(1.-D_coef(ig,jk,indice)))
        C_coef(ig,jk-1,indice) = &
             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
     ENDDO
  ENDDO

!-----------------------------------------------------------------------
! 4)
! Computation of the surface diffusive flux from ground and
! calorific capacity of the ground
!-----------------------------------------------------------------------
  DO ig=1,knon
     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
     pcapcal(ig)  = ztherm_i(ig)* &
          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
     pcapcal(ig)  = pcapcal(ig)/z1s
     pfluxgrd(ig) = pfluxgrd(ig) &
          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
          - lambda * C_coef(ig,1,indice) &
          - ptsrf(ig)) &
          /ptimestep
  ENDDO
    
END SUBROUTINE soil