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