source: LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp_tr.F @ 224

Last change on this file since 224 was 79, checked in by (none), 25 years ago

This commit was manufactured by cvs2svn to create branch 'rel-LF'.

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