1 |
|
|
! |
2 |
|
|
! $Header$ |
3 |
|
|
! |
4 |
|
548379 |
SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, & |
5 |
|
864 |
lon, lat, ptsoil, pcapcal, pfluxgrd) |
6 |
|
|
|
7 |
|
|
USE dimphy |
8 |
|
|
USE mod_phys_lmdz_para |
9 |
|
|
USE indice_sol_mod |
10 |
|
|
USE print_control_mod, ONLY: lunout |
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 |
|
1728 |
REAL, DIMENSION(klon) :: ztherm_i |
88 |
|
1728 |
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 |
✓✓ |
864 |
IF (firstcall) THEN |
112 |
|
|
!----------------------------------------------------------------------- |
113 |
|
|
! ground levels |
114 |
|
|
! grnd=z/l where l is the skin depth of the diurnal cycle: |
115 |
|
|
!----------------------------------------------------------------------- |
116 |
|
|
|
117 |
|
1 |
min_period=1800. ! en secondes |
118 |
|
1 |
dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ. |
119 |
|
|
!$OMP MASTER |
120 |
✓✗ |
1 |
IF (is_mpi_root) THEN |
121 |
|
1 |
OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr) |
122 |
✗✓ |
1 |
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 |
|
1 |
CALL bcast(min_period) |
133 |
|
1 |
CALL bcast(dalph_soil) |
134 |
|
|
|
135 |
|
|
! la premiere couche represente un dixieme de cycle diurne |
136 |
|
1 |
fz1=SQRT(min_period/3.14) |
137 |
|
|
|
138 |
✓✓ |
12 |
DO jk=1,nsoilmx |
139 |
|
11 |
rk1=jk |
140 |
|
11 |
rk2=jk-1 |
141 |
|
12 |
dz2(jk)=fz(rk1)-fz(rk2) |
142 |
|
|
ENDDO |
143 |
✓✓ |
11 |
DO jk=1,nsoilmx-1 |
144 |
|
10 |
rk1=jk+.5 |
145 |
|
10 |
rk2=jk-.5 |
146 |
|
11 |
dz1(jk)=1./(fz(rk1)-fz(rk2)) |
147 |
|
|
ENDDO |
148 |
|
1 |
lambda=fz(.5)*dz1(1) |
149 |
|
1 |
WRITE(lunout,*)'full layers, intermediate layers (seconds)' |
150 |
✓✓ |
12 |
DO jk=1,nsoilmx |
151 |
|
11 |
rk=jk |
152 |
|
11 |
rk1=jk+.5 |
153 |
|
11 |
rk2=jk-.5 |
154 |
|
11 |
WRITE(lunout,*)'fz=', & |
155 |
|
23 |
fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 |
156 |
|
|
ENDDO |
157 |
|
|
|
158 |
|
1 |
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 |
✓✓ |
864 |
IF (indice == is_sic) THEN |
175 |
✓✓ |
62978 |
DO ig = 1, knon |
176 |
|
62978 |
ztherm_i(ig) = inertie_sic |
177 |
|
|
ENDDO |
178 |
✗✓ |
288 |
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 |
✓✓ |
576 |
ELSE IF (indice == is_lic) THEN |
187 |
✓✓ |
44064 |
DO ig = 1, knon |
188 |
|
43776 |
ztherm_i(ig) = inertie_lic |
189 |
✓✓ |
44064 |
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 |
✓✗ |
288 |
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 |
✓✓ |
148896 |
DO ig = 1, knon |
199 |
|
|
! iflag_inertie=0 correspond au cas inertie=constant, comme avant |
200 |
✓✗ |
148608 |
IF (iflag_inertie==0) THEN |
201 |
|
148608 |
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 |
✓✓ |
148896 |
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 |
✓✓ |
10368 |
DO jk=1,nsoilmx |
267 |
|
10368 |
zdz2(jk)=dz2(jk)/ptimestep |
268 |
|
|
ENDDO |
269 |
|
|
|
270 |
✓✓ |
255938 |
DO ig=1,knon |
271 |
|
255074 |
z1s = zdz2(nsoilmx)+dz1(nsoilmx-1) |
272 |
|
|
C_coef(ig,nsoilmx-1,indice)= & |
273 |
|
255074 |
zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s |
274 |
|
255938 |
D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s |
275 |
|
|
ENDDO |
276 |
|
|
|
277 |
✓✓ |
8640 |
DO jk=nsoilmx-1,2,-1 |
278 |
✓✓ |
2304306 |
DO ig=1,knon |
279 |
|
|
z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) & |
280 |
|
2295666 |
*(1.-D_coef(ig,jk,indice))) |
281 |
|
|
C_coef(ig,jk-1,indice)= & |
282 |
|
2295666 |
(ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s |
283 |
|
2303442 |
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 |
✓✓ |
255938 |
DO ig=1,knon |
296 |
|
|
ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/ & |
297 |
|
255938 |
(lambda*(1.-D_coef(ig,1,indice))+1.) |
298 |
|
|
ENDDO |
299 |
|
|
|
300 |
|
|
! Other temperatures |
301 |
✓✓ |
9504 |
DO jk=1,nsoilmx-1 |
302 |
✓✓ |
2560244 |
DO ig=1,knon |
303 |
|
|
ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) & |
304 |
|
2559380 |
*ptsoil(ig,jk) |
305 |
|
|
ENDDO |
306 |
|
|
ENDDO |
307 |
|
|
|
308 |
✓✓ |
864 |
IF (indice == is_sic) THEN |
309 |
✓✓ |
62978 |
DO ig = 1 , knon |
310 |
|
62978 |
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 |
✓✓ |
255938 |
DO ig=1,knon |
320 |
|
255074 |
z1s = zdz2(nsoilmx)+dz1(nsoilmx-1) |
321 |
|
255074 |
C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s |
322 |
|
255938 |
D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s |
323 |
|
|
ENDDO |
324 |
|
|
|
325 |
✓✓ |
8640 |
DO jk=nsoilmx-1,2,-1 |
326 |
✓✓ |
2304306 |
DO ig=1,knon |
327 |
|
|
z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) & |
328 |
|
2295666 |
*(1.-D_coef(ig,jk,indice))) |
329 |
|
|
C_coef(ig,jk-1,indice) = & |
330 |
|
2295666 |
(ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s |
331 |
|
2303442 |
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 |
✓✓ |
255938 |
DO ig=1,knon |
341 |
|
|
pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* & |
342 |
|
255074 |
(C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1)) |
343 |
|
|
pcapcal(ig) = ztherm_i(ig)* & |
344 |
|
255074 |
(dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1)) |
345 |
|
255074 |
z1s = lambda*(1.-D_coef(ig,1,indice))+1. |
346 |
|
255074 |
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 |
|
255938 |
/ptimestep |
352 |
|
|
ENDDO |
353 |
|
|
|
354 |
|
864 |
END SUBROUTINE soil |