1 | ! |
---|
2 | ! $Header$ |
---|
3 | ! |
---|
4 | SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, & |
---|
5 | lon, lat, ptsoil, pcapcal, pfluxgrd) |
---|
6 | |
---|
7 | USE yomcst_mod_h |
---|
8 | USE dimphy |
---|
9 | USE mod_phys_lmdz_para |
---|
10 | USE indice_sol_mod |
---|
11 | USE print_control_mod, ONLY: lunout |
---|
12 | USE dimsoil_mod_h, ONLY: nsoilmx |
---|
13 | |
---|
14 | IMPLICIT NONE |
---|
15 | |
---|
16 | !======================================================================= |
---|
17 | ! |
---|
18 | ! Auteur: Frederic Hourdin 30/01/92 |
---|
19 | ! ------- |
---|
20 | ! |
---|
21 | ! Object: Computation of : the soil temperature evolution |
---|
22 | ! ------- the surfacic heat capacity "Capcal" |
---|
23 | ! the surface conduction flux pcapcal |
---|
24 | ! |
---|
25 | ! Update: 2021/07 : soil thermal inertia, formerly a constant value, |
---|
26 | ! ------ can also be now a function of soil moisture (F Cheruy's idea) |
---|
27 | ! depending on iflag_inertie, read from physiq.def via conf_phys_m.F90 |
---|
28 | ! ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy) |
---|
29 | ! |
---|
30 | ! Method: Implicit time integration |
---|
31 | ! ------- |
---|
32 | ! Consecutive ground temperatures are related by: |
---|
33 | ! T(k+1) = C(k) + D(k)*T(k) (*) |
---|
34 | ! The coefficients C and D are computed at the t-dt time-step. |
---|
35 | ! Routine structure: |
---|
36 | ! 1) C and D coefficients are computed from the old temperature |
---|
37 | ! 2) new temperatures are computed using (*) |
---|
38 | ! 3) C and D coefficients are computed from the new temperature |
---|
39 | ! profile for the t+dt time-step |
---|
40 | ! 4) the coefficients A and B are computed where the diffusive |
---|
41 | ! fluxes at the t+dt time-step is given by |
---|
42 | ! Fdiff = A + B Ts(t+dt) |
---|
43 | ! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt |
---|
44 | ! with F0 = A + B (Ts(t)) |
---|
45 | ! Capcal = B*dt |
---|
46 | ! |
---|
47 | ! Interface: |
---|
48 | ! ---------- |
---|
49 | ! |
---|
50 | ! Arguments: |
---|
51 | ! ---------- |
---|
52 | ! ptimestep physical timestep (s) |
---|
53 | ! indice sub-surface index |
---|
54 | ! snow(klon) snow |
---|
55 | ! ptsrf(klon) surface temperature at time-step t (K) |
---|
56 | ! qsol(klon) soil moisture (kg/m2 or mm) |
---|
57 | ! lon(klon) longitude in radian |
---|
58 | ! lat(klon) latitude in radian |
---|
59 | ! ptsoil(klon,nsoilmx) temperature inside the ground (K) |
---|
60 | ! pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) |
---|
61 | ! pfluxgrd(klon) surface diffusive flux from ground (Wm-2) |
---|
62 | ! |
---|
63 | !======================================================================= |
---|
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 | |
---|
354 | END SUBROUTINE soil |
---|