source: LMDZ5/trunk/libf/phylmd/fisrtilp.F90 @ 1742

Last change on this file since 1742 was 1742, checked in by idelkadi, 11 years ago

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.3 KB
Line 
1!
2! $Id: fisrtilp.F90 1742 2013-04-05 11:49:35Z idelkadi $
3!
4!
5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6     d_t, d_q, d_ql, rneb, radliq, rain, snow,          &
7     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
8     frac_impa, frac_nucl, beta,                        &
9     prfl, psfl, rhcl, zqta, fraca,                     &
10     ztv, zpspsk, ztla, zthl, iflag_cldcon)
11
12  !
13  USE dimphy
14  IMPLICIT none
15  !======================================================================
16  ! Auteur(s): Z.X. Li (LMD/CNRS)
17  ! Date: le 20 mars 1995
18  ! Objet: condensation et precipitation stratiforme.
19  !        schema de nuage
20  !======================================================================
21  !======================================================================
22  !ym include "dimensions.h"
23  !ym include "dimphy.h"
24  include "YOMCST.h"
25  include "tracstoke.h"
26  include "fisrtilp.h"
27  include "iniprint.h"
28
29  !
30  ! Arguments:
31  !
32  REAL dtime ! intervalle du temps (s)
33  REAL paprs(klon,klev+1) ! pression a inter-couche
34  REAL pplay(klon,klev) ! pression au milieu de couche
35  REAL t(klon,klev) ! temperature (K)
36  REAL q(klon,klev) ! humidite specifique (kg/kg)
37  REAL d_t(klon,klev) ! incrementation de la temperature (K)
38  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
39  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
40  REAL rneb(klon,klev) ! fraction nuageuse
41  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
42  REAL rhcl(klon,klev) ! humidite relative en ciel clair
43  REAL rain(klon) ! pluies (mm/s)
44  REAL snow(klon) ! neige (mm/s)
45  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
46  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
47  REAL ztv(klon,klev)
48  REAL zqta(klon,klev),fraca(klon,klev)
49  REAL sigma1(klon,klev),sigma2(klon,klev)
50  REAL qltot(klon,klev),ctot(klon,klev)
51  REAL zpspsk(klon,klev),ztla(klon,klev)
52  REAL zthl(klon,klev)
53
54  logical lognormale(klon)
55
56  !AA
57  ! Coeffients de fraction lessivee : pour OFF-LINE
58  !
59  REAL pfrac_nucl(klon,klev)
60  REAL pfrac_1nucl(klon,klev)
61  REAL pfrac_impa(klon,klev)
62  !
63  ! Fraction d'aerosols lessivee par impaction et par nucleation
64  ! POur ON-LINE
65  !
66  REAL frac_impa(klon,klev)
67  REAL frac_nucl(klon,klev)
68  real zct      ,zcl
69  !AA
70  !
71  ! Options du programme:
72  !
73  REAL seuil_neb ! un nuage existe vraiment au-dela
74  PARAMETER (seuil_neb=0.001)
75
76  INTEGER ninter ! sous-intervals pour la precipitation
77  INTEGER ncoreczq 
78  INTEGER iflag_cldcon
79  PARAMETER (ninter=5)
80  LOGICAL evap_prec ! evaporation de la pluie
81  PARAMETER (evap_prec=.TRUE.)
82  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
83  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
84
85  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
86  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
87  real erf   
88  REAL qcloud(klon)
89  !
90  LOGICAL cpartiel ! condensation partielle
91  PARAMETER (cpartiel=.TRUE.)
92  REAL t_coup
93  PARAMETER (t_coup=234.0)
94  !
95  ! Variables locales:
96  !
97  INTEGER i, k, n, kk
98  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
99  REAL zrfl(klon), zrfln(klon), zqev, zqevt
100  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
101  REAL ztglace, zt(klon)
102  INTEGER nexpo ! exponentiel pour glace/eau
103  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
104  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
105  !
106  LOGICAL appel1er
107  SAVE appel1er
108  !$OMP THREADPRIVATE(appel1er)
109  !
110  !---------------------------------------------------------------
111  !
112  !AA Variables traceurs:
113  !AA  Provisoire !!! Parametres alpha du lessivage
114  !AA  A priori on a 4 scavenging # possibles
115  !
116  REAL a_tr_sca(4)
117  save a_tr_sca
118  !$OMP THREADPRIVATE(a_tr_sca)
119  !
120  ! Variables intermediaires
121  !
122  REAL zalpha_tr
123  REAL zfrac_lessi
124  REAL zprec_cond(klon)
125  !AA
126! RomP >>> 15 nov 2012
127  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
128! RomP <<<
129  REAL zmair, zcpair, zcpeau
130  !     Pour la conversion eau-neige
131  REAL zlh_solid(klon), zm_solid
132  !IM
133  !ym      INTEGER klevm1
134  !---------------------------------------------------------------
135  !
136  ! Fonctions en ligne:
137  !
138  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
139  REAL zzz
140  include "YOETHF.h"
141  include "FCTTRE.h"
142  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
143  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
144  !
145  DATA appel1er /.TRUE./
146  !ym
147  zdelq=0.0
148
149  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
150  IF (appel1er) THEN
151     !
152     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
153     WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
154     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
155     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
156        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
157        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
158        !         CALL abort
159     ENDIF
160     appel1er = .FALSE.
161     !
162     !AA initialiation provisoire
163     a_tr_sca(1) = -0.5
164     a_tr_sca(2) = -0.5
165     a_tr_sca(3) = -0.5
166     a_tr_sca(4) = -0.5
167     !
168     !AA Initialisation a 1 des coefs des fractions lessivees
169     !
170     !cdir collapse
171     DO k = 1, klev
172        DO i = 1, klon
173           pfrac_nucl(i,k)=1.
174           pfrac_1nucl(i,k)=1.
175           pfrac_impa(i,k)=1.
176           beta(i,k)=0.  !RomP initialisation
177        ENDDO
178     ENDDO
179
180  ENDIF          !  test sur appel1er
181  !
182  !MAf Initialisation a 0 de zoliq
183  !      DO i = 1, klon
184  !         zoliq(i)=0.
185  !      ENDDO
186  ! Determiner les nuages froids par leur temperature
187  !  nexpo regle la raideur de la transition eau liquide / eau glace.
188  !
189  ztglace = RTT - 15.0
190  nexpo = 6
191  !cc      nexpo = 1
192  !
193  ! Initialiser les sorties:
194  !
195  !cdir collapse
196  DO k = 1, klev+1
197     DO i = 1, klon
198        prfl(i,k) = 0.0
199        psfl(i,k) = 0.0
200     ENDDO
201  ENDDO
202
203  !cdir collapse
204  DO k = 1, klev
205     DO i = 1, klon
206        d_t(i,k) = 0.0
207        d_q(i,k) = 0.0
208        d_ql(i,k) = 0.0
209        rneb(i,k) = 0.0
210        radliq(i,k) = 0.0
211        frac_nucl(i,k) = 1.
212        frac_impa(i,k) = 1.
213     ENDDO
214  ENDDO
215  DO i = 1, klon
216     rain(i) = 0.0
217     snow(i) = 0.0
218     zoliq(i)=0.
219     !     ENDDO
220     !
221     ! Initialiser le flux de precipitation a zero
222     !
223     !     DO i = 1, klon
224     zrfl(i) = 0.0
225     zneb(i) = seuil_neb
226  ENDDO
227  !
228  !
229  !AA Pour plus de securite
230
231  zalpha_tr   = 0.
232  zfrac_lessi = 0.
233
234  !AA----------------------------------------------------------
235  !
236  ncoreczq=0
237  ! Boucle verticale (du haut vers le bas)
238  !
239  !IM : klevm1
240  !ym      klevm1=klev-1
241  DO k = klev, 1, -1
242     !
243     !AA----------------------------------------------------------
244     !
245     DO i = 1, klon
246        zt(i)=t(i,k)
247        zq(i)=q(i,k)
248     ENDDO
249     !
250     ! Calculer la varition de temp. de l'air du a la chaleur sensible
251     ! transporter par la pluie.
252     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
253     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
254     ! surface.
255     !
256     IF(k.LE.klevm1) THEN         
257        DO i = 1, klon
258           !IM
259           zmair=(paprs(i,k)-paprs(i,k+1))/RG
260           zcpair=RCPD*(1.0+RVTMP2*zq(i))
261           zcpeau=RCPD*RVTMP2
262           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
263                + zmair*zcpair*zt(i) ) &
264                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
265           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
266        ENDDO
267     ENDIF
268     !
269     !
270     ! Calculer l'evaporation de la precipitation
271     !
272
273
274     IF (evap_prec) THEN
275        DO i = 1, klon
276           IF (zrfl(i) .GT.0.) THEN
277              IF (thermcep) THEN
278                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
279                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
280                 zqs(i)=MIN(0.5,zqs(i))
281                 zcor=1./(1.-RETV*zqs(i))
282                 zqs(i)=zqs(i)*zcor
283              ELSE
284                 IF (zt(i) .LT. t_coup) THEN
285                    zqs(i) = qsats(zt(i)) / pplay(i,k)
286                 ELSE
287                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
288                 ENDIF
289              ENDIF
290              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
291              zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
292                   * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
293              zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
294                   * RG*dtime/(paprs(i,k)-paprs(i,k+1))
295              zqev = MIN (zqev, zqevt)
296              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
297                   /RG/dtime
298
299              ! pour la glace, on ré-évapore toute la précip dans la
300              ! couche du dessous
301              ! la glace venant de la couche du dessus est simplement
302              ! dans la couche du dessous.
303
304              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
305
306              zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
307                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
308              zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
309                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
310                   * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
311              zrfl(i) = zrfln(i)
312           ENDIF
313        ENDDO
314     ENDIF
315     !
316     ! Calculer Qs et L/Cp*dQs/dT:
317     !
318     IF (thermcep) THEN
319        DO i = 1, klon
320           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
321           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
322           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
323           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
324           zqs(i) = MIN(0.5,zqs(i))
325           zcor = 1./(1.-RETV*zqs(i))
326           zqs(i) = zqs(i)*zcor
327           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
328        ENDDO
329     ELSE
330        DO i = 1, klon
331           IF (zt(i).LT.t_coup) THEN
332              zqs(i) = qsats(zt(i))/pplay(i,k)
333              zdqs(i) = dqsats(zt(i),zqs(i))
334           ELSE
335              zqs(i) = qsatl(zt(i))/pplay(i,k)
336              zdqs(i) = dqsatl(zt(i),zqs(i))
337           ENDIF
338        ENDDO
339     ENDIF
340     !
341     ! Determiner la condensation partielle et calculer la quantite
342     ! de l'eau condensee:
343     !
344
345     IF (cpartiel) THEN
346
347        !        print*,'Dans partiel k=',k
348        !
349        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
350        !   nuageuse a partir des PDF de Sandrine Bony.
351        !   rneb  : fraction nuageuse
352        !   zqn   : eau totale dans le nuage
353        !   zcond : eau condensee moyenne dans la maille.
354        !  on prend en compte le réchauffement qui diminue la partie
355        ! condensee
356        !
357        !   Version avec les raqts
358
359        if (iflag_pdf.eq.0) then
360
361           do i=1,klon
362              zdelq = min(ratqs(i,k),0.99) * zq(i)
363              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
364              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
365           enddo
366
367        else
368           !
369           !   Version avec les nouvelles PDFs.
370           do i=1,klon
371              if(zq(i).lt.1.e-15) then
372                 ncoreczq=ncoreczq+1
373                 zq(i)=1.e-15
374              endif
375           enddo
376
377           if (iflag_cldcon>=5) then
378
379              call cloudth(klon,klev,k,ztv, &
380                   zq,zqta,fraca, &
381                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
382                   ratqs,zqs,t)
383
384              do i=1,klon
385                 rneb(i,k)=ctot(i,k)
386                 zqn(i)=qcloud(i)
387              enddo
388
389           endif
390
391           if (iflag_cldcon <= 4) then
392              lognormale = .true.
393           elseif (iflag_cldcon >= 6) then
394              ! lognormale en l'absence des thermiques
395              lognormale = fraca(:,k) < 1e-10
396           else
397              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
398              ! bi-gaussienne
399              lognormale = .false.
400           end if
401
402           do i=1,klon
403              if (lognormale(i)) then
404                 zpdf_sig(i)=ratqs(i,k)*zq(i)
405                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
406                 zpdf_delta(i)=log(zq(i)/zqs(i))
407                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
408                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
409                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
410                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
411                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
412                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
413                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
414                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
415              endif
416           enddo
417
418           do i=1,klon
419              if (lognormale(i)) then
420                 if (zpdf_e1(i).lt.1.e-10) then
421                    rneb(i,k)=0.
422                    zqn(i)=zqs(i)
423                 else
424                    rneb(i,k)=0.5*zpdf_e1(i)
425                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
426                 endif
427              endif
428
429           enddo
430
431
432        endif ! iflag_pdf
433
434        DO i=1,klon
435           IF (rneb(i,k) .LE. 0.0) THEN
436              zqn(i) = 0.0
437              rneb(i,k) = 0.0
438              zcond(i) = 0.0
439              rhcl(i,k)=zq(i)/zqs(i)
440           ELSE IF (rneb(i,k) .GE. 1.0) THEN
441              zqn(i) = zq(i)
442              rneb(i,k) = 1.0                 
443              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
444              rhcl(i,k)=1.0
445           ELSE
446              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
447              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
448           ENDIF
449        ENDDO
450        !         do i=1,klon
451        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
452        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
453        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
454        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
455        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
456        !c  la convection.
457        !c  ATTENTION !!! Il va falloir verifier tout ca.
458        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
459        !c           print*,'ZDQS ',zdqs(i)
460        !c--Olivier
461        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
462        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
463        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
464        !c--fin
465        !           ENDDO
466     ELSE
467        DO i = 1, klon
468           IF (zq(i).GT.zqs(i)) THEN
469              rneb(i,k) = 1.0
470           ELSE
471              rneb(i,k) = 0.0
472           ENDIF
473           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
474        ENDDO
475     ENDIF
476     !
477     DO i = 1, klon
478        zq(i) = zq(i) - zcond(i)
479        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
480        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
481     ENDDO
482     !
483     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
484     !
485     DO i = 1, klon
486        IF (rneb(i,k).GT.0.0) THEN
487           zoliq(i) = zcond(i)
488           zrho(i) = pplay(i,k) / zt(i) / RD
489           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
490           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
491           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
492           zfice(i) = zfice(i)**nexpo
493           zneb(i) = MAX(rneb(i,k), seuil_neb)
494           radliq(i,k) = zoliq(i)/REAL(ninter+1)
495        ENDIF
496     ENDDO
497     !
498     DO n = 1, ninter
499        DO i = 1, klon
500           IF (rneb(i,k).GT.0.0) THEN
501              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
502
503              IF (zneb(i).EQ.seuil_neb) THEN
504                 ztot = 0.0
505              ELSE
506                 !  quantite d'eau a eliminer: zchau
507                 !  meme chose pour la glace: zfroi
508                 if (ptconv(i,k)) then
509                    zcl   =cld_lc_con
510                    zct   =1./cld_tau_con
511                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
512                         *fallvc(zrhol(i)) * zfice(i)
513                 else
514                    zcl   =cld_lc_lsc
515                    zct   =1./cld_tau_lsc
516                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
517                         *fallvs(zrhol(i)) * zfice(i)
518                 endif
519                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
520                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
521                 ztot    = zchau    + zfroi
522                 ztot    = MAX(ztot   ,0.0)
523              ENDIF
524              ztot    = MIN(ztot,zoliq(i))
525              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
526              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
527           ENDIF
528        ENDDO
529     ENDDO
530     !
531     DO i = 1, klon
532        IF (rneb(i,k).GT.0.0) THEN
533           d_ql(i,k) = zoliq(i)
534           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
535                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
536        ENDIF
537        IF (zt(i).LT.RTT) THEN
538           psfl(i,k)=zrfl(i)
539        ELSE
540           prfl(i,k)=zrfl(i)
541        ENDIF
542     ENDDO
543     !
544     ! Calculer les tendances de q et de t:
545     !
546     DO i = 1, klon
547        d_q(i,k) = zq(i) - q(i,k)
548        d_t(i,k) = zt(i) - t(i,k)
549     ENDDO
550     !
551     !AA--------------- Calcul du lessivage stratiforme  -------------
552
553     DO i = 1,klon
554        !
555        if(zcond(i).gt.zoliq(i)+1.e-10) then
556         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
557        else
558         beta(i,k) = 0.
559        endif
560        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
561             * (paprs(i,k)-paprs(i,k+1))/RG
562        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
563           !AA lessivage nucleation LMD5 dans la couche elle-meme
564           if (t(i,k) .GE. ztglace) THEN
565              zalpha_tr = a_tr_sca(3)
566           else
567              zalpha_tr = a_tr_sca(4)
568           endif
569           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
570           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
571           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
572           !
573           ! nucleation avec un facteur -1 au lieu de -0.5
574           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
575           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
576        ENDIF
577        !
578     ENDDO      ! boucle sur i
579     !
580     !AA Lessivage par impaction dans les couches en-dessous
581     DO kk = k-1, 1, -1
582        DO i = 1, klon
583           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
584              if (t(i,kk) .GE. ztglace) THEN
585                 zalpha_tr = a_tr_sca(1)
586              else
587                 zalpha_tr = a_tr_sca(2)
588              endif
589              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
590              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
591              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
592           ENDIF
593        ENDDO
594     ENDDO
595     !
596     !AA----------------------------------------------------------
597     !                     FIN DE BOUCLE SUR K   
598  end DO
599  !
600  !AA-----------------------------------------------------------
601  !
602  ! Pluie ou neige au sol selon la temperature de la 1ere couche
603  !
604  DO i = 1, klon
605     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
606        snow(i) = zrfl(i)
607        zlh_solid(i) = RLSTT-RLVTT
608     ELSE
609        rain(i) = zrfl(i)
610        zlh_solid(i) = 0.
611     ENDIF
612  ENDDO
613  !
614  ! For energy conservation : when snow is present, the solification
615  ! latent heat is considered.
616  DO k = 1, klev
617     DO i = 1, klon
618        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
619        zmair=(paprs(i,k)-paprs(i,k+1))/RG
620        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
621        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
622     END DO
623  END DO
624  !
625
626  if (ncoreczq>0) then
627     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
628  endif
629
630END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.