source: LMDZ6/branches/Amaury_dev/libf/phylmd/soil.F90 @ 5112

! $Header$

SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, &
        lon, lat, ptsoil, pcapcal, pfluxgrd)

  USE dimphy
  USE lmdz_phys_para
  USE indice_sol_mod
  USE lmdz_print_control, ONLY: lunout
  USE lmdz_abort_physic, ONLY: abort_physic

  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)

  !=======================================================================
  INCLUDE "YOMCST.h"
  INCLUDE "dimsoil.h"
  INCLUDE "comsoil.h"
  !-----------------------------------------------------------------------
  ! 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