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

Last change on this file since 1550 was 1507, checked in by idelkadi, 15 years ago
  1. Introduction de nouvelles facon de calucler les ratqs pour firstilp quand iflag_cldcon>=7 : ratqs devient pronos

tic avec une relaxation vers un ratqs issu de la convection (ratqsc) quand celle-ci est active et une relaxation ve
rs ratqss dans le cas contraire.
Les constantes de temps sont a 3h en dur dans le code (en attendant mieux)
On met un coeff iflag_cldcon/100 devant ratqsc
On plafone le ratqs a 0.5 (la aussi en dur dans le code en attendant de stabiliser une facon de faire).

  1. Introduction de nouveaux diagnostics concernant la separation entre la contribution des thermiques et le reste d

ans fisrtilp (dq/tlscth, dq/tlscst,plulth, plulst)
Correction de la sortie de la hauteur des thermiques zmax_th

  1. Introduction d'une option iflag_wake=2 pour prendre en compte les precipitations issues de fisrtipl dans le forc

age des wake.

  1. Correction d'un bug sur l'identification des points affectes par les thermiques (ptconvth(:,k)=fm_therm(:,k+1)>0

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