! $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 USE lmdz_comsoil, ONLY: inertie_sol, inertie_sno, inertie_sic, inertie_lic, iflag_sic, iflag_inertie USE lmdz_yomcst 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 "dimsoil.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 lmdz_comsoil 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