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

Contents of /trunk/phylmd/fisrtilp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 103 - (show annotations)
Fri Aug 29 13:00:05 2014 UTC (9 years, 8 months ago) by guez
File size: 16911 byte(s)
Renamed module cvparam to cv_param. Deleted procedure
cv_param. Changed variables of module cv_param into parameters.

In procedures cv_driver, cv_uncompress and cv3_uncompress, removed
some arguments giving dimensions and used module variables klon and
klev instead.

In procedures gradiv2, laplacien_gam and laplacien, changed
declarations of local variables because klevel is not always klev.

Removed code for nudging surface pressure.

Removed arguments pim and pjm of tau2alpha. Added assignment of false
to variable first.

Replaced real argument del of procedures foeew and FOEDE by logical
argument.

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écipitation stratiforme, schéma de
16 ! nuage, schéma de condensation à grande échelle (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, thermcep
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 ! Fonctions en ligne:
117
118 REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace
119 REAL zzz
120
121 fallvc(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_con
122 fallvs(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
123
124 !---------------------------------------------------------------
125
126 zdelq = 0.0
127
128 IF (appel1er) THEN
129 PRINT *, 'fisrtilp, ninter:', ninter
130 PRINT *, 'fisrtilp, evap_prec:', evap_prec
131 PRINT *, 'fisrtilp, cpartiel:', cpartiel
132 IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
133 PRINT *, "fisrtilp : ce n'est pas prévu, voir Z. X. Li", dtime
134 PRINT *, 'Je préfère un sous-intervalle de 6 minutes.'
135 END IF
136 appel1er = .FALSE.
137
138 ! initialiation provisoire
139 a_tr_sca(1) = -0.5
140 a_tr_sca(2) = -0.5
141 a_tr_sca(3) = -0.5
142 a_tr_sca(4) = -0.5
143
144 ! Initialisation a 1 des coefs des fractions lessivees
145 DO k = 1, klev
146 DO i = 1, klon
147 pfrac_nucl(i, k) = 1.
148 pfrac_1nucl(i, k) = 1.
149 pfrac_impa(i, k) = 1.
150 END DO
151 END DO
152 END IF
153
154 ! Initialisation a 0 de zoliq
155 DO i = 1, klon
156 zoliq(i) = 0.
157 END DO
158 ! Determiner les nuages froids par leur temperature
159 ! nexpo regle la raideur de la transition eau liquide / eau glace.
160
161 ztglace = rtt - 15.0
162 nexpo = 6
163
164 ! Initialiser les sorties:
165
166 DO k = 1, klev + 1
167 DO i = 1, klon
168 prfl(i, k) = 0.0
169 psfl(i, k) = 0.0
170 END DO
171 END DO
172
173 DO k = 1, klev
174 DO i = 1, klon
175 d_t(i, k) = 0.0
176 d_q(i, k) = 0.0
177 d_ql(i, k) = 0.0
178 rneb(i, k) = 0.0
179 radliq(i, k) = 0.0
180 frac_nucl(i, k) = 1.
181 frac_impa(i, k) = 1.
182 END DO
183 END DO
184 DO i = 1, klon
185 rain(i) = 0.0
186 snow(i) = 0.0
187 END DO
188
189 ! Initialiser le flux de precipitation a zero
190
191 DO i = 1, klon
192 zrfl(i) = 0.0
193 zneb(i) = seuil_neb
194 END DO
195
196 ! Pour plus de securite
197
198 zalpha_tr = 0.
199 zfrac_lessi = 0.
200
201 loop_vertical: DO k = klev, 1, -1
202 DO i = 1, klon
203 zt(i) = t(i, k)
204 zq(i) = q(i, k)
205 END DO
206
207 ! Calculer la varition de temp. de l'air du a la chaleur sensible
208 ! transporter par la pluie.
209 ! Il resterait a rajouter cet effet de la chaleur sensible sur les
210 ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
211 ! surface.
212
213 DO i = 1, klon
214 IF (k <= klev - 1) THEN
215 zmair = (paprs(i, k)-paprs(i, k+1))/rg
216 zcpair = rcpd*(1.0+rvtmp2*zq(i))
217 zcpeau = rcpd*rvtmp2
218 zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &
219 * zcpeau + zmair * zcpair* zt(i)) &
220 / (zmair * zcpair + zrfl(i) * dtime * zcpeau)
221 END IF
222 END DO
223
224 IF (evap_prec) THEN
225 ! Calculer l'evaporation de la precipitation
226 DO i = 1, klon
227 IF (zrfl(i)>0.) THEN
228 IF (thermcep) THEN
229 zqs(i) = r2es*foeew(zt(i), rtt >= zt(i))/pplay(i, k)
230 zqs(i) = min(0.5, zqs(i))
231 zcor = 1./(1.-retv*zqs(i))
232 zqs(i) = zqs(i)*zcor
233 ELSE
234 IF (zt(i)<t_coup) THEN
235 zqs(i) = qsats(zt(i))/pplay(i, k)
236 ELSE
237 zqs(i) = qsatl(zt(i))/pplay(i, k)
238 END IF
239 END IF
240 zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
241 zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
242 (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg
243 zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtime/ &
244 (paprs(i, k)-paprs(i, k+1))
245 zqev = min(zqev, zqevt)
246 zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime
247
248 ! pour la glace, on réévapore toute la précip dans la
249 ! couche du dessous la glace venant de la couche du
250 ! dessus est simplement dans la couche du dessous.
251
252 IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
253
254 zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
255 k+1)))*dtime
256 zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
257 k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
258 zrfl(i) = zrfln(i)
259 END IF
260 END DO
261 END IF
262
263 ! Calculer Qs et L/Cp*dQs/dT:
264
265 IF (thermcep) THEN
266 DO i = 1, klon
267 zdelta = rtt >= zt(i)
268 zcvm5 = merge(r5ies*rlstt, r5les*rlvtt, zdelta)
269 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
270 zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
271 zqs(i) = min(0.5, zqs(i))
272 zcor = 1./(1.-retv*zqs(i))
273 zqs(i) = zqs(i)*zcor
274 zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
275 END DO
276 ELSE
277 DO i = 1, klon
278 IF (zt(i)<t_coup) THEN
279 zqs(i) = qsats(zt(i))/pplay(i, k)
280 zdqs(i) = dqsats(zt(i), zqs(i))
281 ELSE
282 zqs(i) = qsatl(zt(i))/pplay(i, k)
283 zdqs(i) = dqsatl(zt(i), zqs(i))
284 END IF
285 END DO
286 END IF
287
288 ! Determiner la condensation partielle et calculer la quantite
289 ! de l'eau condensee:
290
291 IF (cpartiel) THEN
292 ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
293 ! nuageuse a partir des PDF de Sandrine Bony.
294 ! rneb : fraction nuageuse
295 ! zqn : eau totale dans le nuage
296 ! zcond : eau condensee moyenne dans la maille.
297
298 ! on prend en compte le réchauffement qui diminue
299 ! la partie condensée
300
301 ! Version avec les ratqs
302
303 IF (iflag_pdf==0) THEN
304 DO i = 1, klon
305 zdelq = min(ratqs(i, k), 0.99)*zq(i)
306 rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
307 zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
308 END DO
309 ELSE
310 ! Version avec les nouvelles PDFs.
311 DO i = 1, klon
312 IF (zq(i) < 1E-15) THEN
313 zq(i) = 1E-15
314 END IF
315 END DO
316 DO i = 1, klon
317 zpdf_sig(i) = ratqs(i, k)*zq(i)
318 zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
319 zpdf_delta(i) = log(zq(i)/zqs(i))
320 zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
321 zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
322 zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
323 zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
324 zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
325 zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
326 zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
327 zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
328 IF (zpdf_e1(i)<1.E-10) THEN
329 rneb(i, k) = 0.
330 zqn(i) = zqs(i)
331 ELSE
332 rneb(i, k) = 0.5*zpdf_e1(i)
333 zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
334 END IF
335 END DO
336 END IF
337
338 DO i = 1, klon
339 IF (rneb(i, k)<=0.0) zqn(i) = 0.0
340 IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
341 rneb(i, k) = max(0., min(1., rneb(i, k)))
342 ! On ne divise pas par 1 + zdqs pour forcer à avoir l'eau
343 ! prédite par la convection. Attention : il va falloir
344 ! verifier tout ca.
345 zcond(i) = max(0., zqn(i)-zqs(i))*rneb(i, k)
346 rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
347 IF (rneb(i, k) <= 0.) rhcl(i, k) = zq(i) / zqs(i)
348 IF (rneb(i, k) >= 1.) rhcl(i, k) = 1.
349 END DO
350 ELSE
351 DO i = 1, klon
352 IF (zq(i)>zqs(i)) THEN
353 rneb(i, k) = 1.0
354 ELSE
355 rneb(i, k) = 0.0
356 END IF
357 zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
358 END DO
359 END IF
360
361 DO i = 1, klon
362 zq(i) = zq(i) - zcond(i)
363 zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
364 END DO
365
366 ! Partager l'eau condensee en precipitation et eau liquide nuageuse
367
368 DO i = 1, klon
369 IF (rneb(i, k)>0.0) THEN
370 zoliq(i) = zcond(i)
371 zrho(i) = pplay(i, k)/zt(i)/rd
372 zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
373 zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
374 zfice(i) = min(max(zfice(i), 0.0), 1.0)
375 zfice(i) = zfice(i)**nexpo
376 zneb(i) = max(rneb(i, k), seuil_neb)
377 radliq(i, k) = zoliq(i)/real(ninter+1)
378 END IF
379 END DO
380
381 DO n = 1, ninter
382 DO i = 1, klon
383 IF (rneb(i, k)>0.0) THEN
384 zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
385
386 IF (ptconv(i, k)) THEN
387 zcl(i) = cld_lc_con
388 zct(i) = 1./cld_tau_con
389 ELSE
390 zcl(i) = cld_lc_lsc
391 zct(i) = 1./cld_tau_lsc
392 END IF
393 ! quantité d'eau à élminier.
394 zchau(i) = zct(i)*dtime/real(ninter)*zoliq(i)* &
395 (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
396 ! meme chose pour la glace.
397 IF (ptconv(i, k)) THEN
398 zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
399 fallvc(zrhol(i))*zfice(i)
400 ELSE
401 zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
402 fallvs(zrhol(i))*zfice(i)
403 END IF
404 ztot(i) = zchau(i) + zfroi(i)
405 IF (zneb(i)==seuil_neb) ztot(i) = 0.0
406 ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
407 zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
408 radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
409 END IF
410 END DO
411 END DO
412
413 DO i = 1, klon
414 IF (rneb(i, k)>0.0) THEN
415 d_ql(i, k) = zoliq(i)
416 zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
417 * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)
418 END IF
419 IF (zt(i)<rtt) THEN
420 psfl(i, k) = zrfl(i)
421 ELSE
422 prfl(i, k) = zrfl(i)
423 END IF
424 END DO
425
426 ! Calculer les tendances de q et de t:
427 DO i = 1, klon
428 d_q(i, k) = zq(i) - q(i, k)
429 d_t(i, k) = zt(i) - t(i, k)
430 END DO
431
432 ! Calcul du lessivage stratiforme
433 DO i = 1, klon
434 zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)* &
435 (paprs(i, k)-paprs(i, k+1))/rg
436 IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
437 ! lessivage nucleation LMD5 dans la couche elle-meme
438 IF (t(i, k)>=ztglace) THEN
439 zalpha_tr = a_tr_sca(3)
440 ELSE
441 zalpha_tr = a_tr_sca(4)
442 END IF
443 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
444 pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
445 frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
446
447 ! nucleation avec un facteur -1 au lieu de -0.5
448 zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
449 pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
450 END IF
451 END DO
452
453 ! Lessivage par impaction dans les couches en-dessous
454 ! boucle sur i
455 DO kk = k - 1, 1, -1
456 DO i = 1, klon
457 IF (rneb(i, k)>0. .AND. zprec_cond(i)>0.) THEN
458 IF (t(i, kk)>=ztglace) THEN
459 zalpha_tr = a_tr_sca(1)
460 ELSE
461 zalpha_tr = a_tr_sca(2)
462 END IF
463 zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
464 pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
465 frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
466 END IF
467 END DO
468 END DO
469 end DO loop_vertical
470
471 ! Pluie ou neige au sol selon la temperature de la 1ere couche
472
473 DO i = 1, klon
474 IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
475 snow(i) = zrfl(i)
476 zlh_solid(i) = rlstt - rlvtt
477 ELSE
478 rain(i) = zrfl(i)
479 zlh_solid(i) = 0.
480 END IF
481 END DO
482
483 ! For energy conservation: when snow is present, the solification
484 ! latent heat is considered.
485 DO k = 1, klev
486 DO i = 1, klon
487 zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))
488 zmair = (paprs(i, k)-paprs(i, k+1))/rg
489 zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime
490 d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)
491 END DO
492 END DO
493
494 END SUBROUTINE fisrtilp
495
496 end module fisrtilp_m

  ViewVC Help
Powered by ViewVC 1.1.21