source: LMDZ4/trunk/libf/phylmd/fisrtilp_tr.F @ 741

Last change on this file since 741 was 524, checked in by lmdzadmin, 20 years ago

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