GCC Code Coverage Report


Directory: ./
File: phys/fisrtilp_tr.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 148 0.0%
Branches: 0 86 0.0%

Line Branch Exec Source
1
2 ! $Id: fisrtilp_tr.F90 2346 2015-08-21 15:13:46Z emillour $
3
4
5 SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, &
6 rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
7 frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical
8 ! properties; aeropt.F)
9
10
11 USE dimphy
12 USE print_control_mod, ONLY: lunout
13 IMPLICIT NONE
14 ! ======================================================================
15 ! Auteur(s): Z.X. Li (LMD/CNRS)
16 ! Date: le 20 mars 1995
17 ! Objet: condensation et precipitation stratiforme.
18 ! schema de nuage
19 ! ======================================================================
20 ! ======================================================================
21 include "YOMCST.h"
22
23 ! Arguments:
24
25 REAL dtime ! intervalle du temps (s)
26 REAL paprs(klon, klev+1) ! pression a inter-couche
27 REAL pplay(klon, klev) ! pression au milieu de couche
28 REAL t(klon, klev) ! temperature (K)
29 REAL q(klon, klev) ! humidite specifique (kg/kg)
30 REAL d_t(klon, klev) ! incrementation de la temperature (K)
31 REAL d_q(klon, klev) ! incrementation de la vapeur d'eau
32 REAL d_ql(klon, klev) ! incrementation de l'eau liquide
33 REAL rneb(klon, klev) ! fraction nuageuse
34 REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements
35 REAL rain(klon) ! pluies (mm/s)
36 REAL snow(klon) ! neige (mm/s)
37 REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
38 REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
39
40 ! jq For aerosol opt properties needed (see aeropt.F)
41 REAL rhcl(klon, klev)
42
43 ! AA
44 ! Coeffients de fraction lessivee : pour OFF-LINE
45
46 REAL pfrac_nucl(klon, klev)
47 REAL pfrac_1nucl(klon, klev)
48 REAL pfrac_impa(klon, klev)
49
50 ! Fraction d'aerosols lessivee par impaction et par nucleation
51 ! POur ON-LINE
52
53 REAL frac_impa(klon, klev)
54 REAL frac_nucl(klon, klev)
55 ! AA
56
57 ! Options du programme:
58
59 REAL seuil_neb ! un nuage existe vraiment au-dela
60 PARAMETER (seuil_neb=0.001)
61 REAL ct ! inverse du temps pour qu'un nuage precipite
62 PARAMETER (ct=1./1800.)
63 REAL cl ! seuil de precipitation
64 PARAMETER (cl=2.6E-4)
65 ! cc PARAMETER (cl=2.3e-4)
66 ! cc PARAMETER (cl=2.0e-4)
67 INTEGER ninter ! sous-intervals pour la precipitation
68 PARAMETER (ninter=5)
69 LOGICAL evap_prec ! evaporation de la pluie
70 PARAMETER (evap_prec=.TRUE.)
71 REAL coef_eva
72 PARAMETER (coef_eva=2.0E-05)
73 LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
74 REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
75 PARAMETER (calcrat=.TRUE.)
76 REAL zx_min, rat_max
77 PARAMETER (zx_min=1.0, rat_max=0.01)
78 REAL zx_max, rat_min
79 PARAMETER (zx_max=0.1, rat_min=0.3)
80 REAL zx
81
82 LOGICAL cpartiel ! condensation partielle
83 PARAMETER (cpartiel=.TRUE.)
84 REAL t_coup
85 PARAMETER (t_coup=234.0)
86
87 ! Variables locales:
88
89 INTEGER i, k, n, kk
90 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
91 REAL zrfl(klon), zrfln(klon), zqev, zqevt
92 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
93 REAL ztglace, zt(klon)
94 INTEGER nexpo ! exponentiel pour glace/eau
95 REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
96 REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
97
98 LOGICAL appel1er
99 SAVE appel1er
100 !$OMP THREADPRIVATE(appel1er)
101
102 ! ---------------------------------------------------------------
103
104 ! AA Variables traceurs:
105 ! AA Provisoire !!! Parametres alpha du lessivage
106 ! AA A priori on a 4 scavenging # possibles
107
108 REAL a_tr_sca(4)
109 SAVE a_tr_sca
110 !$OMP THREADPRIVATE(a_tr_sca)
111
112 ! Variables intermediaires
113
114 REAL zalpha_tr
115 REAL zfrac_lessi
116 REAL zprec_cond(klon)
117 ! AA
118 ! ---------------------------------------------------------------
119
120 ! Fonctions en ligne:
121
122 REAL fallv ! vitesse de chute pour crystaux de glace
123 REAL zzz
124 include "YOETHF.h"
125 include "FCTTRE.h"
126 fallv(zzz) = 3.29/2.0*((zzz)**0.16)
127 ! cc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
128 ! cc fallv (zzz) = 3.29 * ((zzz)**0.16)
129
130 DATA appel1er/.TRUE./
131
132 IF (appel1er) THEN
133
134 WRITE (lunout, *) 'fisrtilp, calcrat:', calcrat
135 WRITE (lunout, *) 'fisrtilp, ninter:', ninter
136 WRITE (lunout, *) 'fisrtilp, evap_prec:', evap_prec
137 WRITE (lunout, *) 'fisrtilp, cpartiel:', cpartiel
138 IF (abs(dtime/real(ninter)-360.0)>0.001) THEN
139 WRITE (lunout, *) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
140 WRITE (lunout, *) 'Je prefere un sous-intervalle de 6 minutes'
141 CALL abort
142 END IF
143 appel1er = .FALSE.
144
145 ! AA initialiation provisoire
146 a_tr_sca(1) = -0.5
147 a_tr_sca(2) = -0.5
148 a_tr_sca(3) = -0.5
149 a_tr_sca(4) = -0.5
150
151 ! AA Initialisation a 1 des coefs des fractions lessivees
152
153 DO k = 1, klev
154 DO i = 1, klon
155 pfrac_nucl(i, k) = 1.
156 pfrac_1nucl(i, k) = 1.
157 pfrac_impa(i, k) = 1.
158 END DO
159 END DO
160
161 END IF ! test sur appel1er
162
163 ! MAf Initialisation a 0 de zoliq
164 DO i = 1, klon
165 zoliq(i) = 0.
166 END DO
167 ! Determiner les nuages froids par leur temperature
168
169 ztglace = rtt - 15.0
170 nexpo = 6
171 ! cc nexpo = 1
172
173 ! Initialiser les sorties:
174
175 DO k = 1, klev + 1
176 DO i = 1, klon
177 prfl(i, k) = 0.0
178 psfl(i, k) = 0.0
179 END DO
180 END DO
181
182 DO k = 1, klev
183 DO i = 1, klon
184 d_t(i, k) = 0.0
185 d_q(i, k) = 0.0
186 d_ql(i, k) = 0.0
187 rneb(i, k) = 0.0
188 radliq(i, k) = 0.0
189 frac_nucl(i, k) = 1.
190 frac_impa(i, k) = 1.
191 END DO
192 END DO
193 DO i = 1, klon
194 rain(i) = 0.0
195 snow(i) = 0.0
196 END DO
197
198 ! Initialiser le flux de precipitation a zero
199
200 DO i = 1, klon
201 zrfl(i) = 0.0
202 zneb(i) = seuil_neb
203 END DO
204
205
206 ! AA Pour plus de securite
207
208 zalpha_tr = 0.
209 zfrac_lessi = 0.
210
211 ! AA----------------------------------------------------------
212
213 ! Boucle verticale (du haut vers le bas)
214
215 DO k = klev, 1, -1
216
217 ! AA----------------------------------------------------------
218
219 DO i = 1, klon
220 zt(i) = t(i, k)
221 zq(i) = q(i, k)
222 END DO
223
224 ! Calculer l'evaporation de la precipitation
225
226 IF (evap_prec) THEN
227 DO i = 1, klon
228 IF (zrfl(i)>0.) THEN
229 IF (thermcep) THEN
230 zdelta = max(0., sign(1.,rtt-zt(i)))
231 zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
232 zqs(i) = min(0.5, zqs(i))
233 zcor = 1./(1.-retv*zqs(i))
234 zqs(i) = zqs(i)*zcor
235 ELSE
236 IF (zt(i)<t_coup) THEN
237 zqs(i) = qsats(zt(i))/pplay(i, k)
238 ELSE
239 zqs(i) = qsatl(zt(i))/pplay(i, k)
240 END IF
241 END IF
242 zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
243 zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
244 (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*zt(i)*rd/rg
245 zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
246 (paprs(i,k)-paprs(i,k+1))
247 zqev = min(zqev, zqevt)
248 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
249 zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
250 k+1)))*dtime
251 zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
252 k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
253 zrfl(i) = zrfln(i)
254 END IF
255 END DO
256 END IF
257
258 ! Calculer Qs et L/Cp*dQs/dT:
259
260 IF (thermcep) THEN
261 DO i = 1, klon
262 zdelta = max(0., sign(1.,rtt-zt(i)))
263 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
264 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
265 zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
266 zqs(i) = min(0.5, zqs(i))
267 zcor = 1./(1.-retv*zqs(i))
268 zqs(i) = zqs(i)*zcor
269 zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
270 END DO
271 ELSE
272 DO i = 1, klon
273 IF (zt(i)<t_coup) THEN
274 zqs(i) = qsats(zt(i))/pplay(i, k)
275 zdqs(i) = dqsats(zt(i), zqs(i))
276 ELSE
277 zqs(i) = qsatl(zt(i))/pplay(i, k)
278 zdqs(i) = dqsatl(zt(i), zqs(i))
279 END IF
280 END DO
281 END IF
282
283 ! Determiner la condensation partielle et calculer la quantite
284 ! de l'eau condensee:
285
286 IF (cpartiel) THEN
287 DO i = 1, klon
288
289 zdelq = ratqs(i, k)*zq(i)
290 rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
291 zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
292 IF (rneb(i,k)<=0.0) zqn(i) = 0.0
293 IF (rneb(i,k)>=1.0) zqn(i) = zq(i)
294 rneb(i, k) = max(0.0, min(1.0,rneb(i,k)))
295 zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)/(1.+zdqs(i))
296
297 ! --Olivier
298 rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
299 IF (rneb(i,k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)
300 IF (rneb(i,k)>=1.0) rhcl(i, k) = 1.0
301 ! --fin
302
303 END DO
304 ELSE
305 DO i = 1, klon
306 IF (zq(i)>zqs(i)) THEN
307 rneb(i, k) = 1.0
308 ELSE
309 rneb(i, k) = 0.0
310 END IF
311 zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
312 END DO
313 END IF
314
315 DO i = 1, klon
316 zq(i) = zq(i) - zcond(i)
317 zt(i) = zt(i) + zcond(i)*rlvtt/rcpd
318 END DO
319
320 ! Partager l'eau condensee en precipitation et eau liquide nuageuse
321
322 DO i = 1, klon
323 IF (rneb(i,k)>0.0) THEN
324 zoliq(i) = zcond(i)
325 zrho(i) = pplay(i, k)/zt(i)/rd
326 zdz(i) = (paprs(i,k)-paprs(i,k+1))/(zrho(i)*rg)
327 zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
328 zfice(i) = min(max(zfice(i),0.0), 1.0)
329 zfice(i) = zfice(i)**nexpo
330 zneb(i) = max(rneb(i,k), seuil_neb)
331 radliq(i, k) = zoliq(i)/real(ninter+1)
332 END IF
333 END DO
334
335 DO n = 1, ninter
336 DO i = 1, klon
337 IF (rneb(i,k)>0.0) THEN
338 zchau(i) = ct*dtime/real(ninter)*zoliq(i)* &
339 (1.0-exp(-(zoliq(i)/zneb(i)/cl)**2))*(1.-zfice(i))
340 zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
341 zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)*fallv(zrhol(i))* &
342 zfice(i)
343 ztot(i) = zchau(i) + zfroi(i)
344 IF (zneb(i)==seuil_neb) ztot(i) = 0.0
345 ztot(i) = min(max(ztot(i),0.0), zoliq(i))
346 zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
347 radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
348 END IF
349 END DO
350 END DO
351
352 DO i = 1, klon
353 IF (rneb(i,k)>0.0) THEN
354 d_ql(i, k) = zoliq(i)
355 zrfl(i) = zrfl(i) + max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k &
356 +1))/(rg*dtime)
357 END IF
358 IF (zt(i)<rtt) THEN
359 psfl(i, k) = zrfl(i)
360 ELSE
361 prfl(i, k) = zrfl(i)
362 END IF
363 END DO
364
365 ! Calculer les tendances de q et de t:
366
367 DO i = 1, klon
368 d_q(i, k) = zq(i) - q(i, k)
369 d_t(i, k) = zt(i) - t(i, k)
370 END DO
371
372 ! AA--------------- Calcul du lessivage stratiforme -------------
373
374 DO i = 1, klon
375
376 zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k+1))/ &
377 rg
378 IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
379 ! AA lessivage nucleation LMD5 dans la couche elle-meme
380 IF (t(i,k)>=ztglace) THEN
381 zalpha_tr = a_tr_sca(3)
382 ELSE
383 zalpha_tr = a_tr_sca(4)
384 END IF
385 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
386 pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
387 frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
388
389 ! nucleation avec un facteur -1 au lieu de -0.5
390 zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
391 pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
392 END IF
393
394 END DO ! boucle sur i
395
396 ! AA Lessivage par impaction dans les couches en-dessous
397 DO kk = k - 1, 1, -1
398 DO i = 1, klon
399 IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
400 IF (t(i,kk)>=ztglace) THEN
401 zalpha_tr = a_tr_sca(1)
402 ELSE
403 zalpha_tr = a_tr_sca(2)
404 END IF
405 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
406 pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
407 frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
408 END IF
409 END DO
410 END DO
411
412 ! AA----------------------------------------------------------
413 ! FIN DE BOUCLE SUR K
414 END DO
415
416 ! AA-----------------------------------------------------------
417
418 ! Pluie ou neige au sol selon la temperature de la 1ere couche
419
420 DO i = 1, klon
421 IF ((t(i,1)+d_t(i,1))<rtt) THEN
422 snow(i) = zrfl(i)
423 ELSE
424 rain(i) = zrfl(i)
425 END IF
426 END DO
427
428 RETURN
429 END SUBROUTINE fisrtilp_tr
430