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

Last change on this file since 17 was 16, checked in by lmdz, 24 years ago

Modif G. Krinner pour acceleration du code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 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(klon)
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
189cAA----------------------------------------------------------
190c
191c Boucle verticale (du haut vers le bas)
192c
193      DO 9999 k = klev, 1, -1
194c
195cAA----------------------------------------------------------
196c
197      DO i = 1, klon
198         zt(i)=t(i,k)
199         zq(i)=q(i,k)
200      ENDDO
201c
202c Calculer l'evaporation de la precipitation
203c
204      IF (evap_prec) THEN
205      DO i = 1, klon
206      IF (zrfl(i) .GT.0.) THEN
207         IF (thermcep) THEN
208           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
209           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
210           zqs(i)=MIN(0.5,zqs(i))
211           zcor=1./(1.-RETV*zqs(i))
212           zqs(i)=zqs(i)*zcor
213         ELSE
214           IF (zt(i) .LT. t_coup) THEN
215              zqs(i) = qsats(zt(i)) / pplay(i,k)
216           ELSE
217              zqs(i) = qsatl(zt(i)) / pplay(i,k)
218           ENDIF
219         ENDIF
220         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
221         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
222     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
223         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
224     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
225         zqev = MIN (zqev, zqevt)
226         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
227     .                            /RG/dtime
228         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
229     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
230         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
231     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
232     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
233         zrfl(i) = zrfln(i)
234      ENDIF
235      ENDDO
236      ENDIF
237c
238c Calculer Qs et L/Cp*dQs/dT:
239c
240      IF (thermcep) THEN
241         DO i = 1, klon
242           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
243           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
244           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
245           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
246           zqs(i) = MIN(0.5,zqs(i))
247           zcor = 1./(1.-RETV*zqs(i))
248           zqs(i) = zqs(i)*zcor
249           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
250         ENDDO
251      ELSE
252         DO i = 1, klon
253            IF (zt(i).LT.t_coup) THEN
254               zqs(i) = qsats(zt(i))/pplay(i,k)
255               zdqs(i) = dqsats(zt(i),zqs(i))
256            ELSE
257               zqs(i) = qsatl(zt(i))/pplay(i,k)
258               zdqs(i) = dqsatl(zt(i),zqs(i))
259            ENDIF
260         ENDDO
261      ENDIF
262c
263c Determiner la condensation partielle et calculer la quantite
264c de l'eau condensee:
265c
266      IF (cpartiel) THEN
267         DO i = 1, klon
268c
269            zx = pplay(i,k)/paprs(i,1)
270            zx = (zx_max-zx)/(zx_max-zx_min)
271            zx = MIN(MAX(zx,0.0),1.0)
272            zx = zx * zx * zx
273            ratqs = zx * (rat_max-rat_min) + rat_min
274            IF (.NOT.calcrat) ratqs=0.2
275c
276            zdelq = ratqs * zq(i)
277            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
278            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
279            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
280            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
281            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
282            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
283         ENDDO
284      ELSE
285         DO i = 1, klon
286            IF (zq(i).GT.zqs(i)) THEN
287               rneb(i,k) = 1.0
288            ELSE
289               rneb(i,k) = 0.0
290            ENDIF
291            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
292         ENDDO
293      ENDIF
294c
295      DO i = 1, klon
296         zq(i) = zq(i) - zcond(i)
297         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
298      ENDDO
299c
300c Partager l'eau condensee en precipitation et eau liquide nuageuse
301c
302      DO i = 1, klon
303      IF (rneb(i,k).GT.0.0) THEN
304         zoliq(i) = zcond(i)
305         zrho(i) = pplay(i,k) / zt(i) / RD
306         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
307         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
308         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
309         zfice(i) = zfice(i)**nexpo
310         zneb(i) = MAX(rneb(i,k), seuil_neb)
311         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
312      ENDIF
313      ENDDO
314c
315      DO n = 1, ninter
316      DO i = 1, klon
317      IF (rneb(i,k).GT.0.0) THEN
318         zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
319     .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
320         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
321         zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
322     .              *fallv(zrhol(i)) * zfice(i)
323         ztot(i) = zchau(i) + zfroi(i)
324         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
325         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
326         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
327         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
328      ENDIF
329      ENDDO
330      ENDDO
331c
332      DO i = 1, klon
333      IF (rneb(i,k).GT.0.0) THEN
334         d_ql(i,k) = zoliq(i)
335         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
336     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
337      ENDIF
338      ENDDO
339c
340c Calculer les tendances de q et de t:
341c
342      DO i = 1, klon
343         d_q(i,k) = zq(i) - q(i,k)
344         d_t(i,k) = zt(i) - t(i,k)
345      ENDDO
346c
347cAA--------------- Calcul du lessivage stratiforme  -------------
348
349      DO i = 1,klon
350c
351         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
352     .                * (paprs(i,k)-paprs(i,k+1))/RG
353         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
354cAA lessivage nucleation LMD5 dans la couche elle-meme
355            if (t(i,k) .GE. ztglace) THEN
356               zalpha_tr = a_tr_sca(3)
357            else
358               zalpha_tr = a_tr_sca(4)
359            endif
360            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
361            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
362            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
363c
364c nucleation avec un facteur -1 au lieu de -0.5
365            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
366            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
367         ENDIF
368c
369      ENDDO      ! boucle sur i
370c
371cAA Lessivage par impaction dans les couches en-dessous
372      DO kk = k-1, 1, -1
373        DO i = 1, klon
374          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
375            if (t(i,kk) .GE. ztglace) THEN
376              zalpha_tr = a_tr_sca(1)
377            else
378              zalpha_tr = a_tr_sca(2)
379            endif
380            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
381            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
382            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
383          ENDIF
384        ENDDO
385      ENDDO
386c
387cAA----------------------------------------------------------
388c                     FIN DE BOUCLE SUR K   
389 9999 CONTINUE
390c
391cAA-----------------------------------------------------------
392c
393c Pluie ou neige au sol selon la temperature de la 1ere couche
394c
395      DO i = 1, klon
396      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
397         snow(i) = zrfl(i)
398      ELSE
399         rain(i) = zrfl(i)
400      ENDIF
401      ENDDO
402c
403      RETURN
404      END
Note: See TracBrowser for help on using the repository browser.