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

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

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

trunk/phylmd/soil.f revision 82 by guez, Wed Mar 5 14:57:53 2014 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
2    
 ! $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, pcapcal, &  
     pfluxgrd)  
   USE dimens_m  
   USE indicesol  
   USE dimphy  
   USE dimsoil  
   USE suphec_m  
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    ! =======================================================================  contains
6    
7    ! Auteur:  Frederic Hourdin     30/01/92    SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
8    ! -------         pfluxgrd)
9    
10    ! objet:  computation of : the soil temperature evolution      ! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19
   ! ------                   the surfacic heat capacity "Capcal"  
   ! the surface conduction flux pcapcal  
   
   
   ! Method: implicit time integration  
   ! -------  
   ! Consecutive ground temperatures are related by:  
   ! T(k+1) = C(k) + D(k)*T(k)  (1)  
   ! the coefficients C and D are computed at the t-dt time-step.  
   ! Routine structure:  
   ! 1)new temperatures are computed  using (1)  
   ! 2)C and D coefficients are computed from the new temperature  
   ! profile for the t+dt time-step  
   ! 3)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,nbsrf)     snow  
   ! ptsrf(klon)          surface temperature at time-step t (K)  
   ! 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)  
   
   ! =======================================================================  
   ! declarations:  
   ! -------------  
   
   
   ! -----------------------------------------------------------------------  
   ! arguments  
   ! ---------  
   
   REAL ptimestep  
   INTEGER indice, knon  
   REAL ptsrf(klon), ptsoil(klon, nsoilmx), snow(klon)  
   REAL pcapcal(klon), pfluxgrd(klon)  
   
   ! -----------------------------------------------------------------------  
   ! local arrays  
   ! ------------  
   
   INTEGER ig, jk  
   ! $$$      REAL zdz2(nsoilmx),z1(klon)  
   REAL zdz2(nsoilmx), z1(klon, nbsrf)  
   REAL min_period, dalph_soil  
   REAL ztherm_i(klon)  
   
   ! local saved variables:  
   ! ----------------------  
   REAL dz1(nsoilmx), dz2(nsoilmx)  
   ! $$$          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./  
   
   ! -----------------------------------------------------------------------  
   ! Depthts:  
   ! --------  
   
   REAL fz, rk, fz1, rk1, rk2  
   
   fz(rk) = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)  
   pfluxgrd(:) = 0.  
   ! calcul de l'inertie thermique a partir de la variable rnat.  
   ! on initialise a iice 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  
11    
12    IF (indice==is_sic) THEN      USE dimens_m
13      DO ig = 1, knon      USE indicesol
14        ztherm_i(ig) = iice      USE dimphy
15        IF (snow(ig)>0.0) ztherm_i(ig) = isno      USE dimsoil
16      END DO      USE suphec_m
17    ELSE IF (indice==is_lic) THEN  
18      DO ig = 1, knon      ! =======================================================================
19        ztherm_i(ig) = iice  
20        IF (snow(ig)>0.0) ztherm_i(ig) = isno      ! Auteur:  Frederic Hourdin     30/01/92
21      END DO      ! -------
22    ELSE IF (indice==is_ter) THEN  
23      DO ig = 1, knon      ! objet:  computation of : the soil temperature evolution
24        ztherm_i(ig) = isol      ! ------                   the surfacic heat capacity "Capcal"
25        IF (snow(ig)>0.0) ztherm_i(ig) = isno      ! the surface conduction flux pcapcal
26      END DO  
27    ELSE IF (indice==is_oce) THEN  
28      DO ig = 1, knon      ! Method: implicit time integration
29        ztherm_i(ig) = iice      ! -------
30      END DO      ! Consecutive ground temperatures are related by:
31    ELSE      ! T(k+1) = C(k) + D(k)*T(k)  (1)
32      PRINT *, 'valeur d indice non prevue', indice      ! the coefficients C and D are computed at the t-dt time-step.
33      STOP 1      ! Routine structure:
34    END IF      ! 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    ! $$$      IF (firstcall) THEN      ! -----------------------------------------------------------------------
63    IF (firstsurf(indice)) THEN      ! 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      ! ground levels      ! local arrays
73      ! 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.  
   
     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, '   dalph', &  
       dalph_soil  
     CLOSE (99)  
 9999 CONTINUE  
74    
75      ! la premiere couche represente un dixieme de cycle diurne      INTEGER ig, jk
76      fz1 = sqrt(min_period/3.14)      ! $$$      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      DO jk = 1, nsoilmx      DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
       rk1 = jk  
       rk2 = jk - 1  
       dz2(jk) = fz(rk1) - fz(rk2)  
     END DO  
     DO jk = 1, nsoilmx - 1  
       rk1 = jk + .5  
       rk2 = jk - .5  
       dz1(jk) = 1./(fz(rk1)-fz(rk2))  
     END DO  
     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  
     END DO  
     ! PB  
     firstsurf(indice) = .FALSE.  
     ! $$$         firstcall =.false.  
94    
95      ! Initialisations:      DATA isol, isno, iice/2000., 2000., 2000./
     ! ----------------  
96    
   ELSE !--not firstcall  
97      ! -----------------------------------------------------------------------      ! -----------------------------------------------------------------------
98      ! Computation of the soil temperatures using the Cgrd and Dgrd      ! Depthts:
99      ! coefficient computed at the previous time-step:      ! --------
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    
     ! surface temperature  
215      DO ig = 1, knon      DO ig = 1, knon
216        ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &         z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
217          indice))+1.)         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      END DO
221    
222      ! other temperatures      DO jk = nsoilmx - 1, 2, -1
223      DO jk = 1, nsoilmx - 1         DO ig = 1, knon
224        DO ig = 1, knon            z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
225          ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &            zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &
226            jk)                 )*z1(ig, indice)
227        END DO            zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
228           END DO
229      END DO      END DO
230    
231    END IF !--not firstcall      ! -----------------------------------------------------------------------
232    ! -----------------------------------------------------------------------      ! computation of the surface diffusive flux from ground and
233    ! Computation of the Cgrd and Dgrd coefficient for the next step:      ! calorific capacity of the ground:
234    ! ---------------------------------------------------------------      ! ---------------------------------
235    
   ! $$$  PB ajout pour cas glace de mer  
   IF (indice==is_sic) THEN  
236      DO ig = 1, knon      DO ig = 1, knon
237        ptsoil(ig, nsoilmx) = rtt - 1.8         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      END DO
   END IF  
245    
246    DO jk = 1, nsoilmx    contains
     zdz2(jk) = dz2(jk)/ptimestep  
   END DO  
   
   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)  
   END DO  
247    
248    DO jk = nsoilmx - 1, 2, -1      real function fz(rk)
249      DO ig = 1, knon        real rk
250        z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))        fz = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
251        zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &      end function fz
         )*z1(ig, indice)  
       zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)  
     END DO  
   END DO  
252    
253    ! -----------------------------------------------------------------------    END SUBROUTINE soil
   ! 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)*(zc(ig,1,indice)+(zd(ig,1, &  
       indice)-1.)*ptsoil(ig,1))  
     pcapcal(ig) = ztherm_i(ig)*(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) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- &  
       lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep  
   END DO  
254    
255    RETURN  end module soil_m
 END SUBROUTINE soil  
Legend:
Removed from v.82  
changed lines
  Added in v.157
  ViewVC Help
Powered by ViewVC 1.1.21