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

Contents of /trunk/Sources/phylmd/fisrtilp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 207 - (show annotations)
Thu Sep 1 10:30:53 2016 UTC (7 years, 8 months ago) by guez
File size: 16423 byte(s)
New philosophy on compiler options.

Removed source code for thermcep = f. (Not used in LMDZ either.)

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

  ViewVC Help
Powered by ViewVC 1.1.21