SUBROUTINE soil(ptimestep, indice, snow, ptsrf, ptsoil, s pcapcal, pfluxgrd) IMPLICIT NONE c======================================================================= c c Auteur: Frederic Hourdin 30/01/92 c ------- c c objet: computation of : the soil temperature evolution c ------ the surfacic heat capacity "Capcal" c the surface conduction flux pcapcal c c c Method: implicit time integration c ------- c Consecutive ground temperatures are related by: c T(k+1) = C(k) + D(k)*T(k) (1) c the coefficients C and D are computed at the t-dt time-step. c Routine structure: c 1)new temperatures are computed using (1) c 2)C and D coefficients are computed from the new temperature c profile for the t+dt time-step c 3)the coefficients A and B are computed where the diffusive c fluxes at the t+dt time-step is given by c Fdiff = A + B Ts(t+dt) c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt c with F0 = A + B (Ts(t)) c Capcal = B*dt c c Interface: c ---------- c c Arguments: c ---------- c ptimestep physical timestep (s) c indice sub-surface index c snow(klon,nbsrf) snow c ptsrf(klon) surface temperature at time-step t (K) c ptsoil(klon,nsoilmx) temperature inside the ground (K) c pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) c pfluxgrd(klon) surface diffusive flux from ground (Wm-2) c c======================================================================= c declarations: c ------------- #include "dimensions.h" #include "dimphy.h" #include "dimsoil.h" #include "indicesol.h" c----------------------------------------------------------------------- c arguments c --------- REAL ptimestep INTEGER indice REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon) REAL pcapcal(klon),pfluxgrd(klon) c----------------------------------------------------------------------- c local arrays c ------------ INTEGER ig,jk REAL zdz2(nsoilmx),z1(klon) REAL min_period,dalph_soil REAL ztherm_i(klon) c local saved variables: c ---------------------- REAL dz1(nsoilmx),dz2(nsoilmx) REAL zc(klon,nsoilmx),zd(klon,nsoilmx) REAL lambda SAVE dz1,dz2,zc,zd,lambda LOGICAL firstcall SAVE firstcall REAL isol,isno,iice SAVE isol,isno,iice DATA firstcall/.true./ DATA isol,isno,iice/2000.,2000.,2000./ c----------------------------------------------------------------------- c Depthts: c -------- REAL fz,rk,fz1,rk1,rk2 fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) c calcul de l'inertie thermique a partir de la variable rnat. c on initialise a iice meme au-dessus d'un point de mer au cas c ou le point de mer devienne point de glace au pas suivant c on corrige si on a un point de terre avec ou sans glace c IF (indice.EQ.is_sic) THEN DO ig = 1, klon ztherm_i(ig) = iice IF (snow(ig).GT.0.0) ztherm_i(ig) = isno ENDDO ELSE IF (indice.EQ.is_lic) THEN DO ig = 1, klon ztherm_i(ig) = iice IF (snow(ig).GT.0.0) ztherm_i(ig) = isno ENDDO ELSE IF (indice.EQ.is_ter) THEN DO ig = 1, klon ztherm_i(ig) = isol IF (snow(ig).GT.0.0) ztherm_i(ig) = isno ENDDO ELSE IF (indice.EQ.is_oce) THEN DO ig = 1, klon ztherm_i(ig) = iice ENDDO ELSE PRINT*, "valeur d indice non prevue", indice CALL abort ENDIF IF (firstcall) THEN c----------------------------------------------------------------------- c ground levels c grnd=z/l where l is the skin depth of the diurnal cycle: c -------------------------------------------------------- min_period=1800. ! en secondes dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ. OPEN(99,file='soil.def',status='old',form='formatted',err=9999) READ(99,*) min_period READ(99,*) dalph_soil PRINT*,'Discretization for the soil model' PRINT*,'First level e-folding depth',min_period, s ' dalph',dalph_soil CLOSE(99) 9999 CONTINUE c 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) PRINT*,'full layers, intermediate layers (seconds)' DO jk=1,nsoilmx rk=jk rk1=jk+.5 rk2=jk-.5 PRINT *,'fz=', . fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 ENDDO firstcall =.false. c Initialisations: c ---------------- ELSE !--not firstcall c----------------------------------------------------------------------- c Computation of the soil temperatures using the Cgrd and Dgrd c coefficient computed at the previous time-step: c ----------------------------------------------- c surface temperature DO ig=1,klon ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ s (lambda*(1.-zd(ig,1))+1.) ENDDO c other temperatures DO jk=1,nsoilmx-1 DO ig=1,klon ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) ENDDO ENDDO ENDIF !--not firstcall c----------------------------------------------------------------------- c Computation of the Cgrd and Dgrd coefficient for the next step: c --------------------------------------------------------------- DO jk=1,nsoilmx zdz2(jk)=dz2(jk)/ptimestep ENDDO DO ig=1,klon z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1) zc(ig,nsoilmx-1)=zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig) zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig) ENDDO DO jk=nsoilmx-1,2,-1 DO ig=1,klon z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) zc(ig,jk-1)= s (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig) zd(ig,jk-1)=dz1(jk-1)*z1(ig) ENDDO ENDDO c----------------------------------------------------------------------- c computation of the surface diffusive flux from ground and c calorific capacity of the ground: c --------------------------------- DO ig=1,klon pfluxgrd(ig)=ztherm_i(ig)*dz1(1)* s (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) pcapcal(ig)=ztherm_i(ig)* s (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1)) z1(ig)=lambda*(1.-zd(ig,1))+1. pcapcal(ig)=pcapcal(ig)/z1(ig) pfluxgrd(ig)=pfluxgrd(ig) s +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) s /ptimestep ENDDO RETURN END