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