My Project
 All Classes Files Functions Variables Macros
fisrtilp.F90
Go to the documentation of this file.
1 !
2 ! $Id: fisrtilp.F90 1746 2013-04-19 13:46:37Z idelkadi $
3 !
4 !
5 SUBROUTINE 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))/(1+zdqs(i))
444  rhcl(i,k)=1.0
445  ELSE
446  zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
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 
630 END SUBROUTINE fisrtilp