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