source: LMDZ.3.3/branches/LF/libf/phylmd/fisrtilp_tr.F @ 467

Last change on this file since 467 was 2, checked in by lmdz, 25 years ago

Initial revision

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