/[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.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC trunk/phylmd/Interface_surf/soil.f revision 254 by guez, Mon Feb 5 10:39:38 2018 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, 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, January 30th, 1992
     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  
   
   
   ! $$$      IF (firstcall) THEN  
   IF (firstsurf(indice)) THEN  
   
     ! -----------------------------------------------------------------------  
     ! ground levels  
     ! 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  
12    
13      ! la premiere couche represente un dixieme de cycle diurne      ! Object: computation of the soil temperature evolution, the heat
14      fz1 = sqrt(min_period/3.14)      ! capacity per unit surface 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        ! Structure of the procedure:
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
29        ! Fdiff = F0 + Soilcap (Ts(t + dt) - Ts(t)) / dt
30        ! with
31        ! F0 = A + B (Ts(t))
32        ! Soilcap = B * dt
33    
34        USE indicesol, only: nbsrf, is_lic, is_oce, is_sic, is_ter
35        USE dimphy, only: klon
36        USE dimsoil, only: nsoilmx
37        USE suphec_m, only: rtt
38    
39        REAL, intent(in):: dtime ! physical timestep (s)
40        INTEGER, intent(in):: nisurf ! sub-surface index
41        REAL, intent(in):: snow(:) ! (knon)
42        REAL, intent(in):: tsurf(:) ! (knon) surface temperature at time-step t (K)
43    
44        real, intent(inout):: tsoil(:, :) ! (knon, nsoilmx)
45        ! temperature inside the ground (K)
46    
47        REAL, intent(out):: soilcap(:) ! (knon)
48        ! specific heat per unit surface (W m-2 s K-1)
49    
50        REAL, intent(out):: soilflux(:) ! (knon)
51        ! surface diffusive flux from ground (W m-2)
52    
53        ! Local:
54    
55        INTEGER knon, ig, jk
56        REAL zdz2(nsoilmx)
57        real z1(size(tsurf), nbsrf) ! (knon, nbsrf)
58        REAL min_period, dalph_soil
59        REAL ztherm_i(size(tsurf)) ! (knon)
60        REAL, save:: dz1(nsoilmx), dz2(nsoilmx)
61        REAL, save:: zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
62        REAL, save:: lambda
63        LOGICAL:: firstsurf(nbsrf) = .TRUE.
64        REAL:: isol = 2000., isno = 2000., iice = 2000.
65    
66        ! Depths:
67        REAL rk, fz1, rk1, rk2
68    
69        !-----------------------------------------------------------------------
70    
71        knon = size(tsurf)
72    
73        ! Calcul de l'inertie thermique. On initialise \`a iice m\^eme
74        ! au-dessus d'un point de mer au cas o\`u le point de mer devienne
75        ! point de glace au pas suivant. On corrige si on a un point de
76        ! terre avec ou sans glace.
77    
78        IF (nisurf==is_sic) THEN
79           DO ig = 1, knon
80              ztherm_i(ig) = iice
81              IF (snow(ig) > 0.0) ztherm_i(ig) = isno
82           END DO
83        ELSE IF (nisurf==is_lic) THEN
84           DO ig = 1, knon
85              ztherm_i(ig) = iice
86              IF (snow(ig) > 0.0) ztherm_i(ig) = isno
87           END DO
88        ELSE IF (nisurf==is_ter) THEN
89           DO ig = 1, knon
90              ztherm_i(ig) = isol
91              IF (snow(ig) > 0.0) ztherm_i(ig) = isno
92           END DO
93        ELSE IF (nisurf==is_oce) THEN
94           DO ig = 1, knon
95              ztherm_i(ig) = iice
96           END DO
97        ELSE
98           PRINT *, 'valeur d indice non prevue', nisurf
99           STOP 1
100        END IF
101    
102        IF (firstsurf(nisurf)) THEN
103           ! ground levels
104           ! grnd=z / l where l is the skin depth of the diurnal cycle:
105    
106           min_period = 1800. ! en secondes
107           dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
108    
109           OPEN(99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
110           READ(99, *) min_period
111           READ(99, *) dalph_soil
112           PRINT *, 'Discretization for the soil model'
113           PRINT *, 'First level e-folding depth', min_period, ' dalph', &
114                dalph_soil
115           CLOSE(99)
116    9999   CONTINUE
117    
118           ! la premiere couche represente un dixieme de cycle diurne
119           fz1 = sqrt(min_period / 3.14)
120    
121           DO jk = 1, nsoilmx
122              rk1 = jk
123              rk2 = jk - 1
124              dz2(jk) = fz(rk1) - fz(rk2)
125           END DO
126           DO jk = 1, nsoilmx - 1
127              rk1 = jk + .5
128              rk2 = jk - .5
129              dz1(jk) = 1. / (fz(rk1) - fz(rk2))
130           END DO
131           lambda = fz(.5) * dz1(1)
132           PRINT *, 'full layers, intermediate layers (seconds)'
133           DO jk = 1, nsoilmx
134              rk = jk
135              rk1 = jk + .5
136              rk2 = jk - .5
137              PRINT *, 'fz=', fz(rk1) * fz(rk2) * 3.14, fz(rk) * fz(rk) * 3.14
138           END DO
139           ! PB
140           firstsurf(nisurf) = .FALSE.
141        ELSE
142           ! Computation of the soil temperatures using the Zc and Zd
143           ! coefficient computed at the previous time-step:
144    
145           ! surface temperature
146           DO ig = 1, knon
147              tsoil(ig, 1) = (lambda * zc(ig, 1, nisurf) + tsurf(ig)) &
148                   / (lambda * (1. - zd(ig, 1, nisurf)) + 1.)
149           END DO
150    
151           ! other temperatures
152           DO jk = 1, nsoilmx - 1
153              DO ig = 1, knon
154                 tsoil(ig, jk + 1) = zc(ig, jk, nisurf) &
155                      + zd(ig, jk, nisurf) * tsoil(ig, jk)
156              END DO
157           END DO
158        END IF
159    
160        ! Computation of the Zc and Zd coefficient for the next step:
161    
162        IF (nisurf==is_sic) THEN
163           DO ig = 1, knon
164              tsoil(ig, nsoilmx) = rtt - 1.8
165           END DO
166        END IF
167    
168      DO jk = 1, nsoilmx      DO jk = 1, nsoilmx
169        rk1 = jk         zdz2(jk) = dz2(jk) / dtime
       rk2 = jk - 1  
       dz2(jk) = fz(rk1) - fz(rk2)  
170      END DO      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.  
   
     ! Initialisations:  
     ! ----------------  
   
   ELSE !--not firstcall  
     ! -----------------------------------------------------------------------  
     ! Computation of the soil temperatures using the Cgrd and Dgrd  
     ! coefficient computed at the previous time-step:  
     ! -----------------------------------------------  
171    
     ! surface temperature  
172      DO ig = 1, knon      DO ig = 1, knon
173        ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &         z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx - 1)
174          indice))+1.)         zc(ig, nsoilmx - 1, nisurf) = zdz2(nsoilmx) * tsoil(ig, nsoilmx) / &
175                z1(ig, nisurf)
176           zd(ig, nsoilmx - 1, nisurf) = dz1(nsoilmx - 1) / z1(ig, nisurf)
177      END DO      END DO
178    
179      ! other temperatures      DO jk = nsoilmx - 1, 2, - 1
180      DO jk = 1, nsoilmx - 1         DO ig = 1, knon
181        DO ig = 1, knon            z1(ig, nisurf) = 1. / (zdz2(jk) + dz1(jk - 1) &
182          ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &                 + dz1(jk) * (1. - zd(ig, jk, nisurf)))
183            jk)            zc(ig, jk - 1, nisurf) = (tsoil(ig, jk) * zdz2(jk) &
184        END DO                 + dz1(jk) * zc(ig, jk, nisurf)) * z1(ig, nisurf)
185              zd(ig, jk - 1, nisurf) = dz1(jk - 1) * z1(ig, nisurf)
186           END DO
187      END DO      END DO
188    
189    END IF !--not firstcall      ! computation of the surface diffusive flux from ground and
190    ! -----------------------------------------------------------------------      ! calorific capacity of the ground:
   ! Computation of the Cgrd and Dgrd coefficient for the next step:  
   ! ---------------------------------------------------------------  
191    
   ! $$$  PB ajout pour cas glace de mer  
   IF (indice==is_sic) THEN  
192      DO ig = 1, knon      DO ig = 1, knon
193        ptsoil(ig, nsoilmx) = rtt - 1.8         soilflux(ig) = ztherm_i(ig) * dz1(1) * (zc(ig, 1, nisurf) + (zd(ig, 1, &
194                nisurf) - 1.) * tsoil(ig, 1))
195           soilcap(ig) = ztherm_i(ig) * (dz2(1) &
196                + dtime * (1. - zd(ig, 1, nisurf)) * dz1(1))
197           z1(ig, nisurf) = lambda * (1. - zd(ig, 1, nisurf)) + 1.
198           soilcap(ig) = soilcap(ig) / z1(ig, nisurf)
199           soilflux(ig) = soilflux(ig) + soilcap(ig) * (tsoil(ig, 1) &
200                * z1(ig, nisurf) - lambda * zc(ig, 1, nisurf) - tsurf(ig)) / dtime
201      END DO      END DO
   END IF  
202    
203    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  
204    
205    DO jk = nsoilmx - 1, 2, -1      pure real function fz(rk)
206      DO ig = 1, knon  
207        z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))        real, intent(in):: rk
208        zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &  
209          )*z1(ig, indice)        !-----------------------------------------
210        zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)  
211      END DO        fz = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
212    END DO  
213        end function fz
214    
215    ! -----------------------------------------------------------------------    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  
216    
217    RETURN  end module soil_m
 END SUBROUTINE soil  
Legend:
Removed from v.82  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21