source: lmdz_wrf/trunk/WRFV3/lmdz/fisrtilp.F90 @ 354

Last change on this file since 354 was 186, checked in by lfita, 10 years ago

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