/[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

revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC revision 101 by guez, Mon Jul 7 17:45:21 2014 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(klon)          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(klon), 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 firstcall, firstsurf(nbsrf)
89        SAVE firstcall, firstsurf
90        REAL isol, isno, iice
91        SAVE isol, isno, iice
92    
93      DO jk = 1, nsoilmx      DATA firstcall/.TRUE./
94        rk1 = jk      DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
       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.  
95    
96      ! Initialisations:      DATA isol, isno, iice/2000., 2000., 2000./
     ! ----------------  
97    
   ELSE !--not firstcall  
98      ! -----------------------------------------------------------------------      ! -----------------------------------------------------------------------
99      ! Computation of the soil temperatures using the Cgrd and Dgrd      ! Depthts:
100      ! coefficient computed at the previous time-step:      ! --------
     ! -----------------------------------------------  
101    
102      ! surface temperature      REAL fz, rk, fz1, rk1, rk2
     DO ig = 1, knon  
       ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &  
         indice))+1.)  
     END DO  
103    
104      ! other temperatures      fz(rk) = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
105      DO jk = 1, nsoilmx - 1      pfluxgrd(:) = 0.
106        DO ig = 1, knon      ! calcul de l'inertie thermique a partir de la variable rnat.
107          ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &      ! on initialise a iice meme au-dessus d'un point de mer au cas
108            jk)      ! ou le point de mer devienne point de glace au pas suivant
109        END DO      ! on corrige si on a un point de terre avec ou sans glace
110      END DO  
111        IF (indice==is_sic) THEN
112           DO ig = 1, knon
113              ztherm_i(ig) = iice
114              IF (snow(ig)>0.0) ztherm_i(ig) = isno
115           END DO
116        ELSE IF (indice==is_lic) THEN
117           DO ig = 1, knon
118              ztherm_i(ig) = iice
119              IF (snow(ig)>0.0) ztherm_i(ig) = isno
120           END DO
121        ELSE IF (indice==is_ter) THEN
122           DO ig = 1, knon
123              ztherm_i(ig) = isol
124              IF (snow(ig)>0.0) ztherm_i(ig) = isno
125           END DO
126        ELSE IF (indice==is_oce) THEN
127           DO ig = 1, knon
128              ztherm_i(ig) = iice
129           END DO
130        ELSE
131           PRINT *, 'valeur d indice non prevue', indice
132           STOP 1
133        END IF
134    
135    
136        ! $$$      IF (firstcall) THEN
137        IF (firstsurf(indice)) THEN
138    
139           ! -----------------------------------------------------------------------
140           ! ground levels
141           ! grnd=z/l where l is the skin depth of the diurnal cycle:
142           ! --------------------------------------------------------
143    
144           min_period = 1800. ! en secondes
145           dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
146    
147           OPEN (99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
148           READ (99, *) min_period
149           READ (99, *) dalph_soil
150           PRINT *, 'Discretization for the soil model'
151           PRINT *, 'First level e-folding depth', min_period, '   dalph', &
152                dalph_soil
153           CLOSE (99)
154    9999   CONTINUE
155    
156           ! la premiere couche represente un dixieme de cycle diurne
157           fz1 = sqrt(min_period/3.14)
158    
159           DO jk = 1, nsoilmx
160              rk1 = jk
161              rk2 = jk - 1
162              dz2(jk) = fz(rk1) - fz(rk2)
163           END DO
164           DO jk = 1, nsoilmx - 1
165              rk1 = jk + .5
166              rk2 = jk - .5
167              dz1(jk) = 1./(fz(rk1)-fz(rk2))
168           END DO
169           lambda = fz(.5)*dz1(1)
170           PRINT *, 'full layers, intermediate layers (seconds)'
171           DO jk = 1, nsoilmx
172              rk = jk
173              rk1 = jk + .5
174              rk2 = jk - .5
175              PRINT *, 'fz=', fz(rk1)*fz(rk2)*3.14, fz(rk)*fz(rk)*3.14
176           END DO
177           ! PB
178           firstsurf(indice) = .FALSE.
179           ! $$$         firstcall =.false.
180    
181           ! Initialisations:
182           ! ----------------
183    
184        ELSE !--not firstcall
185           ! -----------------------------------------------------------------------
186           ! Computation of the soil temperatures using the Cgrd and Dgrd
187           ! coefficient computed at the previous time-step:
188           ! -----------------------------------------------
189    
190           ! surface temperature
191           DO ig = 1, knon
192              ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &
193                   indice))+1.)
194           END DO
195    
196           ! other temperatures
197           DO jk = 1, nsoilmx - 1
198              DO ig = 1, knon
199                 ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &
200                      jk)
201              END DO
202           END DO
203    
204        END IF !--not firstcall
205        ! -----------------------------------------------------------------------
206        ! Computation of the Cgrd and Dgrd coefficient for the next step:
207        ! ---------------------------------------------------------------
208    
209    END IF !--not firstcall      ! $$$  PB ajout pour cas glace de mer
210    ! -----------------------------------------------------------------------      IF (indice==is_sic) THEN
211    ! Computation of the Cgrd and Dgrd coefficient for the next step:         DO ig = 1, knon
212    ! ---------------------------------------------------------------            ptsoil(ig, nsoilmx) = rtt - 1.8
213           END DO
214        END IF
215    
216        DO jk = 1, nsoilmx
217           zdz2(jk) = dz2(jk)/ptimestep
218        END DO
219    
   ! $$$  PB ajout pour cas glace de mer  
   IF (indice==is_sic) THEN  
220      DO ig = 1, knon      DO ig = 1, knon
221        ptsoil(ig, nsoilmx) = rtt - 1.8         z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
222           zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &
223                z1(ig, indice)
224           zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)
225      END DO      END DO
   END IF  
226    
227    DO jk = 1, nsoilmx      DO jk = nsoilmx - 1, 2, -1
228      zdz2(jk) = dz2(jk)/ptimestep         DO ig = 1, knon
229    END DO            z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
230              zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &
231    DO ig = 1, knon                 )*z1(ig, indice)
232      z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)            zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
233      zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &         END DO
234        z1(ig, indice)      END DO
235      zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)  
236    END DO      ! -----------------------------------------------------------------------
237        ! computation of the surface diffusive flux from ground and
238        ! calorific capacity of the ground:
239        ! ---------------------------------
240    
   DO jk = nsoilmx - 1, 2, -1  
241      DO ig = 1, knon      DO ig = 1, knon
242        z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))         pfluxgrd(ig) = ztherm_i(ig)*dz1(1)*(zc(ig,1,indice)+(zd(ig,1, &
243        zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &              indice)-1.)*ptsoil(ig,1))
244          )*z1(ig, indice)         pcapcal(ig) = ztherm_i(ig)*(dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
245        zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)         z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1.
246           pcapcal(ig) = pcapcal(ig)/z1(ig, indice)
247           pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- &
248                lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep
249      END DO      END DO
   END DO  
250    
251    ! -----------------------------------------------------------------------    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  
252    
253    RETURN  end module soil_m
 END SUBROUTINE soil  
Legend:
Removed from v.82  
changed lines
  Added in v.101
  ViewVC Help
Powered by ViewVC 1.1.21