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