GCC Code Coverage Report


Directory: ./
File: phys/diagphy.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 122 0.0%
Branches: 0 26 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 ✗ SUBROUTINE diagphy(airephy, tit, iprt, tops, topl, sols, soll, sens, evap, &
5 rain_fall, snow_fall, ts, d_etp_tot, d_qt_tot, d_ec_tot, fs_bound, &
6 fq_bound)
7 ! ======================================================================
8
9 ! Purpose:
10 ! Compute the thermal flux and the watter mass flux at the atmosphere
11 ! boundaries. Print them and also the atmospheric enthalpy change and
12 ! the atmospheric mass change.
13
14 ! Arguments:
15 ! airephy-------input-R- grid area
16 ! tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
17 ! iprt--------input-I- PRINT level ( <=0 : no PRINT)
18 ! tops(klon)--input-R- SW rad. at TOA (W/m2), positive up.
19 ! topl(klon)--input-R- LW rad. at TOA (W/m2), positive down
20 ! sols(klon)--input-R- Net SW flux above surface (W/m2), positive up
21 ! (i.e. -1 * flux absorbed by the surface)
22 ! soll(klon)--input-R- Net LW flux above surface (W/m2), positive up
23 ! (i.e. flux emited - flux absorbed by the surface)
24 ! sens(klon)--input-R- Sensible Flux at surface (W/m2), positive down
25 ! evap(klon)--input-R- Evaporation + sublimation watter vapour mass flux
26 ! (kg/m2/s), positive up
27 ! rain_fall(klon)
28 ! --input-R- Liquid watter mass flux (kg/m2/s), positive down
29 ! snow_fall(klon)
30 ! --input-R- Solid watter mass flux (kg/m2/s), positive down
31 ! ts(klon)----input-R- Surface temperature (K)
32 ! d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
33 ! change (W/m2)
34 ! d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
35 ! change (kg/m2/s)
36 ! d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
37 ! change (W/m2)
38
39 ! fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
40 ! fq_bound---output-R- Watter mass flux at the atmosphere boundaries
41 ! (kg/m2/s)
42
43 ! J.L. Dufresne, July 2002
44 ! Version prise sur
45 ! ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
46 ! le 25 Novembre 2002.
47 ! ======================================================================
48
49 USE dimphy
50 IMPLICIT NONE
51
52 include "YOMCST.h"
53 include "YOETHF.h"
54
55 ! Input variables
56 REAL airephy(klon)
57 CHARACTER *15 tit
58 INTEGER iprt
59 REAL tops(klon), topl(klon), sols(klon), soll(klon)
60 REAL sens(klon), evap(klon), rain_fall(klon), snow_fall(klon)
61 REAL ts(klon)
62 REAL d_etp_tot, d_qt_tot, d_ec_tot
63 ! Output variables
64 REAL fs_bound, fq_bound
65
66 ! Local variables
67 REAL stops, stopl, ssols, ssoll
68 REAL ssens, sfront, slat
69 REAL airetot, zcpvap, zcwat, zcice
70 REAL rain_fall_tot, snow_fall_tot, evap_tot
71
72 INTEGER i
73
74 INTEGER pas
75 SAVE pas
76 DATA pas/0/
77 !$OMP THREADPRIVATE(pas)
78
79 ✗ pas = pas + 1
80 ✗ stops = 0.
81 ✗ stopl = 0.
82 ✗ ssols = 0.
83 ✗ ssoll = 0.
84 ✗ ssens = 0.
85 sfront = 0.
86 ✗ evap_tot = 0.
87 rain_fall_tot = 0.
88 snow_fall_tot = 0.
89 airetot = 0.
90
91 ! Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
92 ! la glace, on travaille par difference a la chaleur specifique de l'
93 ! air sec. En effet, comme on travaille a niveau de pression donne,
94 ! toute variation de la masse d'un constituant est totalement
95 ! compense par une variation de masse d'air.
96
97 ✗ zcpvap = rcpv - rcpd
98 ✗ zcwat = rcw - rcpd
99 ✗ zcice = rcs - rcpd
100
101 ✗ DO i = 1, klon
102 ✗ stops = stops + tops(i)*airephy(i)
103 ✗ stopl = stopl + topl(i)*airephy(i)
104 ✗ ssols = ssols + sols(i)*airephy(i)
105 ✗ ssoll = ssoll + soll(i)*airephy(i)
106 ✗ ssens = ssens + sens(i)*airephy(i)
107 sfront = sfront + (evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice)* &
108 ✗ ts(i)*airephy(i)
109 ✗ evap_tot = evap_tot + evap(i)*airephy(i)
110 ✗ rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
111 ✗ snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
112 ✗ airetot = airetot + airephy(i)
113 END DO
114 ✗ stops = stops/airetot
115 ✗ stopl = stopl/airetot
116 ✗ ssols = ssols/airetot
117 ✗ ssoll = ssoll/airetot
118 ✗ ssens = ssens/airetot
119 ✗ sfront = sfront/airetot
120 ✗ evap_tot = evap_tot/airetot
121 ✗ rain_fall_tot = rain_fall_tot/airetot
122 ✗ snow_fall_tot = snow_fall_tot/airetot
123
124 ✗ slat = rlvtt*rain_fall_tot + rlstt*snow_fall_tot
125 ! Heat flux at atm. boundaries
126 ✗ fs_bound = stops - stopl - (ssols+ssoll) + ssens + sfront + slat
127 ! Watter flux at atm. boundaries
128 ✗ fq_bound = evap_tot - rain_fall_tot - snow_fall_tot
129
130 ✗ IF (iprt>=1) WRITE (6, 6666) tit, pas, fs_bound, d_etp_tot, fq_bound, &
131 ✗ d_qt_tot
132
133 ✗ IF (iprt>=1) WRITE (6, 6668) tit, pas, d_etp_tot + d_ec_tot - fs_bound, &
134 ✗ d_qt_tot - fq_bound
135
136 ✗ IF (iprt>=2) WRITE (6, 6667) tit, pas, stops, stopl, ssols, ssoll, ssens, &
137 ✗ slat, evap_tot, rain_fall_tot + snow_fall_tot
138
139 ✗ RETURN
140
141 6666 FORMAT ('Phys. Flux Budget ', A15, 1I6, 2F8.2, 2(1PE13.5))
142 6667 FORMAT ('Phys. Boundary Flux ', A15, 1I6, 6F8.2, 2(1PE13.5))
143 6668 FORMAT ('Phys. Total Budget ', A15, 1I6, F8.2, 2(1PE13.5))
144
145 END SUBROUTINE diagphy
146
147 ! ======================================================================
148 ✗ SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
149 ✗ u, v, paprs, pplay, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
150 ! ======================================================================
151
152 ! Purpose:
153 ! Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
154 ! et calcul le flux de chaleur et le flux d'eau necessaire a ces
155 ! changements. Ces valeurs sont moyennees sur la surface de tout
156 ! le globe et sont exprime en W/2 et kg/s/m2
157 ! Outil pour diagnostiquer la conservation de l'energie
158 ! et de la masse dans la physique. Suppose que les niveau de
159 ! pression entre couche ne varie pas entre 2 appels.
160
161 ! Plusieurs de ces diagnostics peuvent etre fait en parallele: les
162 ! bilans sont sauvegardes dans des tableaux indices. On parlera
163 ! "d'indice de diagnostic"
164
165
166 ! ======================================================================
167 ! Arguments:
168 ! airephy-------input-R- grid area
169 ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
170 ! iprt----input-I- PRINT level ( <=1 : no PRINT)
171 ! idiag---input-I- indice dans lequel sera range les nouveaux
172 ! bilans d' entalpie et de masse
173 ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse
174 ! sont compare au bilan de d'enthalpie de masse de
175 ! l'indice numero idiag2
176 ! Cas parriculier : si idiag2=0, pas de comparaison, on
177 ! sort directement les bilans d'enthalpie et de masse
178 ! dtime----input-R- time step (s)
179 ! t--------input-R- temperature (K)
180 ! q--------input-R- vapeur d'eau (kg/kg)
181 ! ql-------input-R- liquid watter (kg/kg)
182 ! qs-------input-R- solid watter (kg/kg)
183 ! u--------input-R- vitesse u
184 ! v--------input-R- vitesse v
185 ! paprs----input-R- pression a intercouche (Pa)
186 ! pplay----input-R- pression au milieu de couche (Pa)
187
188 ! the following total value are computed by UNIT of earth surface
189
190 ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
191 ! change (J/m2) during one time step (dtime) for the whole
192 ! atmosphere (air, watter vapour, liquid and solid)
193 ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the
194 ! total watter (kg/m2) change during one time step (dtime),
195 ! d_qw------output-R- same, for the watter vapour only (kg/m2/s)
196 ! d_ql------output-R- same, for the liquid watter only (kg/m2/s)
197 ! d_qs------output-R- same, for the solid watter only (kg/m2/s)
198 ! d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
199
200 ! other (COMMON...)
201 ! RCPD, RCPV, ....
202
203 ! J.L. Dufresne, July 2002
204 ! ======================================================================
205
206 USE dimphy
207 IMPLICIT NONE
208
209 include "YOMCST.h"
210 include "YOETHF.h"
211
212 ! Input variables
213 REAL airephy(klon)
214 CHARACTER *15 tit
215 INTEGER iprt, idiag, idiag2
216 REAL dtime
217 REAL t(klon, klev), q(klon, klev), ql(klon, klev), qs(klon, klev)
218 REAL u(klon, klev), v(klon, klev)
219 REAL paprs(klon, klev+1), pplay(klon, klev)
220 ! Output variables
221 REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
222
223 ! Local variables
224
225 REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot, &
226 qs_tot, ec_tot
227 ! h_vcol_tot-- total enthalpy of vertical air column
228 ! (air with watter vapour, liquid and solid) (J/m2)
229 ! h_dair_tot-- total enthalpy of dry air (J/m2)
230 ! h_qw_tot---- total enthalpy of watter vapour (J/m2)
231 ! h_ql_tot---- total enthalpy of liquid watter (J/m2)
232 ! h_qs_tot---- total enthalpy of solid watter (J/m2)
233 ! qw_tot------ total mass of watter vapour (kg/m2)
234 ! ql_tot------ total mass of liquid watter (kg/m2)
235 ! qs_tot------ total mass of solid watter (kg/m2)
236 ! ec_tot------ total cinetic energy (kg/m2)
237
238 ✗ REAL zairm(klon, klev) ! layer air mass (kg/m2)
239 ✗ REAL zqw_col(klon)
240 ✗ REAL zql_col(klon)
241 ✗ REAL zqs_col(klon)
242 ✗ REAL zec_col(klon)
243 ✗ REAL zh_dair_col(klon)
244 ✗ REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
245
246 REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
247
248 REAL airetot, zcpvap, zcwat, zcice
249
250 INTEGER i, k
251
252 INTEGER ndiag ! max number of diagnostic in parallel
253 PARAMETER (ndiag=10)
254 INTEGER pas(ndiag)
255 SAVE pas
256 DATA pas/ndiag*0/
257 !$OMP THREADPRIVATE(pas)
258
259 REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag), &
260 h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag), &
261 qs_pre(ndiag), ec_pre(ndiag)
262 SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre, h_qs_pre, qw_pre, ql_pre, &
263 qs_pre, ec_pre
264 !$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
265 !$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
266 ! ======================================================================
267
268 ✗ DO k = 1, klev
269 ✗ DO i = 1, klon
270 ! layer air mass
271 ✗ zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
272 END DO
273 END DO
274
275 ! Reset variables
276 ✗ DO i = 1, klon
277 ✗ zqw_col(i) = 0.
278 ✗ zql_col(i) = 0.
279 ✗ zqs_col(i) = 0.
280 ✗ zec_col(i) = 0.
281 ✗ zh_dair_col(i) = 0.
282 ✗ zh_qw_col(i) = 0.
283 ✗ zh_ql_col(i) = 0.
284 ✗ zh_qs_col(i) = 0.
285 END DO
286
287 ✗ zcpvap = rcpv
288 ✗ zcwat = rcw
289 ✗ zcice = rcs
290
291 ! Compute vertical sum for each atmospheric column
292 ! ================================================
293 ✗ DO k = 1, klev
294 ✗ DO i = 1, klon
295 ! Watter mass
296 ✗ zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
297 ✗ zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
298 ✗ zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
299 ! Cinetic Energy
300 ✗ zec_col(i) = zec_col(i) + 0.5*(u(i,k)**2+v(i,k)**2)*zairm(i, k)
301 ! Air enthalpy
302 zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-q(i,k)-ql(i,k)-qs(i,k))* &
303 ✗ zairm(i, k)*t(i, k)
304 ✗ zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k)
305 zh_ql_col(i) = zh_ql_col(i) + zcwat*ql(i, k)*zairm(i, k)*t(i, k) - &
306 ✗ rlvtt*ql(i, k)*zairm(i, k)
307 zh_qs_col(i) = zh_qs_col(i) + zcice*qs(i, k)*zairm(i, k)*t(i, k) - &
308 ✗ rlstt*qs(i, k)*zairm(i, k)
309
310 END DO
311 END DO
312
313 ! Mean over the planete surface
314 ! =============================
315 qw_tot = 0.
316 ql_tot = 0.
317 qs_tot = 0.
318 ec_tot = 0.
319 h_vcol_tot = 0.
320 h_dair_tot = 0.
321 h_qw_tot = 0.
322 h_ql_tot = 0.
323 h_qs_tot = 0.
324 airetot = 0.
325
326 ✗ DO i = 1, klon
327 ✗ qw_tot = qw_tot + zqw_col(i)*airephy(i)
328 ✗ ql_tot = ql_tot + zql_col(i)*airephy(i)
329 ✗ qs_tot = qs_tot + zqs_col(i)*airephy(i)
330 ✗ ec_tot = ec_tot + zec_col(i)*airephy(i)
331 ✗ h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
332 ✗ h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
333 ✗ h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
334 ✗ h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
335 ✗ airetot = airetot + airephy(i)
336 END DO
337
338 ✗ qw_tot = qw_tot/airetot
339 ✗ ql_tot = ql_tot/airetot
340 ✗ qs_tot = qs_tot/airetot
341 ✗ ec_tot = ec_tot/airetot
342 ✗ h_dair_tot = h_dair_tot/airetot
343 ✗ h_qw_tot = h_qw_tot/airetot
344 ✗ h_ql_tot = h_ql_tot/airetot
345 ✗ h_qs_tot = h_qs_tot/airetot
346
347 ✗ h_vcol_tot = h_dair_tot + h_qw_tot + h_ql_tot + h_qs_tot
348
349 ! Compute the change of the atmospheric state compare to the one
350 ! stored in "idiag2", and convert it in flux. THis computation
351 ! is performed IF idiag2 /= 0 and IF it is not the first CALL
352 ! for "idiag"
353 ! ===================================
354
355 ✗ IF ((idiag2>0) .AND. (pas(idiag2)/=0)) THEN
356 ✗ d_h_vcol = (h_vcol_tot-h_vcol_pre(idiag2))/dtime
357 d_h_dair = (h_dair_tot-h_dair_pre(idiag2))/dtime
358 d_h_qw = (h_qw_tot-h_qw_pre(idiag2))/dtime
359 d_h_ql = (h_ql_tot-h_ql_pre(idiag2))/dtime
360 d_h_qs = (h_qs_tot-h_qs_pre(idiag2))/dtime
361 ✗ d_qw = (qw_tot-qw_pre(idiag2))/dtime
362 ✗ d_ql = (ql_tot-ql_pre(idiag2))/dtime
363 ✗ d_qs = (qs_tot-qs_pre(idiag2))/dtime
364 ✗ d_ec = (ec_tot-ec_pre(idiag2))/dtime
365 ✗ d_qt = d_qw + d_ql + d_qs
366 ELSE
367 ✗ d_h_vcol = 0.
368 d_h_dair = 0.
369 d_h_qw = 0.
370 d_h_ql = 0.
371 d_h_qs = 0.
372 ✗ d_qw = 0.
373 ✗ d_ql = 0.
374 ✗ d_qs = 0.
375 ✗ d_ec = 0.
376 ✗ d_qt = 0.
377 END IF
378
379 ✗ IF (iprt>=2) THEN
380 ✗ WRITE (6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
381 9000 FORMAT ('Phys. Watter Mass Budget (kg/m2/s)', A15, 1I6, 10(1PE14.6))
382 ✗ WRITE (6, 9001) tit, pas(idiag), d_h_vcol
383 9001 FORMAT ('Phys. Enthalpy Budget (W/m2) ', A15, 1I6, 10(F8.2))
384 ✗ WRITE (6, 9002) tit, pas(idiag), d_ec
385 9002 FORMAT ('Phys. Cinetic Energy Budget (W/m2) ', A15, 1I6, 10(F8.2))
386 END IF
387
388 ! Store the new atmospheric state in "idiag"
389
390 ✗ pas(idiag) = pas(idiag) + 1
391 ✗ h_vcol_pre(idiag) = h_vcol_tot
392 ✗ h_dair_pre(idiag) = h_dair_tot
393 ✗ h_qw_pre(idiag) = h_qw_tot
394 ✗ h_ql_pre(idiag) = h_ql_tot
395 ✗ h_qs_pre(idiag) = h_qs_tot
396 ✗ qw_pre(idiag) = qw_tot
397 ✗ ql_pre(idiag) = ql_tot
398 ✗ qs_pre(idiag) = qs_tot
399 ✗ ec_pre(idiag) = ec_tot
400
401 ✗ RETURN
402 END SUBROUTINE diagetpq
403