source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/fisrtilp.F90 @ 2687

Last change on this file since 2687 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.0 KB
Line 
1!
2! $Id: fisrtilp.F90 2160 2014-11-28 15:36:29Z musat $
3!
4!
5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6     d_t, d_q, d_ql, d_qi, 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     iflag_ice_thermo)
12
13  !
14  USE dimphy
15  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
16  IMPLICIT none
17  !======================================================================
18  ! Auteur(s): Z.X. Li (LMD/CNRS)
19  ! Date: le 20 mars 1995
20  ! Objet: condensation et precipitation stratiforme.
21  !        schema de nuage
22  !======================================================================
23  !======================================================================
24  !ym include "dimensions.h"
25  !ym include "dimphy.h"
26  include "YOMCST.h"
27  include "tracstoke.h"
28  include "fisrtilp.h"
29  include "nuage.h" ! JBM (3/14)
30  include "iniprint.h"
31
32  !
33  ! Arguments:
34  !
35  REAL dtime ! intervalle du temps (s)
36  REAL paprs(klon,klev+1) ! pression a inter-couche
37  REAL pplay(klon,klev) ! pression au milieu de couche
38  REAL t(klon,klev) ! temperature (K)
39  REAL q(klon,klev) ! humidite specifique (kg/kg)
40  REAL d_t(klon,klev) ! incrementation de la temperature (K)
41  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
42  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
43  REAL d_qi(klon,klev) ! incrementation de l'eau glace
44  REAL rneb(klon,klev) ! fraction nuageuse
45  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
46  REAL rhcl(klon,klev) ! humidite relative en ciel clair
47  REAL rain(klon) ! pluies (mm/s)
48  REAL snow(klon) ! neige (mm/s)
49  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
50  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
51  REAL ztv(klon,klev)
52  REAL zqta(klon,klev),fraca(klon,klev)
53  REAL sigma1(klon,klev),sigma2(klon,klev)
54  REAL qltot(klon,klev),ctot(klon,klev)
55  REAL zpspsk(klon,klev),ztla(klon,klev)
56  REAL zthl(klon,klev)
57  REAL ztfondue, qsl, qsi
58
59  logical lognormale(klon)
60  logical ice_thermo
61
62  !AA
63  ! Coeffients de fraction lessivee : pour OFF-LINE
64  !
65  REAL pfrac_nucl(klon,klev)
66  REAL pfrac_1nucl(klon,klev)
67  REAL pfrac_impa(klon,klev)
68  !
69  ! Fraction d'aerosols lessivee par impaction et par nucleation
70  ! POur ON-LINE
71  !
72  REAL frac_impa(klon,klev)
73  REAL frac_nucl(klon,klev)
74  real zct      ,zcl
75  !AA
76  !
77  ! Options du programme:
78  !
79  REAL seuil_neb ! un nuage existe vraiment au-dela
80  PARAMETER (seuil_neb=0.001)
81
82  INTEGER ninter ! sous-intervals pour la precipitation
83  INTEGER ncoreczq 
84  INTEGER iflag_cldcon
85  INTEGER iflag_ice_thermo
86  PARAMETER (ninter=5)
87  LOGICAL evap_prec ! evaporation de la pluie
88  PARAMETER (evap_prec=.TRUE.)
89  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
90  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
91
92  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
93  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
94  real erf   
95  REAL qcloud(klon)
96  !
97  LOGICAL cpartiel ! condensation partielle
98  PARAMETER (cpartiel=.TRUE.)
99  REAL t_coup
100  PARAMETER (t_coup=234.0)
101  !
102  ! Variables locales:
103  !
104  INTEGER i, k, n, kk
105  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
106  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
107  LOGICAL convergence(klon)
108  REAL DDT0
109  PARAMETER (DDT0=.01)
110  INTEGER n_i(klon), iter
111  REAL cste
112 
113  REAL zrfl(klon), zrfln(klon), zqev, zqevt
114  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
115  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
116  REAL zoliqp(klon), zoliqi(klon)
117  REAL zt(klon)
118! JBM (3/14) nexpo is replaced by exposant_glace
119! REAL nexpo ! exponentiel pour glace/eau
120! INTEGER, PARAMETER :: nexpo=6
121  INTEGER exposant_glace_old
122  REAL t_glace_min_old
123  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
124  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
125  REAL zmelt, zpluie, zice, zcondold
126  PARAMETER (ztfondue=278.15)
127  REAL dzfice(klon)
128  !
129  LOGICAL appel1er
130  SAVE appel1er
131  !$OMP THREADPRIVATE(appel1er)
132  !
133  !---------------------------------------------------------------
134  !
135  !AA Variables traceurs:
136  !AA  Provisoire !!! Parametres alpha du lessivage
137  !AA  A priori on a 4 scavenging # possibles
138  !
139  REAL a_tr_sca(4)
140  save a_tr_sca
141  !$OMP THREADPRIVATE(a_tr_sca)
142  !
143  ! Variables intermediaires
144  !
145  REAL zalpha_tr
146  REAL zfrac_lessi
147  REAL zprec_cond(klon)
148  !AA
149! RomP >>> 15 nov 2012
150  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
151! RomP <<<
152  REAL zmair, zcpair, zcpeau
153  !     Pour la conversion eau-neige
154  REAL zlh_solid(klon), zm_solid
155  !IM
156  !ym      INTEGER klevm1
157  !---------------------------------------------------------------
158  !
159  ! Fonctions en ligne:
160  !
161  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
162  REAL zzz
163  include "YOETHF.h"
164  include "FCTTRE.h"
165  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
166  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
167  !
168  DATA appel1er /.TRUE./
169  !ym
170!CR: pour iflag_ice_thermo=2, on active que la convection
171!  ice_thermo = iflag_ice_thermo .GE. 1
172   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
173  zdelq=0.0
174
175  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
176  IF (appel1er) THEN
177     !
178     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
179     WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
180     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
181     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
182        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
183        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
184        !         CALL abort
185     ENDIF
186     appel1er = .FALSE.
187     !
188     !AA initialiation provisoire
189     a_tr_sca(1) = -0.5
190     a_tr_sca(2) = -0.5
191     a_tr_sca(3) = -0.5
192     a_tr_sca(4) = -0.5
193     !
194     !AA Initialisation a 1 des coefs des fractions lessivees
195     !
196     !cdir collapse
197     DO k = 1, klev
198        DO i = 1, klon
199           pfrac_nucl(i,k)=1.
200           pfrac_1nucl(i,k)=1.
201           pfrac_impa(i,k)=1.
202           beta(i,k)=0.  !RomP initialisation
203        ENDDO
204     ENDDO
205
206  ENDIF          !  test sur appel1er
207  !
208  !MAf Initialisation a 0 de zoliq
209  !      DO i = 1, klon
210  !         zoliq(i)=0.
211  !      ENDDO
212  ! Determiner les nuages froids par leur temperature
213  !  nexpo regle la raideur de la transition eau liquide / eau glace.
214  !
215!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
216  IF (iflag_t_glace.EQ.0) THEN
217!   ztglace = RTT - 15.0
218    t_glace_min_old = RTT - 15.0
219    !AJ<
220    IF (ice_thermo) THEN
221!     nexpo = 2
222      exposant_glace_old = 2
223    ELSE
224!     nexpo = 6
225      exposant_glace_old = 6
226    ENDIF
227   
228  ENDIF
229 
230!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
231!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
232!>AJ
233  !cc      nexpo = 1
234  !
235  ! Initialiser les sorties:
236  !
237  !cdir collapse
238  DO k = 1, klev+1
239     DO i = 1, klon
240        prfl(i,k) = 0.0
241        psfl(i,k) = 0.0
242     ENDDO
243  ENDDO
244
245  !cdir collapse
246  DO k = 1, klev
247     DO i = 1, klon
248        d_t(i,k) = 0.0
249        d_q(i,k) = 0.0
250        d_ql(i,k) = 0.0
251        d_qi(i,k) = 0.0
252        rneb(i,k) = 0.0
253        radliq(i,k) = 0.0
254        frac_nucl(i,k) = 1.
255        frac_impa(i,k) = 1.
256     ENDDO
257  ENDDO
258  DO i = 1, klon
259     rain(i) = 0.0
260     snow(i) = 0.0
261     zoliq(i)=0.
262     !     ENDDO
263     !
264     ! Initialiser le flux de precipitation a zero
265     !
266     !     DO i = 1, klon
267     zrfl(i) = 0.0
268     zifl(i) = 0.0
269     zneb(i) = seuil_neb
270  ENDDO
271  !
272  !
273  !AA Pour plus de securite
274
275  zalpha_tr   = 0.
276  zfrac_lessi = 0.
277
278  !AA----------------------------------------------------------
279  !
280  ncoreczq=0
281  ! Boucle verticale (du haut vers le bas)
282  !
283  !IM : klevm1
284  !ym      klevm1=klev-1
285  DO k = klev, 1, -1
286     !
287     !AA----------------------------------------------------------
288     !
289     DO i = 1, klon
290        zt(i)=t(i,k)
291        zq(i)=q(i,k)
292     ENDDO
293     !
294     ! Calculer la varition de temp. de l'air du a la chaleur sensible
295     ! transporter par la pluie.
296     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
297     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
298     ! surface.
299     !
300     IF(k.LE.klevm1) THEN         
301        DO i = 1, klon
302           !IM
303           zmair=(paprs(i,k)-paprs(i,k+1))/RG
304           zcpair=RCPD*(1.0+RVTMP2*zq(i))
305           zcpeau=RCPD*RVTMP2
306           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
307                + zmair*zcpair*zt(i) ) &
308                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
309           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
310        ENDDO
311     ENDIF
312     !
313     !
314     ! Calculer l'evaporation de la precipitation
315     !
316
317
318     ! Calculer l'evaporation de la precipitation
319     !
320
321
322     IF (evap_prec) THEN
323        DO i = 1, klon
324!AJ<
325!!           IF (zrfl(i) .GT.0.) THEN
326           IF (zrfl(i)+zifl(i).GT.0.) THEN
327!>AJ
328              IF (thermcep) THEN
329                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
330                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
331                 zqs(i)=MIN(0.5,zqs(i))
332                 zcor=1./(1.-RETV*zqs(i))
333                 zqs(i)=zqs(i)*zcor
334              ELSE
335                 IF (zt(i) .LT. t_coup) THEN
336                    zqs(i) = qsats(zt(i)) / pplay(i,k)
337                 ELSE
338                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
339                 ENDIF
340              ENDIF
341           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
342        ENDDO
343!AJ<
344       IF (.NOT. ice_thermo) THEN
345        DO i = 1, klon
346!AJ<
347!!           IF (zrfl(i) .GT.0.) THEN
348           IF (zrfl(i)+zifl(i).GT.0.) THEN
349!>AJ
350                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
351                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
352                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
353                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
354                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
355                zqev = MIN (zqev, zqevt)
356                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
357                     /RG/dtime
358       
359                ! pour la glace, on ré-évapore toute la précip dans la
360                ! couche du dessous
361                ! la glace venant de la couche du dessus est simplement
362                ! dans la couche du dessous.
363       
364                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
365       
366                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
367                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
368                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
369                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
370                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
371                zrfl(i) = zrfln(i)
372                zifl(i) = 0.
373           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
374        ENDDO
375!
376       ELSE ! (.NOT. ice_thermo)
377!
378        DO i = 1, klon
379!AJ<
380!!           IF (zrfl(i) .GT.0.) THEN
381           IF (zrfl(i)+zifl(i).GT.0.) THEN
382!>AJ
383     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384     ! Modification de l'évaporation avec la glace
385     ! Différentiation entre précipitation liquide et solide
386     ! On suppose que coef_evai=2*coef_eva
387     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
388     
389         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
390     !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
391
392     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
393     ! On différencie qsat pour l'eau et la glace
394     ! Si zdelta=1. --> glace
395     ! Si zdelta=0. --> eau liquide
396     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397         
398         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
399         qsl= MIN(0.5,qsl)
400         zcor= 1./(1.-RETV*qsl)
401         qsl= qsl*zcor
402         
403         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
404              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
405         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
406              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
407
408         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
409         qsi= MIN(0.5,qsi)
410         zcor= 1./(1.-RETV*qsi)
411         qsi= qsi*zcor
412
413         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
414              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
415         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
416              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
417
418             
419     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420     ! Vérification sur l'évaporation
421     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
422     
423         IF (zqevt+zqevti.GT.zqev0) THEN
424                  zqev=zqev0*zqevt/(zqevt+zqevti)
425                  zqevi=zqev0*zqevti/(zqevt+zqevti)
426             
427         ELSE
428             IF (zqevt+zqevti.GT.0.) THEN
429                  zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
430                  zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
431             ELSE
432             zqev=0.
433             zqevi=0.
434             ENDIF
435         ENDIF
436     
437         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
438                                 /RG/dtime)
439         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
440                                 /RG/dtime)
441         
442     ! Pour la glace, on révapore toute la précip dans la couche du dessous
443     ! la glace venant de la couche du dessus est simplement dans la couche
444     ! du dessous.
445     
446     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
447!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
448         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
449                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
450         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
451                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
452                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
453                  + (zifln(i)-zifl(i)) &
454                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
455                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
456     
457         zrfl(i) = zrfln(i)
458         zifl(i) = zifln(i)
459
460!CR ATTENTION: deplacement de la fonte de la glace
461           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
462           zmelt = MIN(MAX(zmelt,0.),1.)
463           zrfl(i)=zrfl(i)+zmelt*zifl(i)
464           zifl(i)=zifl(i)*(1.-zmelt)
465!           print*,zt(i),'octavio1'
466           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
467                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
468!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
469!fin CR
470
471
472
473           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
474        ENDDO
475
476      ENDIF ! (.NOT. ice_thermo)
477     
478     ENDIF ! (evap_prec)
479     !
480     ! Calculer Qs et L/Cp*dQs/dT:
481     !
482     IF (thermcep) THEN
483        DO i = 1, klon
484           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
485           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
486           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
487           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
488           zqs(i) = MIN(0.5,zqs(i))
489           zcor = 1./(1.-RETV*zqs(i))
490           zqs(i) = zqs(i)*zcor
491           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
492        ENDDO
493     ELSE
494        DO i = 1, klon
495           IF (zt(i).LT.t_coup) THEN
496              zqs(i) = qsats(zt(i))/pplay(i,k)
497              zdqs(i) = dqsats(zt(i),zqs(i))
498           ELSE
499              zqs(i) = qsatl(zt(i))/pplay(i,k)
500              zdqs(i) = dqsatl(zt(i),zqs(i))
501           ENDIF
502        ENDDO
503     ENDIF
504     !
505     ! Determiner la condensation partielle et calculer la quantite
506     ! de l'eau condensee:
507     !
508!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
509!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
510!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
511!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
512!         stop
513!       endif
514
515     IF (cpartiel) THEN
516
517        !        print*,'Dans partiel k=',k
518        !
519        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
520        !   nuageuse a partir des PDF de Sandrine Bony.
521        !   rneb  : fraction nuageuse
522        !   zqn   : eau totale dans le nuage
523        !   zcond : eau condensee moyenne dans la maille.
524        !  on prend en compte le réchauffement qui diminue la partie
525        ! condensee
526        !
527        !   Version avec les raqts
528
529        if (iflag_pdf.eq.0) then
530
531           do i=1,klon
532              zdelq = min(ratqs(i,k),0.99) * zq(i)
533              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
534              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
535           enddo
536
537        else
538           !
539           !   Version avec les nouvelles PDFs.
540           do i=1,klon
541              if(zq(i).lt.1.e-15) then
542                 ncoreczq=ncoreczq+1
543                 zq(i)=1.e-15
544              endif
545           enddo
546
547           if (iflag_cldcon>=5) then
548
549              call cloudth(klon,klev,k,ztv, &
550                   zq,zqta,fraca, &
551                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
552                   ratqs,zqs,t)
553
554              do i=1,klon
555                 rneb(i,k)=ctot(i,k)
556                 zqn(i)=qcloud(i)
557              enddo
558
559           endif
560
561           if (iflag_cldcon <= 4) then
562              lognormale = .true.
563           elseif (iflag_cldcon >= 6) then
564              ! lognormale en l'absence des thermiques
565              lognormale = fraca(:,k) < 1e-10
566           else
567              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
568              ! bi-gaussienne
569              lognormale = .false.
570           end if
571
572!CR: variation de qsat avec T en présence de glace ou non
573!initialisations
574           do i=1,klon
575              DT(i) = 0.
576              n_i(i)=0
577              Tbef(i)=zt(i)
578              qlbef(i)=0.
579           enddo
580
581
582!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
583           if (iflag_fisrtilp_qsat.ge.0) then
584             do iter=1,iflag_fisrtilp_qsat+1
585               
586               do i=1,klon
587!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
588                 convergence(i)=abs(DT(i)).gt.DDT0
589                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
590                 Tbef(i)=Tbef(i)+DT(i)
591                 if (.not.ice_thermo) then   
592                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
593                 else
594                    if (iflag_t_glace.eq.0) then
595                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
596                    else if (iflag_t_glace.eq.1) then
597                    zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
598                    endif
599                 endif
600                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
601                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
602                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
603                 zqs(i) = MIN(0.5,zqs(i))
604                 zcor = 1./(1.-RETV*zqs(i))
605                 zqs(i) = zqs(i)*zcor
606                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
607                 zpdf_sig(i)=ratqs(i,k)*zq(i)
608                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
609                 zpdf_delta(i)=log(zq(i)/zqs(i))
610                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
611                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
612                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
613                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
614                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
615                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
616                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
617                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
618             
619                 if (zpdf_e1(i).lt.1.e-10) then
620                    rneb(i,k)=0.
621                    zqn(i)=zqs(i)
622                 else
623                    rneb(i,k)=0.5*zpdf_e1(i)
624                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
625                 endif
626
627                 endif !convergence
628               enddo ! boucle en i
629
630                 if (.not. ice_thermo) then
631
632                 do i=1,klon
633                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
634
635                 qlbef(i)=max(0.,zqn(i)-zqs(i))
636                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
637                 denom=1.+rneb(i,k)*zdqs(i)
638                 DT(i)=num/denom
639                 n_i(i)=n_i(i)+1
640                 endif
641                 enddo
642
643                 else
644
645!calcul de la fraction de glace
646!CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
647!                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
648!                                     t_glace_max, exposant_glace)
649!                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
650!                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
651!                 zfice(i) = zfice(i)**nexpo
652                 if (iflag_t_glace.eq.1) then
653                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
654                 endif
655
656                 do i=1,klon
657                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
658                 
659                 if (iflag_t_glace.eq.0) then
660                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
661                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
662                    zfice(i) = zfice(i)**exposant_glace_old
663                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
664                 endif
665                 
666                 if (iflag_t_glace.eq.1) then
667                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
668                 endif
669               
670                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
671                    dzfice(i)=0.
672                 endif
673
674                 if (zfice(i).lt.1) then
675                    cste=RLVTT
676                 else
677                    cste=RLSTT
678                 endif
679
680                 qlbef(i)=max(0.,zqn(i)-zqs(i))
681                 num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
682                 denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
683                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
684                 DT(i)=num/denom
685                 n_i(i)=n_i(i)+1
686
687                 endif ! fin convergence
688                 enddo ! fin boucle i
689
690                 endif !ice_thermo
691
692!                 endif
693!               enddo
694 
695         
696             enddo
697           endif
698
699
700        endif ! iflag_pdf
701
702
703!        if (iflag_fisrtilp_qsat.eq.-1) then
704!CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes:
705       IF (1.eq.0) THEN
706       DO i=1,klon
707           IF (rneb(i,k) .LE. 0.0) THEN
708              zqn(i) = 0.0
709              rneb(i,k) = 0.0
710              zcond(i) = 0.0
711              rhcl(i,k)=zq(i)/zqs(i)
712           ELSE IF (rneb(i,k) .GE. 1.0) THEN
713              zqn(i) = zq(i)
714              rneb(i,k) = 1.0                 
715              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
716              rhcl(i,k)=1.0
717           ELSE
718              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
719              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
720           ENDIF
721        ENDDO
722        ENDIF
723
724!        ELSE
725
726        DO i=1,klon
727           IF (rneb(i,k) .LE. 0.0) THEN
728              zqn(i) = 0.0
729              rneb(i,k) = 0.0
730              zcond(i) = 0.0
731              rhcl(i,k)=zq(i)/zqs(i)
732           ELSE IF (rneb(i,k) .GE. 1.0) THEN
733              zqn(i) = zq(i)
734              rneb(i,k) = 1.0
735              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
736              rhcl(i,k)=1.0
737           ELSE
738              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
739              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
740           ENDIF
741        ENDDO
742
743
744!        ENDIF
745
746        !         do i=1,klon
747        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
748        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
749        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
750        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
751        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
752        !c  la convection.
753        !c  ATTENTION !!! Il va falloir verifier tout ca.
754        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
755        !c           print*,'ZDQS ',zdqs(i)
756        !c--Olivier
757        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
758        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
759        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
760        !c--fin
761        !           ENDDO
762     ELSE
763        DO i = 1, klon
764           IF (zq(i).GT.zqs(i)) THEN
765              rneb(i,k) = 1.0
766           ELSE
767              rneb(i,k) = 0.0
768           ENDIF
769           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
770        ENDDO
771     ENDIF
772     !
773     DO i = 1, klon
774        zq(i) = zq(i) - zcond(i)
775        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
776     ENDDO
777!AJ<
778     IF (.NOT. ice_thermo) THEN
779        if (iflag_fisrtilp_qsat.lt.1) then
780           DO i = 1, klon
781             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
782           ENDDO
783        else if (iflag_fisrtilp_qsat.gt.0) then
784           DO i= 1, klon
785             if (lognormale(i)) then
786             zt(i)=Tbef(i)
787             else
788             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
789             endif
790           ENDDO
791        endif
792     ELSE
793         if (iflag_t_glace.eq.1) then
794            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
795         endif
796         if (iflag_fisrtilp_qsat.lt.1) then
797           DO i = 1, klon
798! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
799!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
800!                                     t_glace_max, exposant_glace)
801              if (iflag_t_glace.eq.0) then
802                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
803                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
804                    zfice(i) = zfice(i)**exposant_glace_old
805              endif
806              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
807                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
808           ENDDO
809         else
810           DO i=1, klon
811              if (lognormale(i)) then
812                 zt(i)=Tbef(i)
813              else
814! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
815!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
816!                                     t_glace_max, exposant_glace)
817              if (iflag_t_glace.eq.0) then
818                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
819                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
820                    zfice(i) = zfice(i)**exposant_glace_old
821              endif
822              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
823                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
824              endif
825           ENDDO
826         endif
827!         print*,zt(i),zrfl(i),zifl(i),'temp1'
828       ENDIF
829!>AJ
830     !
831     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
832     !
833     DO i = 1, klon
834        IF (rneb(i,k).GT.0.0) THEN
835           zoliq(i) = zcond(i)
836           zrho(i) = pplay(i,k) / zt(i) / RD
837           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
838        ENDIF
839     ENDDO
840!AJ<
841     IF (.NOT. ice_thermo) THEN
842       IF (iflag_t_glace.EQ.0) THEN
843         DO i = 1, klon
844            IF (rneb(i,k).GT.0.0) THEN
845               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
846               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
847               zfice(i) = zfice(i)**exposant_glace_old
848!              zfice(i) = zfice(i)**nexpo
849         !!      zfice(i)=0.
850            ENDIF
851         ENDDO
852       ELSE ! of IF (iflag_t_glace.EQ.0)
853         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
854!         DO i = 1, klon
855!            IF (rneb(i,k).GT.0.0) THEN
856! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
857!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
858!                                     t_glace_max, exposant_glace)
859!            ENDIF
860!         ENDDO
861       ENDIF
862     ENDIF
863     DO i = 1, klon
864        IF (rneb(i,k).GT.0.0) THEN
865           zneb(i) = MAX(rneb(i,k), seuil_neb)
866     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
867!           print*,zt(i),'fractionglace'
868!>AJ
869           radliq(i,k) = zoliq(i)/REAL(ninter+1)
870        ENDIF
871     ENDDO
872     !
873     DO n = 1, ninter
874        DO i = 1, klon
875           IF (rneb(i,k).GT.0.0) THEN
876              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
877              ! Initialization of zpluie and zice:
878              zpluie=0
879              zice=0
880              IF (zneb(i).EQ.seuil_neb) THEN
881                 ztot = 0.0
882              ELSE
883                 !  quantite d'eau a eliminer: zchau
884                 !  meme chose pour la glace: zfroi
885                 if (ptconv(i,k)) then
886                    zcl   =cld_lc_con
887                    zct   =1./cld_tau_con
888                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
889                         *fallvc(zrhol(i)) * zfice(i)
890                 else
891                    zcl   =cld_lc_lsc
892                    zct   =1./cld_tau_lsc
893                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
894                         *fallvs(zrhol(i)) * zfice(i)
895                 endif
896                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
897                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
898!AJ<
899                 IF (.NOT. ice_thermo) THEN
900                   ztot    = zchau + zfroi
901                 ELSE
902                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
903                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
904                   ztot    = zpluie    + zice
905                 ENDIF
906!>AJ
907                 ztot    = MAX(ztot   ,0.0)
908              ENDIF
909              ztot    = MIN(ztot,zoliq(i))
910!AJ<
911     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
912     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
913              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
914              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
915              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
916!>AJ
917              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
918           ENDIF
919        ENDDO
920     ENDDO
921     !
922         IF (.NOT. ice_thermo) THEN
923       DO i = 1, klon
924         IF (rneb(i,k).GT.0.0) THEN
925           d_ql(i,k) = zoliq(i)
926           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
927                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
928         ENDIF
929       ENDDO
930     ELSE
931       DO i = 1, klon
932         IF (rneb(i,k).GT.0.0) THEN
933!CR on prend en compte la phase glace
934           if (.not.ice_thermo) then
935           d_ql(i,k) = zoliq(i)
936           d_qi(i,k) = 0.
937           else
938           d_ql(i,k) = (1-zfice(i))*zoliq(i)
939           d_qi(i,k) = zfice(i)*zoliq(i)
940           endif
941!AJ<
942           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
943               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
944           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
945                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
946     !      zrfl(i) = zrfl(i)+  zpluie                         &
947     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
948     !      zifl(i) = zifl(i)+  zice                    &
949     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
950
951         ENDIF                     
952       ENDDO
953     ENDIF
954
955!CR: la fonte est faite au debut
956!      IF (ice_thermo) THEN
957!       DO i = 1, klon
958!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
959!           zmelt = MIN(MAX(zmelt,0.),1.)
960!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
961!           zifl(i)=zifl(i)*(1.-zmelt)
962!           print*,zt(i),'octavio1'
963!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
964!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
965!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
966!       ENDDO
967!     ENDIF
968
969       
970     IF (.NOT. ice_thermo) THEN
971       DO i = 1, klon
972         IF (zt(i).LT.RTT) THEN
973           psfl(i,k)=zrfl(i)
974         ELSE
975           prfl(i,k)=zrfl(i)
976         ENDIF
977       ENDDO
978     ELSE
979     ! JAM*************************************************
980     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
981     ! *****************************************************
982
983       DO i = 1, klon
984     !   IF (zt(i).LT.RTT) THEN
985           psfl(i,k)=zifl(i)
986     !   ELSE
987           prfl(i,k)=zrfl(i)
988     !   ENDIF
989!>AJ
990       ENDDO
991     ENDIF
992     !
993     !
994     ! Calculer les tendances de q et de t:
995     !
996     DO i = 1, klon
997        d_q(i,k) = zq(i) - q(i,k)
998        d_t(i,k) = zt(i) - t(i,k)
999     ENDDO
1000     !
1001     !AA--------------- Calcul du lessivage stratiforme  -------------
1002
1003     DO i = 1,klon
1004        !
1005        if(zcond(i).gt.zoliq(i)+1.e-10) then
1006         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1007        else
1008         beta(i,k) = 0.
1009        endif
1010        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1011             * (paprs(i,k)-paprs(i,k+1))/RG
1012        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1013           !AA lessivage nucleation LMD5 dans la couche elle-meme
1014          IF (iflag_t_glace.EQ.0) THEN
1015           if (t(i,k) .GE. t_glace_min_old) THEN
1016              zalpha_tr = a_tr_sca(3)
1017           else
1018              zalpha_tr = a_tr_sca(4)
1019           endif
1020          ELSE ! of IF (iflag_t_glace.EQ.0)
1021           if (t(i,k) .GE. t_glace_min) THEN
1022              zalpha_tr = a_tr_sca(3)
1023           else
1024              zalpha_tr = a_tr_sca(4)
1025           endif
1026          ENDIF
1027           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1028           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1029           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1030           !
1031           ! nucleation avec un facteur -1 au lieu de -0.5
1032           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1033           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1034        ENDIF
1035        !
1036     ENDDO      ! boucle sur i
1037     !
1038     !AA Lessivage par impaction dans les couches en-dessous
1039     DO kk = k-1, 1, -1
1040        DO i = 1, klon
1041           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1042             IF (iflag_t_glace.EQ.0) THEN
1043              if (t(i,kk) .GE. t_glace_min_old) THEN
1044                 zalpha_tr = a_tr_sca(1)
1045              else
1046                 zalpha_tr = a_tr_sca(2)
1047              endif
1048             ELSE ! of IF (iflag_t_glace.EQ.0)
1049              if (t(i,kk) .GE. t_glace_min) THEN
1050                 zalpha_tr = a_tr_sca(1)
1051              else
1052                 zalpha_tr = a_tr_sca(2)
1053              endif
1054             ENDIF
1055              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1056              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1057              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1058           ENDIF
1059        ENDDO
1060     ENDDO
1061     !
1062     !AA----------------------------------------------------------
1063     !                     FIN DE BOUCLE SUR K   
1064  end DO
1065  !
1066  !AA-----------------------------------------------------------
1067  !
1068  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1069  !
1070!CR: si la thermo de la glace est active, on calcule zifl directement
1071  IF (.NOT.ice_thermo) THEN
1072  DO i = 1, klon
1073     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1074!AJ<
1075!        snow(i) = zrfl(i)
1076        snow(i) = zrfl(i)+zifl(i)
1077!>AJ
1078        zlh_solid(i) = RLSTT-RLVTT
1079     ELSE
1080        rain(i) = zrfl(i)
1081        zlh_solid(i) = 0.
1082     ENDIF
1083  ENDDO
1084
1085  ELSE
1086     DO i = 1, klon
1087        snow(i) = zifl(i)
1088        rain(i) = zrfl(i)
1089     ENDDO
1090   
1091   ENDIF
1092  !
1093  ! For energy conservation : when snow is present, the solification
1094  ! latent heat is considered.
1095!CR: si thermo de la glace, neige deja prise en compte
1096  IF (.not.ice_thermo) THEN
1097  DO k = 1, klev
1098     DO i = 1, klon
1099        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
1100        zmair=(paprs(i,k)-paprs(i,k+1))/RG
1101        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1102        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
1103     END DO
1104  END DO
1105  ENDIF
1106  !
1107
1108  if (ncoreczq>0) then
1109     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1110  endif
1111
1112END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.