GCC Code Coverage Report

Directory: ./
File: phys/soil.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 78 100 78.0%
Branches: 53 80 66.2%
Line Branch Exec Source
1 !
2 ! $Header$
3 !
4 911247 SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, &
5 1440 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 2880 REAL, DIMENSION(klon) :: ztherm_i
88 2880 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
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1439 times.
 1440 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/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 IF (is_mpi_root) THEN
121 1 OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 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
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 1 times.
 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
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 1 times.
 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
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 1 times.
 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
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 960 times.
 1440 IF (indice == is_sic) THEN
175
2/2
✓ Branch 0 taken 104835 times.
✓ Branch 1 taken 480 times.
 105315 DO ig = 1, knon
176 105315 ztherm_i(ig) = inertie_sic
177 ENDDO
178
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 480 times.
 480 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
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 480 times.
 960 ELSE IF (indice == is_lic) THEN
187
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 72960 times.
 73440 DO ig = 1, knon
188 72960 ztherm_i(ig) = inertie_lic
189
2/2
✓ Branch 0 taken 55058 times.
✓ Branch 1 taken 17902 times.
 73440 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
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
 480 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
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 247680 times.
 248160 DO ig = 1, knon
199 ! iflag_inertie=0 correspond au cas inertie=constant, comme avant
200
1/2
✓ Branch 0 taken 247680 times.
✗ Branch 1 not taken.
 247680 IF (iflag_inertie==0) THEN
201 247680 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
2/2
✓ Branch 0 taken 76279 times.
✓ Branch 1 taken 171401 times.
 248160 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
2/2
✓ Branch 0 taken 15840 times.
✓ Branch 1 taken 1440 times.
 17280 DO jk=1,nsoilmx
267 17280 zdz2(jk)=dz2(jk)/ptimestep
268 ENDDO
269
270
2/2
✓ Branch 0 taken 425475 times.
✓ Branch 1 taken 1440 times.
 426915 DO ig=1,knon
271 425475 z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
272 C_coef(ig,nsoilmx-1,indice)= &
273 425475 zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
274 426915 D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
275 ENDDO
276
277
2/2
✓ Branch 0 taken 12960 times.
✓ Branch 1 taken 1440 times.
 14400 DO jk=nsoilmx-1,2,-1
278
2/2
✓ Branch 0 taken 3829275 times.
✓ Branch 1 taken 12960 times.
 3843675 DO ig=1,knon
279 z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
280 3829275 *(1.-D_coef(ig,jk,indice)))
281 C_coef(ig,jk-1,indice)= &
282 3829275 (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
283 3842235 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
2/2
✓ Branch 0 taken 425475 times.
✓ Branch 1 taken 1440 times.
 426915 DO ig=1,knon
296 ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/ &
297 426915 (lambda*(1.-D_coef(ig,1,indice))+1.)
298 ENDDO
299
300 ! Other temperatures
301
2/2
✓ Branch 0 taken 14400 times.
✓ Branch 1 taken 1440 times.
 15840 DO jk=1,nsoilmx-1
302
2/2
✓ Branch 0 taken 4254750 times.
✓ Branch 1 taken 14400 times.
 4270590 DO ig=1,knon
303 ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
304 4269150 *ptsoil(ig,jk)
305 ENDDO
306 ENDDO
307
308
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 960 times.
 1440 IF (indice == is_sic) THEN
309
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 104835 times.
 105315 DO ig = 1 , knon
310 105315 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
2/2
✓ Branch 0 taken 425475 times.
✓ Branch 1 taken 1440 times.
 426915 DO ig=1,knon
320 425475 z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
321 425475 C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
322 426915 D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
323 ENDDO
324
325
2/2
✓ Branch 0 taken 12960 times.
✓ Branch 1 taken 1440 times.
 14400 DO jk=nsoilmx-1,2,-1
326
2/2
✓ Branch 0 taken 3829275 times.
✓ Branch 1 taken 12960 times.
 3843675 DO ig=1,knon
327 z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
328 3829275 *(1.-D_coef(ig,jk,indice)))
329 C_coef(ig,jk-1,indice) = &
330 3829275 (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
331 3842235 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
2/2
✓ Branch 0 taken 425475 times.
✓ Branch 1 taken 1440 times.
 426915 DO ig=1,knon
341 pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
342 425475 (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
343 pcapcal(ig) = ztherm_i(ig)* &
344 425475 (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
345 425475 z1s = lambda*(1.-D_coef(ig,1,indice))+1.
346 425475 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 426915 /ptimestep
352 ENDDO
353
354 1440 END SUBROUTINE soil
355