source: lmdz_wrf/WRFV3/lmdz/fisrtilp.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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