source: LMDZ5/branches/LMDZ5-dev032011/libf/phylmd/fisrtilp.F90 @ 5080

Last change on this file since 5080 was 1492, checked in by Laurent Fairhead, 14 years ago

Merge of development branch LMDZ5V2.0-dev r1455:r1491 into the trunk.
Validation made locally: restart files are strictly equal between the HEAD of the trunk
and r1491 of LMDZ5V2.0-dev

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