! ! $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