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

Last change on this file since 3133 was 517, checked in by lmdzadmin, 21 years ago

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