source: LMDZ5/trunk/libf/phylmd/fisrtilp.F90 @ 2353

Last change on this file since 2353 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

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