/[lmdze]/trunk/Sources/phylmd/Interface_surf/soil.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Interface_surf/soil.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/soil.f revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/Sources/phylmd/Interface_surf/soil.f revision 157 by guez, Mon Jul 20 16:01:49 2015 UTC
# Line 1  Line 1 
1  !  module soil_m
 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/soil.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $  
 !  
       SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil,  
      s          pcapcal, pfluxgrd)  
       use dimens_m  
       use indicesol  
       use dimphy  
       use dimsoil  
       use SUPHEC_M  
       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   -------------  
   
   
 c-----------------------------------------------------------------------  
 c  arguments  
 c  ---------  
   
       REAL ptimestep  
       INTEGER indice, knon  
       REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)  
       REAL pcapcal(klon),pfluxgrd(klon)  
   
 c-----------------------------------------------------------------------  
 c  local arrays  
 c  ------------  
   
       INTEGER ig,jk  
 c$$$      REAL zdz2(nsoilmx),z1(klon)  
       REAL zdz2(nsoilmx),z1(klon,nbsrf)  
       REAL min_period,dalph_soil  
       REAL ztherm_i(klon)  
   
 c   local saved variables:  
 c   ----------------------  
       REAL dz1(nsoilmx),dz2(nsoilmx)  
 c$$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)  
       REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)  
       REAL lambda  
       SAVE dz1,dz2,zc,zd,lambda  
       LOGICAL firstcall, firstsurf(nbsrf)  
       SAVE firstcall, firstsurf  
       REAL isol,isno,iice  
       SAVE isol,isno,iice  
   
       DATA firstcall/.true./  
       DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.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.)  
       pfluxgrd(:) = 0.  
 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, knon  
             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, knon  
             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, knon  
             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, knon  
             ztherm_i(ig)   = iice  
          ENDDO  
       ELSE  
          PRINT*, "valeur d indice non prevue", indice  
          stop 1  
       ENDIF  
   
   
 c$$$      IF (firstcall) THEN  
       IF (firstsurf(indice)) 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  
 C PB  
          firstsurf(indice) = .FALSE.  
 c$$$         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,knon  
             ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/  
      s      (lambda*(1.-zd(ig,1,indice))+1.)  
          ENDDO  
   
 c   other temperatures  
          DO jk=1,nsoilmx-1  
             DO ig=1,knon  
                ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)  
      $            *ptsoil(ig,jk)  
             ENDDO  
          ENDDO  
   
       ENDIF !--not firstcall  
 c-----------------------------------------------------------------------  
 c   Computation of the Cgrd and Dgrd coefficient for the next step:  
 c   ---------------------------------------------------------------  
   
 c$$$  PB ajout pour cas glace de mer  
       IF (indice .EQ. is_sic) THEN  
           DO ig = 1 , knon  
             ptsoil(ig,nsoilmx) = RTT - 1.8  
           END DO  
       ENDIF  
   
       DO jk=1,nsoilmx  
          zdz2(jk)=dz2(jk)/ptimestep  
       ENDDO  
   
       DO ig=1,knon  
          z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)  
          zc(ig,nsoilmx-1,indice)=  
      $       zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)  
          zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)  
       ENDDO  
   
       DO jk=nsoilmx-1,2,-1  
          DO ig=1,knon  
             z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)  
      $         *(1.-zd(ig,jk,indice)))  
             zc(ig,jk-1,indice)=  
      s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice))  
      $          *z1(ig,indice)  
             zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)  
          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,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))  
          pcapcal(ig)=ztherm_i(ig)*  
      s   (dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))  
          z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.  
          pcapcal(ig)=pcapcal(ig)/z1(ig,indice)  
          pfluxgrd(ig) = pfluxgrd(ig)  
      s   + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)  
      $       - lambda * zc(ig,1,indice)  
      $       - ptsrf(ig))  
      s   /ptimestep  
       ENDDO  
2    
3        RETURN    IMPLICIT NONE
4        END  
5    contains
6    
7      SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
8           pfluxgrd)
9    
10        ! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19
11    
12        USE dimens_m
13        USE indicesol
14        USE dimphy
15        USE dimsoil
16        USE suphec_m
17    
18        ! =======================================================================
19    
20        ! Auteur:  Frederic Hourdin     30/01/92
21        ! -------
22    
23        ! objet:  computation of : the soil temperature evolution
24        ! ------                   the surfacic heat capacity "Capcal"
25        ! the surface conduction flux pcapcal
26    
27    
28        ! Method: implicit time integration
29        ! -------
30        ! Consecutive ground temperatures are related by:
31        ! T(k+1) = C(k) + D(k)*T(k)  (1)
32        ! the coefficients C and D are computed at the t-dt time-step.
33        ! Routine structure:
34        ! 1)new temperatures are computed  using (1)
35        ! 2)C and D coefficients are computed from the new temperature
36        ! profile for the t+dt time-step
37        ! 3)the coefficients A and B are computed where the diffusive
38        ! fluxes at the t+dt time-step is given by
39        ! Fdiff = A + B Ts(t+dt)
40        ! or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
41        ! with F0 = A + B (Ts(t))
42        ! Capcal = B*dt
43    
44        ! Interface:
45        ! ----------
46    
47        ! Arguments:
48        ! ----------
49        ! ptimestep            physical timestep (s)
50        ! indice               sub-surface index
51        ! snow(klon,nbsrf)     snow
52        ! ptsrf(knon)          surface temperature at time-step t (K)
53        ! ptsoil(klon,nsoilmx) temperature inside the ground (K)
54        ! pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
55        ! pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
56    
57        ! =======================================================================
58        ! declarations:
59        ! -------------
60    
61    
62        ! -----------------------------------------------------------------------
63        ! arguments
64        ! ---------
65    
66        REAL ptimestep
67        INTEGER indice, knon
68        REAL ptsrf(knon), ptsoil(klon, nsoilmx), snow(klon)
69        REAL pcapcal(klon), pfluxgrd(klon)
70    
71        ! -----------------------------------------------------------------------
72        ! local arrays
73        ! ------------
74    
75        INTEGER ig, jk
76        ! $$$      REAL zdz2(nsoilmx),z1(klon)
77        REAL zdz2(nsoilmx), z1(klon, nbsrf)
78        REAL min_period, dalph_soil
79        REAL ztherm_i(klon)
80    
81        ! local saved variables:
82        ! ----------------------
83        REAL dz1(nsoilmx), dz2(nsoilmx)
84        ! $$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
85        REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
86        REAL lambda
87        SAVE dz1, dz2, zc, zd, lambda
88        LOGICAL firstsurf(nbsrf)
89        SAVE firstsurf
90        REAL isol, isno, iice
91        SAVE isol, isno, iice
92    
93        DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
94    
95        DATA isol, isno, iice/2000., 2000., 2000./
96    
97        ! -----------------------------------------------------------------------
98        ! Depthts:
99        ! --------
100    
101        REAL rk, fz1, rk1, rk2
102    
103        pfluxgrd(:) = 0.
104        ! calcul de l'inertie thermique a partir de la variable rnat.
105        ! on initialise a iice meme au-dessus d'un point de mer au cas
106        ! ou le point de mer devienne point de glace au pas suivant
107        ! on corrige si on a un point de terre avec ou sans glace
108    
109        IF (indice==is_sic) THEN
110           DO ig = 1, knon
111              ztherm_i(ig) = iice
112              IF (snow(ig)>0.0) ztherm_i(ig) = isno
113           END DO
114        ELSE IF (indice==is_lic) THEN
115           DO ig = 1, knon
116              ztherm_i(ig) = iice
117              IF (snow(ig)>0.0) ztherm_i(ig) = isno
118           END DO
119        ELSE IF (indice==is_ter) THEN
120           DO ig = 1, knon
121              ztherm_i(ig) = isol
122              IF (snow(ig)>0.0) ztherm_i(ig) = isno
123           END DO
124        ELSE IF (indice==is_oce) THEN
125           DO ig = 1, knon
126              ztherm_i(ig) = iice
127           END DO
128        ELSE
129           PRINT *, 'valeur d indice non prevue', indice
130           STOP 1
131        END IF
132    
133        IF (firstsurf(indice)) THEN
134    
135           ! -----------------------------------------------------------------------
136           ! ground levels
137           ! grnd=z/l where l is the skin depth of the diurnal cycle:
138           ! --------------------------------------------------------
139    
140           min_period = 1800. ! en secondes
141           dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
142    
143           OPEN (99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
144           READ (99, *) min_period
145           READ (99, *) dalph_soil
146           PRINT *, 'Discretization for the soil model'
147           PRINT *, 'First level e-folding depth', min_period, '   dalph', &
148                dalph_soil
149           CLOSE (99)
150    9999   CONTINUE
151    
152           ! la premiere couche represente un dixieme de cycle diurne
153           fz1 = sqrt(min_period/3.14)
154    
155           DO jk = 1, nsoilmx
156              rk1 = jk
157              rk2 = jk - 1
158              dz2(jk) = fz(rk1) - fz(rk2)
159           END DO
160           DO jk = 1, nsoilmx - 1
161              rk1 = jk + .5
162              rk2 = jk - .5
163              dz1(jk) = 1./(fz(rk1)-fz(rk2))
164           END DO
165           lambda = fz(.5)*dz1(1)
166           PRINT *, 'full layers, intermediate layers (seconds)'
167           DO jk = 1, nsoilmx
168              rk = jk
169              rk1 = jk + .5
170              rk2 = jk - .5
171              PRINT *, 'fz=', fz(rk1)*fz(rk2)*3.14, fz(rk)*fz(rk)*3.14
172           END DO
173           ! PB
174           firstsurf(indice) = .FALSE.
175    
176           ! Initialisations:
177           ! ----------------
178    
179        ELSE
180           ! -----------------------------------------------------------------------
181           ! Computation of the soil temperatures using the Cgrd and Dgrd
182           ! coefficient computed at the previous time-step:
183           ! -----------------------------------------------
184    
185           ! surface temperature
186           DO ig = 1, knon
187              ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &
188                   indice))+1.)
189           END DO
190    
191           ! other temperatures
192           DO jk = 1, nsoilmx - 1
193              DO ig = 1, knon
194                 ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &
195                      jk)
196              END DO
197           END DO
198    
199        END IF
200        ! -----------------------------------------------------------------------
201        ! Computation of the Cgrd and Dgrd coefficient for the next step:
202        ! ---------------------------------------------------------------
203    
204        ! $$$  PB ajout pour cas glace de mer
205        IF (indice==is_sic) THEN
206           DO ig = 1, knon
207              ptsoil(ig, nsoilmx) = rtt - 1.8
208           END DO
209        END IF
210    
211        DO jk = 1, nsoilmx
212           zdz2(jk) = dz2(jk)/ptimestep
213        END DO
214    
215        DO ig = 1, knon
216           z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
217           zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &
218                z1(ig, indice)
219           zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)
220        END DO
221    
222        DO jk = nsoilmx - 1, 2, -1
223           DO ig = 1, knon
224              z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
225              zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &
226                   )*z1(ig, indice)
227              zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
228           END DO
229        END DO
230    
231        ! -----------------------------------------------------------------------
232        ! computation of the surface diffusive flux from ground and
233        ! calorific capacity of the ground:
234        ! ---------------------------------
235    
236        DO ig = 1, knon
237           pfluxgrd(ig) = ztherm_i(ig)*dz1(1)*(zc(ig,1,indice)+(zd(ig,1, &
238                indice)-1.)*ptsoil(ig,1))
239           pcapcal(ig) = ztherm_i(ig)*(dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
240           z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1.
241           pcapcal(ig) = pcapcal(ig)/z1(ig, indice)
242           pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- &
243                lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep
244        END DO
245    
246      contains
247    
248        real function fz(rk)
249          real rk
250          fz = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
251        end function fz
252    
253      END SUBROUTINE soil
254    
255    end module soil_m
Legend:
Removed from v.38  
changed lines
  Added in v.157
  ViewVC Help
Powered by ViewVC 1.1.21