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