! ! $Header$ ! SUBROUTINE soil_hetero(ptimestep, indice, knon, snow, ptsrf, qsol, & lon, lat, ptsoil, pcapcal, pfluxgrd, ztherm_i, conv_ratio) USE yomcst_mod_h, ONLY: RTT, RPI USE dimsoil_mod_h, ONLY: nsoilmx USE comsoil_mod_h USE compbl_mod_h USE dimpft_mod_h USE dimphy USE mod_phys_lmdz_para USE indice_sol_mod USE print_control_mod, ONLY: lunout USE phys_state_var_mod, ONLY: alpha_soil_tersrf, period_tersrf USE surf_param_mod, ONLY: eff_surf_param 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) ! 2025/04 : A. Maison, adapting the routine for heterogeneous continental sub-surfaces ! ! 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) ! ztherm_i(klon) soil thermal inertia (J.m-2.K.s-1/2) ! conv_ratio(klon) ratio to convert soil depths in meters (-) ! !======================================================================= !----------------------------------------------------------------------- ! 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), INTENT(IN) :: ztherm_i REAL, DIMENSION(klon), INTENT(IN) :: conv_ratio 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, DIMENSION(nsoilmx) :: zdz2 REAL :: z1s 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*(alpha_soil_tersrf**rk-1.)/(alpha_soil_tersrf-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: !----------------------------------------------------------------------- ! la premiere couche represente un dixieme de cycle diurne fz1=SQRT(period_tersrf/RPI) 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,*) 'surface index:', indice 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)*RPI,fz(rk)*fz(rk)*RPI ENDDO WRITE(lunout,*)'full layers, intermediate layers (meters)' DO jk=1,nsoilmx rk=jk rk2=jk-.5 WRITE(lunout,*)'fz=', & fz(rk2)*conv_ratio, fz(rk)*conv_ratio ENDDO firstcall =.FALSE. END IF !----------------------------------------------------------------------- ! 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_hetero