/[lmdze]/trunk/phylmd/fisrtilp.f
ViewVC logotype

Contents of /trunk/phylmd/fisrtilp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 298 - (show annotations)
Thu Jul 26 16:45:51 2018 UTC (5 years, 9 months ago) by guez
File size: 16315 byte(s)
Use directly dtphys from module comconst when possible instead of
having it trickle down through procedure arguments.

1 module fisrtilp_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE fisrtilp(paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, d_ql, rneb, &
8 radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
9 frac_nucl, prfl, psfl, rhcl)
10
11 ! From phylmd/fisrtilp.F, version 1.2, 2004/11/09 16:55:40
12 ! First author: Z. X. Li (LMD/CNRS), 20 mars 1995
13
14 ! Objet : condensation et pr\'ecipitation stratiforme, sch\'ema de
15 ! nuage, sch\'ema de condensation \`a grande \'echelle (pluie).
16
17 USE numer_rec_95, ONLY: nr_erf
18
19 use comconst, only: dtphys
20 USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
21 cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
22 USE dimphy, ONLY: klev, klon
23 USE fcttre, ONLY: foede, foeew
24 USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
25 USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
26
27 REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche
28 REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
29 REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
30 REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
31 LOGICAL, INTENT (IN):: ptconv(klon, klev)
32
33 REAL, INTENT (IN):: ratqs(klon, klev)
34 ! determine la largeur de distribution de vapeur
35
36 REAL, INTENT (out):: d_t(klon, klev) ! incrementation de la temperature (K)
37 REAL, INTENT (out):: d_q(klon, klev) ! incrementation de la vapeur d'eau
38 REAL, INTENT (out):: d_ql(klon, klev) ! incrementation de l'eau liquide
39 REAL, INTENT (out):: rneb(klon, klev) ! fraction nuageuse
40
41 REAL, INTENT (out):: radliq(klon, klev)
42 ! eau liquide utilisee dans rayonnement
43
44 REAL, INTENT (out):: rain(klon) ! pluies (mm/s)
45 REAL, INTENT (out):: snow(klon) ! neige (mm/s)
46
47 ! Coeffients de fraction lessivee :
48 REAL, INTENT (inout):: pfrac_impa(klon, klev)
49 REAL, INTENT (inout):: pfrac_nucl(klon, klev)
50 REAL, INTENT (inout):: pfrac_1nucl(klon, klev)
51
52 ! Fraction d'aerosols lessivee par impaction
53 REAL, INTENT (out):: frac_impa(klon, klev)
54
55 ! Fraction d'aerosols lessivee par nucleation
56 REAL, INTENT (out):: frac_nucl(klon, klev)
57
58 REAL, INTENT (out):: prfl(klon, klev+1)
59 ! flux d'eau precipitante aux interfaces (kg/m2/s)
60
61 REAL, INTENT (out):: psfl(klon, klev+1)
62 ! flux d'eau precipitante aux interfaces (kg/m2/s)
63
64 REAL, INTENT (out):: rhcl(klon, klev) ! humidite relative en ciel clair
65
66 ! Local:
67
68 REAL zct(klon), zcl(klon)
69
70 ! Options du programme:
71
72 REAL seuil_neb ! un nuage existe vraiment au-dela
73 PARAMETER (seuil_neb=0.001)
74
75 INTEGER ninter ! sous-intervals pour la precipitation
76 PARAMETER (ninter=5)
77 LOGICAL evap_prec ! evaporation de la pluie
78 PARAMETER (evap_prec=.TRUE.)
79 REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
80 REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
81
82 LOGICAL cpartiel ! condensation partielle
83 PARAMETER (cpartiel=.TRUE.)
84 REAL t_coup
85 PARAMETER (t_coup=234.0)
86
87 INTEGER i, k, n, kk
88 REAL zqs(klon), zdqs(klon), zcor, zcvm5
89 logical zdelta
90 REAL zrfl(klon), zrfln(klon), zqev, zqevt
91 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
92 REAL ztglace, zt(klon)
93 INTEGER nexpo ! exponentiel pour glace/eau
94 REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
95 REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
96
97 LOGICAL:: appel1er = .TRUE.
98
99 ! Variables traceurs:
100 ! Provisoire !!! Parametres alpha du lessivage
101 ! A priori on a 4 scavenging numbers possibles
102
103 REAL, save:: a_tr_sca(4)
104
105 ! Variables intermediaires
106
107 REAL zalpha_tr
108 REAL zfrac_lessi
109 REAL zprec_cond(klon)
110 REAL zmair, zcpair, zcpeau
111 ! Pour la conversion eau-neige
112 REAL zlh_solid(klon), zm_solid
113
114 !---------------------------------------------------------------
115
116 zdelq = 0.0
117
118 IF (appel1er) THEN
119 PRINT *, 'fisrtilp, ninter:', ninter
120 PRINT *, 'fisrtilp, evap_prec:', evap_prec
121 PRINT *, 'fisrtilp, cpartiel:', cpartiel
122 IF (abs(dtphys / real(ninter) - 360.) > 0.001) THEN
123 PRINT *, "fisrtilp : ce n'est pas pr\'evu, voir Z. X. Li", dtphys
124 PRINT *, "Je pr\'ef\`ere un sous-intervalle de 6 minutes."
125 END IF
126 appel1er = .FALSE.
127
128 ! initialiation provisoire
129 a_tr_sca(1) = -0.5
130 a_tr_sca(2) = -0.5
131 a_tr_sca(3) = -0.5
132 a_tr_sca(4) = -0.5
133
134 ! Initialisation a 1 des coefs des fractions lessivees
135 DO k = 1, klev
136 DO i = 1, klon
137 pfrac_nucl(i, k) = 1.
138 pfrac_1nucl(i, k) = 1.
139 pfrac_impa(i, k) = 1.
140 END DO
141 END DO
142 END IF
143
144 ! Initialisation a 0 de zoliq
145 DO i = 1, klon
146 zoliq(i) = 0.
147 END DO
148 ! Determiner les nuages froids par leur temperature
149 ! nexpo regle la raideur de la transition eau liquide / eau glace.
150
151 ztglace = rtt - 15.0
152 nexpo = 6
153
154 ! Initialiser les sorties:
155
156 DO k = 1, klev + 1
157 DO i = 1, klon
158 prfl(i, k) = 0.0
159 psfl(i, k) = 0.0
160 END DO
161 END DO
162
163 DO k = 1, klev
164 DO i = 1, klon
165 d_t(i, k) = 0.0
166 d_q(i, k) = 0.0
167 d_ql(i, k) = 0.0
168 rneb(i, k) = 0.0
169 radliq(i, k) = 0.0
170 frac_nucl(i, k) = 1.
171 frac_impa(i, k) = 1.
172 END DO
173 END DO
174 DO i = 1, klon
175 rain(i) = 0.0
176 snow(i) = 0.0
177 END DO
178
179 ! Initialiser le flux de precipitation a zero
180
181 DO i = 1, klon
182 zrfl(i) = 0.0
183 zneb(i) = seuil_neb
184 END DO
185
186 ! Pour plus de securite
187
188 zalpha_tr = 0.
189 zfrac_lessi = 0.
190
191 loop_vertical: DO k = klev, 1, -1
192 DO i = 1, klon
193 zt(i) = t(i, k)
194 zq(i) = q(i, k)
195 END DO
196
197 ! Calculer la varition de temp. de l'air du a la chaleur sensible
198 ! transporter par la pluie.
199 ! Il resterait a rajouter cet effet de la chaleur sensible sur les
200 ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
201 ! surface.
202
203 DO i = 1, klon
204 IF (k <= klev - 1) THEN
205 zmair = (paprs(i, k)-paprs(i, k+1))/rg
206 zcpair = rcpd*(1.0+rvtmp2*zq(i))
207 zcpeau = rcpd*rvtmp2
208 zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtphys &
209 * zcpeau + zmair * zcpair* zt(i)) &
210 / (zmair * zcpair + zrfl(i) * dtphys * zcpeau)
211 END IF
212 END DO
213
214 IF (evap_prec) THEN
215 ! Calculer l'evaporation de la precipitation
216 DO i = 1, klon
217 IF (zrfl(i)>0.) THEN
218 zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
219 zqs(i) = min(0.5, zqs(i))
220 zcor = 1./(1.-retv*zqs(i))
221 zqs(i) = zqs(i)*zcor
222 zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
223 zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
224 (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg
225 zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtphys/ &
226 (paprs(i, k)-paprs(i, k+1))
227 zqev = min(zqev, zqevt)
228 zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtphys
229
230 ! pour la glace, on r\'e\'evapore toute la pr\'ecip dans la
231 ! couche du dessous la glace venant de la couche du
232 ! dessus est simplement dans la couche du dessous.
233
234 IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
235
236 zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
237 k+1)))*dtphys
238 zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
239 k+1)))*dtphys*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
240 zrfl(i) = zrfln(i)
241 END IF
242 END DO
243 END IF
244
245 ! Calculer Qs et L/Cp*dQs/dT:
246
247 DO i = 1, klon
248 zdelta = rtt >= zt(i)
249 zcvm5 = merge(r5ies*rlstt, r5les*rlvtt, zdelta)
250 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
251 zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
252 zqs(i) = min(0.5, zqs(i))
253 zcor = 1./(1.-retv*zqs(i))
254 zqs(i) = zqs(i)*zcor
255 zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
256 END DO
257
258 ! Determiner la condensation partielle et calculer la quantite
259 ! de l'eau condensee:
260
261 IF (cpartiel) THEN
262 ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
263 ! nuageuse a partir des PDF de Sandrine Bony.
264 ! rneb : fraction nuageuse
265 ! zqn : eau totale dans le nuage
266 ! zcond : eau condensee moyenne dans la maille.
267
268 ! on prend en compte le r\'echauffement qui diminue
269 ! la partie condens\'ee
270
271 ! Version avec les ratqs
272
273 IF (iflag_pdf==0) THEN
274 DO i = 1, klon
275 zdelq = min(ratqs(i, k), 0.99)*zq(i)
276 rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
277 zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
278 END DO
279 ELSE
280 ! Version avec les nouvelles PDFs.
281 DO i = 1, klon
282 IF (zq(i) < 1E-15) THEN
283 zq(i) = 1E-15
284 END IF
285 END DO
286 DO i = 1, klon
287 zpdf_sig(i) = ratqs(i, k)*zq(i)
288 zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
289 zpdf_delta(i) = log(zq(i)/zqs(i))
290 zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
291 zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
292 zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
293 zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
294 zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
295 zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
296 zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
297 zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
298 IF (zpdf_e1(i)<1.E-10) THEN
299 rneb(i, k) = 0.
300 zqn(i) = zqs(i)
301 ELSE
302 rneb(i, k) = 0.5*zpdf_e1(i)
303 zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
304 END IF
305 END DO
306 END IF
307
308 DO i = 1, klon
309 IF (rneb(i, k)<=0.0) zqn(i) = 0.0
310 IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
311 rneb(i, k) = max(0., min(1., rneb(i, k)))
312 ! On ne divise pas par 1 + zdqs pour forcer \`a avoir l'eau
313 ! pr\'edite par la convection. Attention : il va falloir
314 ! verifier tout ca.
315 zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
316 rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
317 IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
318 IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
319 END DO
320 ELSE
321 DO i = 1, klon
322 IF (zq(i)>zqs(i)) THEN
323 rneb(i, k) = 1.0
324 ELSE
325 rneb(i, k) = 0.0
326 END IF
327 zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
328 END DO
329 END IF
330
331 DO i = 1, klon
332 zq(i) = zq(i) - zcond(i)
333 zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
334 END DO
335
336 ! Partager l'eau condensee en precipitation et eau liquide nuageuse
337
338 DO i = 1, klon
339 IF (rneb(i, k)>0.0) THEN
340 zoliq(i) = zcond(i)
341 zrho(i) = pplay(i, k)/zt(i)/rd
342 zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
343 zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
344 zfice(i) = min(max(zfice(i), 0.0), 1.0)
345 zfice(i) = zfice(i)**nexpo
346 zneb(i) = max(rneb(i, k), seuil_neb)
347 radliq(i, k) = zoliq(i)/real(ninter+1)
348 END IF
349 END DO
350
351 DO n = 1, ninter
352 DO i = 1, klon
353 IF (rneb(i, k)>0.0) THEN
354 zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
355
356 IF (ptconv(i, k)) THEN
357 zcl(i) = cld_lc_con
358 zct(i) = 1./cld_tau_con
359 ELSE
360 zcl(i) = cld_lc_lsc
361 zct(i) = 1./cld_tau_lsc
362 END IF
363 ! quantit\'e d'eau \`a \'eliminer
364 zchau(i) = zct(i)*dtphys/real(ninter)*zoliq(i)* &
365 (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
366 ! m\^eme chose pour la glace
367 IF (ptconv(i, k)) THEN
368 zfroi(i) = dtphys/real(ninter)/zdz(i)*zoliq(i)* &
369 fallvc(zrhol(i))*zfice(i)
370 ELSE
371 zfroi(i) = dtphys/real(ninter)/zdz(i)*zoliq(i)* &
372 fallvs(zrhol(i))*zfice(i)
373 END IF
374 ztot(i) = zchau(i) + zfroi(i)
375 IF (zneb(i)==seuil_neb) ztot(i) = 0.0
376 ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
377 zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
378 radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
379 END IF
380 END DO
381 END DO
382
383 DO i = 1, klon
384 IF (rneb(i, k)>0.0) THEN
385 d_ql(i, k) = zoliq(i)
386 zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
387 * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtphys)
388 END IF
389 IF (zt(i)<rtt) THEN
390 psfl(i, k) = zrfl(i)
391 ELSE
392 prfl(i, k) = zrfl(i)
393 END IF
394 END DO
395
396 ! Calculer les tendances de q et de t :
397 DO i = 1, klon
398 d_q(i, k) = zq(i) - q(i, k)
399 d_t(i, k) = zt(i) - t(i, k)
400 END DO
401
402 ! Calcul du lessivage stratiforme
403 DO i = 1, klon
404 zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) &
405 * (paprs(i, k)-paprs(i, k+1))/rg
406 IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
407 ! lessivage nucleation LMD5 dans la couche elle-meme
408 IF (t(i, k)>=ztglace) THEN
409 zalpha_tr = a_tr_sca(3)
410 ELSE
411 zalpha_tr = a_tr_sca(4)
412 END IF
413 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
414 pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
415 frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
416
417 ! nucleation avec un facteur -1 au lieu de -0.5
418 zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
419 pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
420 END IF
421 END DO
422
423 ! Lessivage par impaction dans les couches en-dessous
424 ! boucle sur i
425 DO kk = k - 1, 1, -1
426 DO i = 1, klon
427 IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
428 IF (t(i, kk)>=ztglace) THEN
429 zalpha_tr = a_tr_sca(1)
430 ELSE
431 zalpha_tr = a_tr_sca(2)
432 END IF
433 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
434 pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
435 frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
436 END IF
437 END DO
438 END DO
439 end DO loop_vertical
440
441 ! Pluie ou neige au sol selon la temperature de la 1ere couche
442
443 DO i = 1, klon
444 IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
445 snow(i) = zrfl(i)
446 zlh_solid(i) = rlstt - rlvtt
447 ELSE
448 rain(i) = zrfl(i)
449 zlh_solid(i) = 0.
450 END IF
451 END DO
452
453 ! For energy conservation: when snow is present, the solification
454 ! latent heat is considered.
455 DO k = 1, klev
456 DO i = 1, klon
457 zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))
458 zmair = (paprs(i, k)-paprs(i, k+1))/rg
459 zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtphys
460 d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)
461 END DO
462 END DO
463
464 contains
465
466 ! vitesse de chute pour cristaux de glace
467
468 REAL function fallvs(zzz)
469 REAL, intent(in):: zzz
470 fallvs = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
471 end function fallvs
472
473 !********************************************************
474
475 real function fallvc(zzz)
476 REAL, intent(in):: zzz
477 fallvc = 3.29/2.0*((zzz)**0.16)*ffallv_con
478 end function fallvc
479
480 END SUBROUTINE fisrtilp
481
482 end module fisrtilp_m

  ViewVC Help
Powered by ViewVC 1.1.21