source: LMDZ6/branches/Amaury_dev/libf/phylmd/soil.F90 @ 5134

Last change on this file since 5134 was 5112, checked in by abarral, 4 months ago

Rename modules in phy_common from *_mod > lmdz_*
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.0 KB
Line 
1! $Header$
2
3SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, &
4        lon, lat, ptsoil, pcapcal, pfluxgrd)
5
6  USE dimphy
7  USE lmdz_phys_para
8  USE indice_sol_mod
9  USE lmdz_print_control, ONLY: lunout
10  USE lmdz_abort_physic, ONLY: abort_physic
11
12  IMPLICIT NONE
13
14  !=======================================================================
15
16  !   Auteur:  Frederic Hourdin     30/01/92
17  !   -------
18
19  !   Object:  Computation of : the soil temperature evolution
20  !   -------                   the surfacic heat capacity "Capcal"
21  !                            the surface conduction flux pcapcal
22
23  !   Update: 2021/07 : soil thermal inertia, formerly a constant value,
24  !   ------   can also be now a function of soil moisture (F Cheruy's idea)
25  !            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
26  !            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
27
28  !   Method: Implicit time integration
29  !   -------
30  !   Consecutive ground temperatures are related by:
31  !           T(k+1) = C(k) + D(k)*T(k)  (*)
32  !   The coefficients C and D are computed at the t-dt time-step.
33  !   Routine structure:
34  !   1) C and D coefficients are computed from the old temperature
35  !   2) new temperatures are computed using (*)
36  !   3) C and D coefficients are computed from the new temperature
37  !      profile for the t+dt time-step
38  !   4) the coefficients A and B are computed where the diffusive
39  !      fluxes at the t+dt time-step is given by
40  !             Fdiff = A + B Ts(t+dt)
41  !      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
42  !             with F0 = A + B (Ts(t))
43  !                 Capcal = B*dt
44
45  !   Interface:
46  !   ----------
47
48  !   Arguments:
49  !   ----------
50  !   ptimestep            physical timestep (s)
51  !   indice               sub-surface index
52  !   snow(klon)           snow
53  !   ptsrf(klon)          surface temperature at time-step t (K)
54  !   qsol(klon)           soil moisture (kg/m2 or mm)
55  !   lon(klon)            longitude in radian
56  !   lat(klon)            latitude in radian
57  !   ptsoil(klon,nsoilmx) temperature inside the ground (K)
58  !   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
59  !   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
60
61  !=======================================================================
62  INCLUDE "YOMCST.h"
63  INCLUDE "dimsoil.h"
64  INCLUDE "comsoil.h"
65  !-----------------------------------------------------------------------
66  ! Arguments
67  ! ---------
68  REAL, INTENT(IN) :: ptimestep
69  INTEGER, INTENT(IN) :: indice, knon !, knindex
70  REAL, DIMENSION(klon), INTENT(IN) :: snow
71  REAL, DIMENSION(klon), INTENT(IN) :: ptsrf
72  REAL, DIMENSION(klon), INTENT(IN) :: qsol
73  REAL, DIMENSION(klon), INTENT(IN) :: lon
74  REAL, DIMENSION(klon), INTENT(IN) :: lat
75
76  REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: ptsoil
77  REAL, DIMENSION(klon), INTENT(OUT) :: pcapcal
78  REAL, DIMENSION(klon), INTENT(OUT) :: pfluxgrd
79
80  !-----------------------------------------------------------------------
81  ! Local variables
82  ! ---------------
83  INTEGER :: ig, jk, ierr
84  REAL :: min_period, dalph_soil
85  REAL, DIMENSION(nsoilmx) :: zdz2
86  REAL :: z1s
87  REAL, DIMENSION(klon) :: ztherm_i
88  REAL, DIMENSION(klon, nsoilmx, nbsrf) :: C_coef, D_coef
89
90  ! Local saved variables
91  ! ---------------------
92  REAL, SAVE :: lambda
93  !$OMP THREADPRIVATE(lambda)
94  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
95  !$OMP THREADPRIVATE(dz1,dz2)
96  LOGICAL, SAVE :: firstcall = .TRUE.
97  !$OMP THREADPRIVATE(firstcall)
98
99  !-----------------------------------------------------------------------
100  !   Depthts:
101  !   --------
102  REAL fz, rk, fz1, rk1, rk2
103  fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
104
105
106  !-----------------------------------------------------------------------
107  ! Calculation of some constants
108  ! NB! These constants do not depend on the sub-surfaces
109  !-----------------------------------------------------------------------
110
111  IF (firstcall) THEN
112    !-----------------------------------------------------------------------
113    !   ground levels
114    !   grnd=z/l where l is the skin depth of the diurnal cycle:
115    !-----------------------------------------------------------------------
116
117    min_period = 1800. ! en secondes
118    dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
119    !$OMP MASTER
120    IF (is_mpi_root) THEN
121      OPEN(99, file = 'soil.def', status = 'old', form = 'formatted', iostat = ierr)
122      IF (ierr == 0) THEN ! Read file only if it exists
123        READ(99, *) min_period
124        READ(99, *) dalph_soil
125        WRITE(lunout, *)'Discretization for the soil model'
126        WRITE(lunout, *)'First level e-folding depth', min_period, &
127                '   dalph', dalph_soil
128        CLOSE(99)
129      END IF
130    ENDIF
131    !$OMP END MASTER
132    CALL bcast(min_period)
133    CALL bcast(dalph_soil)
134
135    !   la premiere couche represente un dixieme de cycle diurne
136    fz1 = SQRT(min_period / 3.14)
137
138    DO jk = 1, nsoilmx
139      rk1 = jk
140      rk2 = jk - 1
141      dz2(jk) = fz(rk1) - fz(rk2)
142    ENDDO
143    DO jk = 1, nsoilmx - 1
144      rk1 = jk + .5
145      rk2 = jk - .5
146      dz1(jk) = 1. / (fz(rk1) - fz(rk2))
147    ENDDO
148    lambda = fz(.5) * dz1(1)
149    WRITE(lunout, *)'full layers, intermediate layers (seconds)'
150    DO jk = 1, nsoilmx
151      rk = jk
152      rk1 = jk + .5
153      rk2 = jk - .5
154      WRITE(lunout, *)'fz=', &
155              fz(rk1) * fz(rk2) * 3.14, fz(rk) * fz(rk) * 3.14
156    ENDDO
157
158    firstcall = .FALSE.
159  END IF
160
161
162  !-----------------------------------------------------------------------
163  !   Calcul de l'inertie thermique a partir de la variable rnat.
164  !   on initialise a inertie_sic meme au-dessus d'un point de mer au cas
165  !   ou le point de mer devienne point de glace au pas suivant
166  !   on corrige si on a un point de terre avec ou sans glace
167
168  !   iophys can be used to write the ztherm_i variable in a phys.nc file
169  !   and check the results; to do so, add "CALL iophys_ini" in physiq_mod
170  !   and add knindex to the list of inputs in all the calls to soil.F90
171  !   (and to soil.F90 itself !)
172  !-----------------------------------------------------------------------
173
174  IF (indice == is_sic) THEN
175    DO ig = 1, knon
176      ztherm_i(ig) = inertie_sic
177    ENDDO
178    IF (iflag_sic == 0) THEN
179      DO ig = 1, knon
180        IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno
181      ENDDO
182      !      Otherwise sea-ice keeps the same inertia, even when covered by snow
183    ENDIF
184    !    CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', &
185    !      knon, knindex, ztherm_i)
186  ELSE IF (indice == is_lic) THEN
187    DO ig = 1, knon
188      ztherm_i(ig) = inertie_lic
189      IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno
190    ENDDO
191    !    CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', &
192    !      knon, knindex, ztherm_i)
193  ELSE IF (indice == is_ter) THEN
194
195    ! La relation entre l'inertie thermique du sol et qsol change d'apres
196    !   iflag_inertie, defini dans physiq.def, et appele via comsoil.h
197
198    DO ig = 1, knon
199      ! iflag_inertie=0 correspond au cas inertie=constant, comme avant
200      IF (iflag_inertie==0) THEN
201        ztherm_i(ig) = inertie_sol
202      ELSE IF (iflag_inertie == 1) THEN
203        ! I = a_qsol * qsol + b  modele lineaire deduit d'une
204        ! regression lineaire I = a_mrsos * mrsos + b obtenue sur
205        ! sorties MO d'une simulation LMDZOR(CMIP6) sur l'annee 2000
206        ! sur tous les points avec frac_snow=0
207        ! Difference entre qsol et mrsos prise en compte par un
208        ! facteur d'echelle sur le coefficient directeur de regression:
209        ! fact = 35./150. = mrsos_max/qsol_max
210        ! et a_qsol = a_mrsos * fact (car a = dI/dHumidite)
211        ztherm_i(ig) = 30.0 * 35.0 / 150.0 * qsol(ig) + 770.0
212        ! AS : pour qsol entre 0 - 150, on a I entre 770 - 1820
213      ELSE IF (iflag_inertie == 2) THEN
214        ! deux regressions lineaires, sur les memes sorties,
215        ! distinguant le type de sol : sable ou autre (limons/argile)
216        ! Implementation simple : regression type "sable" seulement pour
217        ! Sahara, defini par une "boite" lat/lon (NB : en radians !! )
218        IF (lon(ig)>-0.35 .AND. lon(ig)<0.70 .AND. lat(ig)>0.17 .AND. lat(ig)<0.52) THEN
219          ! Valeurs theoriquement entre 728 et 2373 ; qsol valeurs basses
220          ztherm_i(ig) = 47. * 35.0 / 150.0 * qsol(ig) + 728.  ! boite type "sable" pour Sahara
221        ELSE
222          ! Valeurs theoriquement entre 550 et 1940 ; qsol valeurs moyennes et hautes
223          ztherm_i(ig) = 41. * 35.0 / 150.0 * qsol(ig) + 505.
224        ENDIF
225      ELSE IF (iflag_inertie == 3) THEN
226        ! AS : idee a tester :
227        ! si la relation doit etre une droite,
228        ! definissons-la en fonction des valeurs min et max de qsol (0:150),
229        ! et de l'inertie (900 : 2000 ou 2400 ; choix ici: 2000)
230        ! I = I_min + qsol * (I_max - I_min)/(qsol_max - qsol_min)
231        ztherm_i(ig) = 900. + qsol(ig) * (2000. - 900.) / 150.
232      ELSE
233        WRITE (lunout, *) "Le choix iflag_inertie = ", iflag_inertie, " n'est pas defini. Veuillez choisir un entier entre 0 et 3"
234      ENDIF
235
236      ! Fin de l'introduction de la relation entre l'inertie thermique du sol et qsol
237      !-------------------------------------------
238      !AS : donc le moindre flocon de neige sur un point de grid
239      ! fait que l'inertie du point passe a la valeur pour neige !
240      IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno
241
242    ENDDO
243    !    CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', &
244    !      knon, knindex, ztherm_i)
245  ELSE IF (indice == is_oce) THEN
246    DO ig = 1, knon
247      !       This is just in case, but SST should be used by the model anyway
248      ztherm_i(ig) = inertie_sic
249    ENDDO
250    !    CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', &
251    !      knon, knindex, ztherm_i)
252  ELSE
253    WRITE(lunout, *) "valeur d indice non prevue", indice
254    CALL abort_physic("soil", "", 1)
255  ENDIF
256
257
258  !-----------------------------------------------------------------------
259  ! 1)
260  ! Calculation of Cgrf and Dgrd coefficients using soil temperature from
261  ! previous time step.
262
263  ! These variables are recalculated on the local compressed grid instead
264  ! of saved in restart file.
265  !-----------------------------------------------------------------------
266  DO jk = 1, nsoilmx
267    zdz2(jk) = dz2(jk) / ptimestep
268  ENDDO
269
270  DO ig = 1, knon
271    z1s = zdz2(nsoilmx) + dz1(nsoilmx - 1)
272    C_coef(ig, nsoilmx - 1, indice) = &
273            zdz2(nsoilmx) * ptsoil(ig, nsoilmx) / z1s
274    D_coef(ig, nsoilmx - 1, indice) = dz1(nsoilmx - 1) / z1s
275  ENDDO
276
277  DO jk = nsoilmx - 1, 2, -1
278    DO ig = 1, knon
279      z1s = 1. / (zdz2(jk) + dz1(jk - 1) + dz1(jk) &
280              * (1. - D_coef(ig, jk, indice)))
281      C_coef(ig, jk - 1, indice) = &
282              (ptsoil(ig, jk) * zdz2(jk) + dz1(jk) * C_coef(ig, jk, indice)) * z1s
283      D_coef(ig, jk - 1, indice) = dz1(jk - 1) * z1s
284    ENDDO
285  ENDDO
286
287  !-----------------------------------------------------------------------
288  ! 2)
289  ! Computation of the soil temperatures using the Cgrd and Dgrd
290  ! coefficient computed above
291
292  !-----------------------------------------------------------------------
293
294  !    Surface temperature
295  DO ig = 1, knon
296    ptsoil(ig, 1) = (lambda * C_coef(ig, 1, indice) + ptsrf(ig)) / &
297            (lambda * (1. - D_coef(ig, 1, indice)) + 1.)
298  ENDDO
299
300  !   Other temperatures
301  DO jk = 1, nsoilmx - 1
302    DO ig = 1, knon
303      ptsoil(ig, jk + 1) = C_coef(ig, jk, indice) + D_coef(ig, jk, indice) &
304              * ptsoil(ig, jk)
305    ENDDO
306  ENDDO
307
308  IF (indice == is_sic) THEN
309    DO ig = 1, knon
310      ptsoil(ig, nsoilmx) = RTT - 1.8
311    END DO
312  ENDIF
313
314  !-----------------------------------------------------------------------
315  ! 3)
316  ! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
317  ! temperature
318  !-----------------------------------------------------------------------
319  DO ig = 1, knon
320    z1s = zdz2(nsoilmx) + dz1(nsoilmx - 1)
321    C_coef(ig, nsoilmx - 1, indice) = zdz2(nsoilmx) * ptsoil(ig, nsoilmx) / z1s
322    D_coef(ig, nsoilmx - 1, indice) = dz1(nsoilmx - 1) / z1s
323  ENDDO
324
325  DO jk = nsoilmx - 1, 2, -1
326    DO ig = 1, knon
327      z1s = 1. / (zdz2(jk) + dz1(jk - 1) + dz1(jk) &
328              * (1. - D_coef(ig, jk, indice)))
329      C_coef(ig, jk - 1, indice) = &
330              (ptsoil(ig, jk) * zdz2(jk) + dz1(jk) * C_coef(ig, jk, indice)) * z1s
331      D_coef(ig, jk - 1, indice) = dz1(jk - 1) * z1s
332    ENDDO
333  ENDDO
334
335  !-----------------------------------------------------------------------
336  ! 4)
337  ! Computation of the surface diffusive flux from ground and
338  ! calorific capacity of the ground
339  !-----------------------------------------------------------------------
340  DO ig = 1, knon
341    pfluxgrd(ig) = ztherm_i(ig) * dz1(1) * &
342            (C_coef(ig, 1, indice) + (D_coef(ig, 1, indice) - 1.) * ptsoil(ig, 1))
343    pcapcal(ig) = ztherm_i(ig) * &
344            (dz2(1) + ptimestep * (1. - D_coef(ig, 1, indice)) * dz1(1))
345    z1s = lambda * (1. - D_coef(ig, 1, indice)) + 1.
346    pcapcal(ig) = pcapcal(ig) / z1s
347    pfluxgrd(ig) = pfluxgrd(ig) &
348            + pcapcal(ig) * (ptsoil(ig, 1) * z1s &
349                    - lambda * C_coef(ig, 1, indice) &
350                    - ptsrf(ig)) &
351                    / ptimestep
352  ENDDO
353
354END SUBROUTINE soil
Note: See TracBrowser for help on using the repository browser.