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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 221 - (show annotations)
Thu Apr 20 14:44:47 2017 UTC (7 years ago) by guez
File size: 16334 byte(s)
clcdrag is no longer used in LMDZ. Replaced by cdrag in LMDZ. In cdrag
in LMDZ, zxli is a symbolic constant, false. So removed case zxli true
in LMDZE.

read_sst is called zero (if no ocean point on the whole planet) time or
once per call of physiq. If mod(itap - 1, lmt_pas) == 0 then we have
advanced in time of lmt_pas and deja_lu is necessarily false.

qsat[sl] and dqsat[sl] were never called.

Added output of qsurf in histins, following LMDZ.

Last dummy argument dtime of phystokenc is always the same as first
dummy argument pdtphys, removed dtime.

Removed make rules for nag_xref95, since it does not exist any longer.

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

  ViewVC Help
Powered by ViewVC 1.1.21