/[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 134 by guez, Wed Apr 29 15:47:56 2015 UTC trunk/phylmd/Interface_surf/soil.f revision 254 by guez, Mon Feb 5 10:39:38 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(dtime, 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 indicesol, only: nbsrf, is_lic, is_oce, is_sic, is_ter
35      ! Arguments:      USE dimphy, only: klon
36      ! ----------      USE dimsoil, only: nsoilmx
37      ! ptimestep            physical timestep (s)      USE suphec_m, only: rtt
38      ! indice               sub-surface index  
39      ! snow(klon,nbsrf)     snow      REAL, intent(in):: dtime ! physical timestep (s)
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 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:  
     ! --------  
65    
66        ! Depths:
67      REAL rk, fz1, rk1, rk2      REAL rk, fz1, rk1, rk2
68    
69      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  
70    
71      IF (indice==is_sic) THEN      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         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(nisurf)) THEN
     ! $$$      IF (firstcall) THEN  
     IF (firstsurf(indice)) 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 163  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.
141         ! $$$         firstcall =.false.      ELSE
142           ! Computation of the soil temperatures using the Zc and Zd
        ! Initialisations:  
        ! ----------------  
   
     ELSE !--not firstcall  
        ! -----------------------------------------------------------------------  
        ! 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
159    
160      END IF !--not firstcall      ! Computation of the Zc and Zd coefficient for the next step:
     ! -----------------------------------------------------------------------  
     ! Computation of the Cgrd and Dgrd coefficient for the next step:  
     ! ---------------------------------------------------------------  
161    
162      ! $$$  PB ajout pour cas glace de mer      IF (nisurf==is_sic) THEN
     IF (indice==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) / dtime
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.              + dtime * (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)) / dtime
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.134  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21