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

Last change on this file since 34 was 27, checked in by lmdz, 25 years ago

modif de physiq et fistrlp.F pour recuperer les flux d'eau precipitante (cf Olivier)

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