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

Last change on this file since 346 was 207, checked in by lmdz, 24 years ago

petit detail
LF

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