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