source: LMDZ.3.3/trunk/libf/phylmd/fisrtilp_tr.F @ 518

Last change on this file since 518 was 347, checked in by lmdz, 22 years ago

Inclusion humidite relative en ciel clair. O. Boucher
LF

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