! ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/soil.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $ ! SUBROUTINE soil(ptimestep, knon, ptsrf, ptsoil, s pcapcal, pfluxgrd) 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 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 ------------- use dimphy, only: klon IMPLICIT NONE include "YOMCST.h" include "dimsoil.h" include "clesphys.h" c----------------------------------------------------------------------- c arguments c --------- REAL, intent(IN) :: ptimestep INTEGER, intent(IN) :: knon REAL, intent(IN) :: ptsrf(klon) REAL, intent(OUT) :: ptsoil(klon,nsoilmx) REAL, intent(OUT) :: 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,SAVE :: dz1(nsoilmx),dz2(nsoilmx) REAL,allocatable,save :: zc(:,:),zd(:,:) REAL,SAVE :: lambda LOGICAL,SAVE :: firstcall=.true. c----------------------------------------------------------------------- c Depths: c ------- REAL fz,rk,fz1,rk1,rk2 fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) pfluxgrd(:) = 0. ! on Venus thermal inertia is assumed constant over the globe DO ig = 1, knon ztherm_i(ig) = inertie ENDDO IF (firstcall) THEN allocate(zc(klon,nsoilmx),zd(klon,nsoilmx)) c----------------------------------------------------------------------- c ground levels c grnd=z/l where l is the skin depth of the diurnal cycle: c -------------------------------------------------------- c VENUS : A REVOIR !!!! min_period=20000. ! in seconds dalph_soil=2. ! ratio between successive layer sizes 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 The first soil layer depth, based on min_period 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 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,knon 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,knon 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 firstcall =.false. ENDIF !--not firstcall c----------------------------------------------------------------------- c Computation of the soil temperatures using the Cgrd and Dgrd c coefficient computed at the previous time-step: c ----------------------------------------------- c temperature in the first soil layer DO ig=1,knon ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ s (lambda*(1.-zd(ig,1))+1.) ENDDO c temperatures in the other soil layers DO jk=1,nsoilmx-1 DO ig=1,knon ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) ENDDO ENDDO 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,knon 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,knon 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,knon 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 END