source: LMDZ6/branches/Ocean_skin/libf/phylmd/fisrtilp.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago

Merge revisions 3427:3600 of trunk into branch Ocean_skin

  • 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: 57.6 KB
Line 
1!
2! $Id: fisrtilp.F90 3605 2019-11-21 15:43:45Z lguez $
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  USE cloudth_mod
18  USE ioipsl_getin_p_mod, ONLY : getin_p
19  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
20  USE phys_local_var_mod, ONLY: rneblsvol
21  ! flag to include modifications to ensure energy conservation (if flag >0)
22  USE add_phys_tend_mod, only : fl_cor_ebil
23  IMPLICIT none
24  !======================================================================
25  ! Auteur(s): Z.X. Li (LMD/CNRS)
26  ! Date: le 20 mars 1995
27  ! Objet: condensation et precipitation stratiforme.
28  !        schema de nuage
29  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
30  !             et ilp (il pleut, L. Li)
31  ! Principales parties:
32  ! P0> Thermalisation des precipitations venant de la couche du dessus
33  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
34  ! P2> Formation du nuage (en k)
35  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
36  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
37  ! les valeurs de T et Q initiales
38  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
39  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
40  ! P2.B> Nuage "tout ou rien"
41  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
42  ! P3> Formation de la precipitation (en k)
43  !======================================================================
44  ! JLD:
45  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
46  ! * fl_cor_ebil doit etre > 0 ;
47  !   fl_cor_ebil= 0 pour reproduire anciens bugs
48  !======================================================================
49  include "YOMCST.h"
50  include "fisrtilp.h"
51  include "nuage.h" ! JBM (3/14)
52
53  !
54  ! Principaux inputs:
55  !
56  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
57  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
58  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
59  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
60  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
61  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
62  INTEGER,                         INTENT(IN)   :: iflag_cld_th
63  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
64  !
65  ! Inputs lies aux thermiques
66  !
67  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
68  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca
69  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
70  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zthl
71  !
72  !  Input/output
73  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
74  !
75  ! Principaux outputs:
76  !
77  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
78  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
79  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
80  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
81  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb ! fraction nuageuse
82  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
83  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
84  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
85  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
86  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
87  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl
88
89  !AA
90  ! Coeffients de fraction lessivee : pour OFF-LINE
91  !
92  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_nucl
93  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_1nucl
94  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_impa
95  !
96  ! Fraction d'aerosols lessivee par impaction et par nucleation
97  ! POur ON-LINE
98  !
99  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
100  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
101  !AA
102  ! --------------------------------------------------------------------------------
103  !
104  ! Options du programme:
105  !
106  REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
107  !$OMP THREADPRIVATE(seuil_neb)
108
109
110  INTEGER ninter ! sous-intervals pour la precipitation
111  PARAMETER (ninter=5)
112  INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
113  !$OMP THREADPRIVATE(iflag_evap_prec)
114  !
115  LOGICAL cpartiel ! condensation partielle
116  PARAMETER (cpartiel=.TRUE.)
117  REAL t_coup
118  PARAMETER (t_coup=234.0)
119  REAL DDT0
120  PARAMETER (DDT0=.01)
121  REAL ztfondue
122  PARAMETER (ztfondue=278.15)
123  ! --------------------------------------------------------------------------------
124  !
125  ! Variables locales:
126  !
127  INTEGER i, k, n, kk
128  INTEGER,save::itap=0
129  !$OMP THREADPRIVATE(itap)
130
131  REAL qsl, qsi
132  real zct      ,zcl
133  INTEGER ncoreczq 
134  REAL ctot(klon,klev)
135  REAL ctot_vol(klon,klev)
136  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
137  REAL zdqsdT_raw(klon)
138  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
139
140  logical lognormale(klon)
141  logical ice_thermo
142  LOGICAL convergence(klon)
143  INTEGER n_i(klon), iter
144  REAL cste
145
146  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
147  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
148  real erf   
149  REAL qcloud(klon)
150 
151  REAL zrfl(klon), zrfln(klon), zqev, zqevt
152  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
153  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
154  REAL zoliqp(klon), zoliqi(klon)
155  REAL zt(klon)
156! JBM (3/14) nexpo is replaced by exposant_glace
157! REAL nexpo ! exponentiel pour glace/eau
158! INTEGER, PARAMETER :: nexpo=6
159  INTEGER exposant_glace_old
160  REAL t_glace_min_old
161  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
162  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
163  REAL zmelt, zpluie, zice
164  REAL dzfice(klon)
165  REAL zsolid
166!!!!
167!  Variables pour Bergeron
168  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
169  REAL zqpreci(klon), zqprecl(klon)
170! Variable pour conservation enegie des precipitations
171  REAL zmqc(klon)
172  !
173  LOGICAL appel1er
174  SAVE appel1er
175  !$OMP THREADPRIVATE(appel1er)
176  !
177! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
178! iflag_oldbug_fisrtilp=1 ajoute le BUG
179  INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
180  !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
181  !---------------------------------------------------------------
182  !
183  !AA Variables traceurs:
184  !AA  Provisoire !!! Parametres alpha du lessivage
185  !AA  A priori on a 4 scavenging # possibles
186  !
187  REAL a_tr_sca(4)
188  save a_tr_sca
189  !$OMP THREADPRIVATE(a_tr_sca)
190  !
191  ! Variables intermediaires
192  !
193  REAL zalpha_tr
194  REAL zfrac_lessi
195  REAL zprec_cond(klon)
196  !AA
197! RomP >>> 15 nov 2012
198  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
199! RomP <<<
200  REAL zmair(klon), zcpair, zcpeau
201  !     Pour la conversion eau-neige
202  REAL zlh_solid(klon), zm_solid
203  !---------------------------------------------------------------
204  !
205  ! Fonctions en ligne:
206  !
207  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
208                     ! (Heymsfield & Donner, 1990)
209  REAL zzz
210
211  include "YOETHF.h"
212  include "FCTTRE.h"
213  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
214  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
215  !
216  DATA appel1er /.TRUE./
217  !ym
218!CR: pour iflag_ice_thermo=2, on active que la convection
219!  ice_thermo = iflag_ice_thermo .GE. 1
220
221  itap=itap+1
222  znebprecip(:)=0.
223
224  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
225  zdelq=0.0
226  ctot_vol(1:klon,1:klev)=0.0
227  rneblsvol(1:klon,1:klev)=0.0
228
229  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
230  IF (appel1er) THEN
231     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
232     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
233     CALL getin_p('seuil_neb',seuil_neb)
234     write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp
235     !
236     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
237     WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
238     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
239     
240     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
241        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
242        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
243        !         CALL abort
244     ENDIF
245     appel1er = .FALSE.
246     !
247     !AA initialiation provisoire
248     a_tr_sca(1) = -0.5
249     a_tr_sca(2) = -0.5
250     a_tr_sca(3) = -0.5
251     a_tr_sca(4) = -0.5
252     !
253     !AA Initialisation a 1 des coefs des fractions lessivees
254     !
255     !cdir collapse
256     DO k = 1, klev
257        DO i = 1, klon
258           pfrac_nucl(i,k)=1.
259           pfrac_1nucl(i,k)=1.
260           pfrac_impa(i,k)=1.
261           beta(i,k)=0.  !RomP initialisation
262        ENDDO
263     ENDDO
264
265  ENDIF          !  test sur appel1er
266  !
267  !MAf Initialisation a 0 de zoliq
268  !      DO i = 1, klon
269  !         zoliq(i)=0.
270  !      ENDDO
271  ! Determiner les nuages froids par leur temperature
272  !  nexpo regle la raideur de la transition eau liquide / eau glace.
273  !
274!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
275  IF (iflag_t_glace.EQ.0) THEN
276!   ztglace = RTT - 15.0
277    t_glace_min_old = RTT - 15.0
278    !AJ<
279    IF (ice_thermo) THEN
280!     nexpo = 2
281      exposant_glace_old = 2
282    ELSE
283!     nexpo = 6
284      exposant_glace_old = 6
285    ENDIF
286   
287  ENDIF
288 
289!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
290!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
291!>AJ
292  !cc      nexpo = 1
293  !
294  ! Initialiser les sorties:
295  !
296  !cdir collapse
297  DO k = 1, klev+1
298     DO i = 1, klon
299        prfl(i,k) = 0.0
300        psfl(i,k) = 0.0
301     ENDDO
302  ENDDO
303
304  !cdir collapse
305  DO k = 1, klev
306     DO i = 1, klon
307        d_t(i,k) = 0.0
308        d_q(i,k) = 0.0
309        d_ql(i,k) = 0.0
310        d_qi(i,k) = 0.0
311        rneb(i,k) = 0.0
312        radliq(i,k) = 0.0
313        frac_nucl(i,k) = 1.
314        frac_impa(i,k) = 1.
315     ENDDO
316  ENDDO
317  DO i = 1, klon
318     rain(i) = 0.0
319     snow(i) = 0.0
320     zoliq(i)=0.
321     !     ENDDO
322     !
323     ! Initialiser le flux de precipitation a zero
324     !
325     !     DO i = 1, klon
326     zrfl(i) = 0.0
327     zifl(i) = 0.0
328     zneb(i) = seuil_neb
329  ENDDO
330  !
331  !
332  !AA Pour plus de securite
333
334  zalpha_tr   = 0.
335  zfrac_lessi = 0.
336
337  !AA==================================================================
338  !
339  ncoreczq=0
340  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
341  !
342  DO k = klev, 1, -1
343     !
344     !AA===============================================================
345     !
346     ! Initialisation temperature et vapeur
347     DO i = 1, klon
348        zt(i)=t(i,k)
349        zq(i)=q(i,k)
350     ENDDO
351     !
352     ! ----------------------------------------------------------------
353     ! P0> Thermalisation des precipitations venant de la couche du dessus
354     ! ----------------------------------------------------------------
355     ! Calculer la varition de temp. de l'air du a la chaleur sensible
356     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
357     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
358     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
359     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
360     ! de l'enthalpie  de la couche avec la temperature
361     ! Variables calculees ou modifiees:
362     !   -  zt: temperature de la cocuhe
363     !   - zmqc: masse de precip qui doit etre thermalisee
364     !
365     IF(k.LE.klevm1) THEN         
366        DO i = 1, klon
367           !IM
368           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
369           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
370           zcpair=RCPD*(1.0+RVTMP2*zq(i))
371           zcpeau=RCPD*RVTMP2
372         if (fl_cor_ebil .GT. 0) then
373           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
374           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
375           ! derniere couche
376           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
377           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
378           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
379                 / (zcpair + zmqc(i)*zcpeau)
380         else ! si on maintient les anciennes erreurs
381           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
382                + zmair(i)*zcpair*zt(i) ) &
383                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
384         end if
385        ENDDO
386     ELSE  ! IF(k.LE.klevm1)
387        DO i = 1, klon
388           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
389           zmqc(i) = 0.
390        ENDDO
391     ENDIF ! end IF(k.LE.klevm1)
392     !
393     ! ----------------------------------------------------------------
394     ! P1> Calcul de l'evaporation de la precipitation
395     ! ----------------------------------------------------------------
396     ! On evapore une partie des precipitations venant de la maille du dessus.
397     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
398     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
399     ! Variables calculees ou modifiees:
400     !   - zrfl et zifl: flux de precip liquide et glace
401     !   - zq, zt: humidite et temperature de la cocuhe
402     !   - zmqc: masse de precip qui doit etre thermalisee
403     !
404     IF (iflag_evap_prec>=1) THEN
405        DO i = 1, klon
406!          S'il y a des precipitations
407           IF (zrfl(i)+zifl(i).GT.0.) THEN
408              ! Calcul du qsat
409              IF (thermcep) THEN
410                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
411                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
412                 zqs(i)=MIN(0.5,zqs(i))
413                 zcor=1./(1.-RETV*zqs(i))
414                 zqs(i)=zqs(i)*zcor
415              ELSE
416                 IF (zt(i) .LT. t_coup) THEN
417                    zqs(i) = qsats(zt(i)) / pplay(i,k)
418                 ELSE
419                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
420                 ENDIF
421              ENDIF
422           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
423        ENDDO
424!AJ<
425
426       IF (.NOT. ice_thermo) THEN
427        DO i = 1, klon
428!          S'il y a des precipitations
429           IF (zrfl(i)+zifl(i).GT.0.) THEN
430                ! Evap max pour ne pas saturer la fraction sous le nuage
431                ! Evap max jusqu'à atteindre la saturation dans la partie
432                ! de la maille qui est sous le nuage de la couche du dessus
433                !!! On ne tient compte de cette fraction que sous une seule
434                !!! couche sous le nuage
435                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
436             ! Ajout de la prise en compte des precip a thermiser
437             ! avec petite reecriture
438             if  (fl_cor_ebil .GT. 0) then ! nouveau
439                ! Calcul de l'evaporation du flux de precip herite
440                !   d'au-dessus
441                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
442                     * zmair(i)/pplay(i,k)*zt(i)*RD
443                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
444
445                ! Seuil pour ne pas saturer la fraction sous le nuage
446                zqev = MIN (zqev, zqevt)
447                ! Nouveau flux de precip
448                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
449                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
450                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
451                  zrfln(i)=0.
452                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
453                END IF
454                ! Nouvelle vapeur
455                zq(i) = zq(i) + zqev
456                zmqc(i) = zmqc(i)-zqev
457                ! Nouvelle temperature (chaleur latente)
458                zt(i) = zt(i) - zqev &
459                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
460!!JLD debut de partie a supprimer a terme
461            else ! if  (fl_cor_ebil .GT. 0)
462                ! Calcul de l'evaporation du flux de precip herite
463                !   d'au-dessus
464                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
465                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
466                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
467                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
468                ! Seuil pour ne pas saturer la fraction sous le nuage
469                zqev = MIN (zqev, zqevt)
470                ! Nouveau flux de precip
471                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
472                     /RG/dtime
473                ! Aucun flux liquide pour T < t_coup
474                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
475                ! Nouvelle vapeur
476                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
477                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
478                ! Nouvelle temperature (chaleur latente)
479                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
480                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
481                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
482              end if ! if  (fl_cor_ebil .GT. 0)
483!!JLD fin de partie a supprimer a terme
484                zrfl(i) = zrfln(i)
485                zifl(i) = 0.
486           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
487        ENDDO
488!
489       ELSE ! (.NOT. ice_thermo)
490!      ================================
491!      Avec thermodynamique de la glace
492!      ================================
493        DO i = 1, klon
494!AJ<
495!        S'il y a des precipitations
496         IF (zrfl(i)+zifl(i).GT.0.) THEN
497
498         IF (iflag_evap_prec==1) THEN
499            znebprecip(i)=zneb(i)
500         ELSE
501            znebprecip(i)=MAX(zneb(i),znebprecip(i))
502         ENDIF
503     
504        ! Evap max pour ne pas saturer la fraction sous le nuage
505         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
506
507         !JAM
508         ! On differencie qsat pour l'eau et la glace
509         ! Si zdelta=1. --> glace
510         ! Si zdelta=0. --> eau liquide
511       
512         ! Calcul du qsat par rapport a l'eau liquide
513         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
514         qsl= MIN(0.5,qsl)
515         zcor= 1./(1.-RETV*qsl)
516         qsl= qsl*zcor
517         
518         ! Calcul de l'evaporation du flux de precip venant du dessus
519         ! Formulation en racine du flux de precip
520         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
521         IF (iflag_evap_prec==3) THEN
522         zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
523              *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
524              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
525         ELSE
526         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
527              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
528         ENDIF
529
530
531         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
532              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
533         
534         ! Calcul du qsat par rapport a la glace
535         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
536         qsi= MIN(0.5,qsi)
537         zcor= 1./(1.-RETV*qsi)
538         qsi= qsi*zcor
539
540         ! Calcul de la sublimation du flux de precip solide herite
541         !   d'au-dessus
542         IF (iflag_evap_prec==3) THEN
543         zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
544              *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
545              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
546         ELSE
547         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
548              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
549         ENDIF
550         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
551              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
552
553        !JAM
554        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
555        ! la fraction de la couche sous le nuage sinon on repartit zqev0
556        ! en conservant la proportion liquide / glace
557     
558         IF (zqevt+zqevti.GT.zqev0) THEN
559            zqev=zqev0*zqevt/(zqevt+zqevti)
560            zqevi=zqev0*zqevti/(zqevt+zqevti)
561         ELSE
562!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
563!       liquides et solides meme si on ne sature pas la couche.
564!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
565!            zqev=zqevt
566!            zqevi=zqevti
567             IF (zqevt+zqevti.GT.0.) THEN
568                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
569                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
570             ELSE
571             zqev=0.
572             zqevi=0.
573             ENDIF
574         ENDIF
575         ! Nouveaux flux de precip liquide et solide
576         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
577                                 /RG/dtime)
578         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
579                                 /RG/dtime)
580
581         ! Mise a jour de la vapeur, temperature et flux de precip
582         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
583                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
584       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
585         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
586                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
587         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
588                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
589                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
590                  + (zifln(i)-zifl(i)) &
591                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
592                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
593       else ! sans correction thermalisation des precips
594         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
595                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
596                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
597                  + (zifln(i)-zifl(i)) &
598                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
599                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
600       end if
601        ! Nouvelles vaeleurs des precips liquides et solides
602         zrfl(i) = zrfln(i)
603         zifl(i) = zifln(i)
604!        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
605
606!CR ATTENTION: deplacement de la fonte de la glace
607!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
608!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
609!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
610           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
611           ! fraction de la precip solide qui est fondue
612           zmelt = MIN(MAX(zmelt,0.),1.)
613           ! Fusion de la glace
614           zrfl(i)=zrfl(i)+zmelt*zifl(i)
615           if (fl_cor_ebil .LE. 0) then
616             ! the following line should not be here. Indeed, if zifl is modified
617             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
618             ! and therefore the change in temperature computed below is wrong
619             zifl(i)=zifl(i)*(1.-zmelt)
620           end if
621           ! Chaleur latente de fusion
622        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
623           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
624                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
625        else ! sans correction thermalisation des precips
626           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
627                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
628        end if
629           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
630             zifl(i)=zifl(i)*(1.-zmelt)
631           end if
632
633           ELSE
634              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
635              ! sous nuageuse utilisee pour la pluie.
636              znebprecip(i)=0.
637           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
638        ENDDO
639
640      ENDIF ! (.NOT. ice_thermo)
641     
642     ! ----------------------------------------------------------------
643     ! Fin evaporation de la precipitation
644     ! ----------------------------------------------------------------
645     ENDIF ! (iflag_evap_prec>=1)
646     !
647     ! Calculer Qs et L/Cp*dQs/dT:
648     !
649     IF (thermcep) THEN
650        DO i = 1, klon
651           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
652           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
653       if  (fl_cor_ebil .GT. 0) then ! nouveau
654           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
655       else   
656           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
657       endif
658           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
659           zqs(i) = MIN(0.5,zqs(i))
660           zcor = 1./(1.-RETV*zqs(i))
661           zqs(i) = zqs(i)*zcor
662           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
663           zdqsdT_raw(i) = zdqs(i)*  &
664         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
665        ENDDO
666     ELSE
667        DO i = 1, klon
668           IF (zt(i).LT.t_coup) THEN
669              zqs(i) = qsats(zt(i))/pplay(i,k)
670              zdqs(i) = dqsats(zt(i),zqs(i))
671           ELSE
672              zqs(i) = qsatl(zt(i))/pplay(i,k)
673              zdqs(i) = dqsatl(zt(i),zqs(i))
674           ENDIF
675        ENDDO
676     ENDIF
677     !
678     ! Determiner la condensation partielle et calculer la quantite
679     ! de l'eau condensee:
680     !
681!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
682!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
683!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
684!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
685!         stop
686!       endif
687
688     ! ----------------------------------------------------------------
689     ! P2> Formation du nuage
690     ! ----------------------------------------------------------------
691     ! Variables calculees:
692     !   rneb  : fraction nuageuse
693     !   zcond : eau condensee moyenne dans la maille.
694     !   rhcl: humidite relative ciel-clair
695     !   zt : temperature de la maille
696     ! ----------------------------------------------------------------
697     !
698     IF (cpartiel) THEN
699        ! -------------------------
700        ! P2.A> Nuage fractionnaire
701        ! -------------------------
702        !
703        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
704        !   nuageuse a partir des PDF de Sandrine Bony.
705        !   rneb  : fraction nuageuse
706        !   zqn   : eau totale dans le nuage
707        !   zcond : eau condensee moyenne dans la maille.
708        !  on prend en compte le réchauffement qui diminue la partie
709        ! condensee
710        !
711        !   Version avec les raqts
712
713        ! ----------------------------------------------------------------
714        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
715        ! ----------------------------------------------------------------
716        if (iflag_pdf.eq.0) then
717
718           ! version creneau de (Li, 1998)
719           do i=1,klon
720              zdelq = min(ratqs(i,k),0.99) * zq(i)
721              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
722              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
723           enddo
724
725        else !  if (iflag_pdf.eq.0)
726           ! ----------------------------------------------------------------
727           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
728           ! les valeurs de T et Q initiales
729           ! ----------------------------------------------------------------
730           do i=1,klon
731              if(zq(i).lt.1.e-15) then
732                 ncoreczq=ncoreczq+1
733                 zq(i)=1.e-15
734              endif
735           enddo
736
737           if (iflag_cld_th>=5) then
738
739              if (iflag_cloudth_vert<=2) then
740               call cloudth(klon,klev,k,ztv, &
741                   zq,zqta,fraca, &
742                   qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
743                   ratqs,zqs,t)
744              elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
745               call cloudth_v3(klon,klev,k,ztv, &
746                   zq,zqta,fraca, &
747                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
748                   ratqs,zqs,t)
749              !----------------------------------
750              !Version these Jean Jouhaud, Decembre 2018
751              !----------------------------------             
752              elseif (iflag_cloudth_vert==6) then
753               call cloudth_v6(klon,klev,k,ztv, &
754                   zq,zqta,fraca, &
755                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
756                   ratqs,zqs,t)
757
758              endif
759              do i=1,klon
760                 rneb(i,k)=ctot(i,k)
761                 rneblsvol(i,k)=ctot_vol(i,k)
762                 zqn(i)=qcloud(i)
763              enddo
764
765           endif
766
767           if (iflag_cld_th <= 4) then
768              lognormale = .true.
769           elseif (iflag_cld_th >= 6) then
770              ! lognormale en l'absence des thermiques
771              lognormale = fraca(:,k) < 1e-10
772           else
773              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
774              ! bi-gaussienne
775              lognormale = .false.
776           end if
777
778!CR: variation de qsat avec T en presence de glace ou non
779!initialisations
780           do i=1,klon
781              DT(i) = 0.
782              n_i(i)=0
783              Tbef(i)=zt(i)
784              qlbef(i)=0.
785           enddo
786
787        ! ----------------------------------------------------------------
788        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
789        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
790        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
791        ! qui en dependent. Ce changement de temperature est provisoire, et
792        ! la valeur definitive sera calcule plus tard.
793        ! Variables calculees:
794        !   rneb : nebulosite
795        !   zcond: eau condensee en moyenne dans la maille
796        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
797        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
798        ! T sur qsat
799        ! ----------------------------------------------------------------
800
801!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
802           if (iflag_fisrtilp_qsat.ge.0) then
803             ! Iteration pour condensation avec variation de qsat(T)
804             ! -----------------------------------------------------
805             do iter=1,iflag_fisrtilp_qsat+1
806               
807               do i=1,klon
808!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
809                 ! !! convergence = .true. tant que l'on n'a pas converge !!
810                 ! ------------------------------
811                 convergence(i)=abs(DT(i)).gt.DDT0
812                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
813                 ! si on n'a pas converge
814                 !
815                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
816                 ! ---------------------------------------------------------------
817                 ! Variables calculees:
818                 ! rneb : nebulosite
819                 ! zqn : eau condensee, dans le nuage (in cloud water content)
820                 ! zcond: eau condensee en moyenne dans la maille
821                 ! rhcl: humidite relative ciel-clair
822                 !
823                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
824                 if (.not.ice_thermo) then   
825                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
826                 else
827                    if (iflag_t_glace.eq.0) then
828                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
829                    else if (iflag_t_glace.ge.1) then
830                       if (iflag_oldbug_fisrtilp.EQ.0) then
831                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
832                       else
833                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
834                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
835                       endif
836                    endif
837                 endif
838                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
839                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
840               if (fl_cor_ebil .GT. 0) then
841                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
842               else
843                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
844               end if
845                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
846                 zqs(i) = MIN(0.5,zqs(i))
847                 zcor = 1./(1.-RETV*zqs(i))
848                 zqs(i) = zqs(i)*zcor
849                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
850                 zpdf_sig(i)=ratqs(i,k)*zq(i)
851                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
852                 zpdf_delta(i)=log(zq(i)/zqs(i))
853                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
854                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
855                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
856                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
857                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
858                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
859                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
860                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
861             
862                 if (zpdf_e1(i).lt.1.e-10) then
863                    rneb(i,k)=0.
864                    zqn(i)=zqs(i)
865                 else
866                    rneb(i,k)=0.5*zpdf_e1(i)
867                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
868                 endif
869
870                 ! If vertical heterogeneity, change fraction by volume as well
871                 if (iflag_cloudth_vert>=3) then
872                   ctot_vol(i,k)=rneb(i,k)
873                   rneblsvol(i,k)=ctot_vol(i,k)
874                 endif
875
876                 endif !convergence
877
878               enddo ! boucle en i
879
880                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
881                 !         due a la condensation.
882                 ! ---------------------------------------------------------------
883                 ! Variables calculees:
884                 ! DT : variation de temperature due a la condensation
885
886                 if (.not. ice_thermo) then
887                 ! --------------------------
888                 do i=1,klon
889                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
890
891                 qlbef(i)=max(0.,zqn(i)-zqs(i))
892               if (fl_cor_ebil .GT. 0) then
893                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
894               else
895                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
896               end if
897                 denom=1.+rneb(i,k)*zdqs(i)
898                 DT(i)=num/denom
899                 n_i(i)=n_i(i)+1
900                 endif
901                 enddo
902
903                 else ! if (.not. ice_thermo)
904                 ! --------------------------
905                 if (iflag_t_glace.ge.1) then
906                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
907                 endif
908
909                 do i=1,klon
910                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
911                 
912                 if (iflag_t_glace.eq.0) then
913                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
914                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
915                    zfice(i) = zfice(i)**exposant_glace_old
916                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
917          &                     / (t_glace_min_old - RTT)
918                 endif
919                 
920                 if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
921                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
922          &                    / (t_glace_min - t_glace_max)
923                 endif
924               
925                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
926                    dzfice(i)=0.
927                 endif
928
929                 if (zfice(i).lt.1) then
930                    cste=RLVTT
931                 else
932                    cste=RLSTT
933                 endif
934
935                 qlbef(i)=max(0.,zqn(i)-zqs(i))
936               if (fl_cor_ebil .GT. 0) then
937                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
938           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
939                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
940                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
941           &               *qlbef(i)*dzfice(i)
942               else
943                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
944           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
945                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
946                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
947               end if
948                 DT(i)=num/denom
949                 n_i(i)=n_i(i)+1
950
951                 endif ! fin convergence
952                 enddo ! fin boucle i
953
954                 endif !ice_thermo
955
956             enddo ! iter=1,iflag_fisrtilp_qsat+1
957             ! Fin d'iteration pour condensation avec variation de qsat(T)
958             ! -----------------------------------------------------------
959           endif !  if (iflag_fisrtilp_qsat.ge.0)
960     ! ----------------------------------------------------------------
961     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
962     ! ----------------------------------------------------------------
963
964        endif ! iflag_pdf
965
966!        if (iflag_fisrtilp_qsat.eq.-1) then
967!------------------------------------------
968!CR: ATTENTION option fausse mais a existe:
969! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
970! activer les lignes suivantes:
971       IF (1.eq.0) THEN
972       DO i=1,klon
973           IF (rneb(i,k) .LE. 0.0) THEN
974              zqn(i) = 0.0
975              rneb(i,k) = 0.0
976              zcond(i) = 0.0
977              rhcl(i,k)=zq(i)/zqs(i)
978           ELSE IF (rneb(i,k) .GE. 1.0) THEN
979              zqn(i) = zq(i)
980              rneb(i,k) = 1.0                 
981              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
982              rhcl(i,k)=1.0
983           ELSE
984              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
985              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
986           ENDIF
987       ENDDO
988       ENDIF
989!------------------------------------------
990
991!        ELSE
992        ! ----------------------------------------------------------------
993        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
994        ! Variables calculees:
995        !   rneb : nebulosite
996        !   zcond: eau condensee en moyenne dans la maille
997        !   zq : eau vapeur dans la maille
998        !   zt : temperature de la maille
999        !   rhcl: humidite relative ciel-clair
1000        ! ----------------------------------------------------------------
1001        !
1002        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
1003        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
1004        ! et de l'humidite relative ciel-clair (rhcl)
1005        DO i=1,klon
1006           IF (rneb(i,k) .LE. 0.0) THEN
1007              zqn(i) = 0.0
1008              rneb(i,k) = 0.0
1009              zcond(i) = 0.0
1010              rhcl(i,k)=zq(i)/zqs(i)
1011           ELSE IF (rneb(i,k) .GE. 1.0) THEN
1012              zqn(i) = zq(i)
1013              rneb(i,k) = 1.0
1014              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1015              rhcl(i,k)=1.0
1016           ELSE
1017              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1018              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1019           ENDIF
1020        ENDDO
1021        ! If vertical heterogeneity, change fraction by volume as well
1022        if (iflag_cloudth_vert>=3) then
1023          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1024          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1025        endif
1026
1027!        ENDIF
1028
1029     ELSE ! de IF (cpartiel)
1030        ! -------------------------
1031        ! P2.B> Nuage "tout ou rien"
1032        ! -------------------------
1033        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
1034        DO i = 1, klon
1035           IF (zq(i).GT.zqs(i)) THEN
1036              rneb(i,k) = 1.0
1037           ELSE
1038              rneb(i,k) = 0.0
1039           ENDIF
1040           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1041        ENDDO
1042     ENDIF ! de IF (cpartiel)
1043     !
1044     ! Mise a jour vapeur d'eau
1045     ! -------------------------
1046     DO i = 1, klon
1047        zq(i) = zq(i) - zcond(i)
1048        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1049     ENDDO
1050!AJ<
1051     ! ------------------------------------
1052     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
1053     ! -------------------------------------
1054     ! Variable calcule:
1055     !   zt : temperature de la maille
1056     !
1057     IF (.NOT. ice_thermo) THEN
1058        if (iflag_fisrtilp_qsat.lt.1) then
1059           DO i = 1, klon
1060             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1061           ENDDO
1062        else if (iflag_fisrtilp_qsat.gt.0) then
1063           DO i= 1, klon
1064    if (fl_cor_ebil .GT. 0) then
1065             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1066    else
1067             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1068    end if
1069           ENDDO
1070        endif
1071     ELSE
1072         if (iflag_t_glace.ge.1) then
1073            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1074         endif
1075         if (iflag_fisrtilp_qsat.lt.1) then
1076           DO i = 1, klon
1077! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1078!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1079!                                     t_glace_max, exposant_glace)
1080              if (iflag_t_glace.eq.0) then
1081                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1082                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1083                    zfice(i) = zfice(i)**exposant_glace_old
1084              endif
1085              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1086                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1087           ENDDO
1088         else
1089           DO i=1, klon
1090! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1091!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1092!                                     t_glace_max, exposant_glace)
1093              if (iflag_t_glace.eq.0) then
1094                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1095                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1096                    zfice(i) = zfice(i)**exposant_glace_old
1097              endif
1098        if (fl_cor_ebil .GT. 0) then
1099              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1100           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1101                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1102        else
1103              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1104                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1105        end if
1106           ENDDO
1107         endif
1108!         print*,zt(i),zrfl(i),zifl(i),'temp1'
1109       ENDIF
1110!>AJ
1111
1112     ! ----------------------------------------------------------------
1113     ! P3> Formation des precipitations
1114     ! ----------------------------------------------------------------
1115     !
1116     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
1117     !
1118
1119     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
1120     DO i = 1, klon
1121        IF (rneb(i,k).GT.0.0) THEN
1122           zoliq(i) = zcond(i)
1123           zrho(i) = pplay(i,k) / zt(i) / RD
1124           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
1125        ENDIF
1126     ENDDO
1127!AJ<
1128     IF (.NOT. ice_thermo) THEN
1129       IF (iflag_t_glace.EQ.0) THEN
1130         DO i = 1, klon
1131            IF (rneb(i,k).GT.0.0) THEN
1132               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1133               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1134               zfice(i) = zfice(i)**exposant_glace_old
1135!              zfice(i) = zfice(i)**nexpo
1136         !!      zfice(i)=0.
1137            ENDIF
1138         ENDDO
1139       ELSE ! of IF (iflag_t_glace.EQ.0)
1140         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1141!         DO i = 1, klon
1142!            IF (rneb(i,k).GT.0.0) THEN
1143! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1144!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1145!                                     t_glace_max, exposant_glace)
1146!            ENDIF
1147!         ENDDO
1148       ENDIF
1149     ENDIF
1150
1151     ! Calcul de radliq (eau condensee pour le rayonnement)
1152     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1153     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1154     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1155     ! transmise au rayonnement;
1156     ! ----------------------------------------------------------------
1157     DO i = 1, klon
1158        IF (rneb(i,k).GT.0.0) THEN
1159           zneb(i) = MAX(rneb(i,k), seuil_neb)
1160     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
1161!           print*,zt(i),'fractionglace'
1162!>AJ
1163           radliq(i,k) = zoliq(i)/REAL(ninter+1)
1164        ENDIF
1165     ENDDO
1166     !
1167     DO n = 1, ninter
1168        DO i = 1, klon
1169           IF (rneb(i,k).GT.0.0) THEN
1170              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
1171              ! Initialization of zpluie and zice:
1172              zpluie=0
1173              zice=0
1174              IF (zneb(i).EQ.seuil_neb) THEN
1175                 ztot = 0.0
1176              ELSE
1177                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1178                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
1179                 if (ptconv(i,k)) then
1180                    zcl   =cld_lc_con
1181                    zct   =1./cld_tau_con
1182                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1183                         *fallvc(zrhol(i)) * zfice(i)
1184                 else
1185                    zcl   =cld_lc_lsc
1186                    zct   =1./cld_tau_lsc
1187                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1188                         *fallvs(zrhol(i)) * zfice(i)
1189                 endif
1190
1191                 ! si l'heterogeneite verticale est active, on utilise
1192                 ! la fraction volumique "vraie" plutot que la fraction
1193                 ! surfacique modifiee, qui est plus grande et reduit
1194                 ! sinon l'eau in-cloud de facon artificielle
1195                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1196                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1197                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
1198                 else
1199                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1200                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
1201                 endif
1202!AJ<
1203                 IF (.NOT. ice_thermo) THEN
1204                   ztot    = zchau + zfroi
1205                 ELSE
1206                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1207                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1208                   ztot    = zpluie    + zice
1209                 ENDIF
1210!>AJ
1211                 ztot    = MAX(ztot   ,0.0)
1212              ENDIF
1213              ztot    = MIN(ztot,zoliq(i))
1214!AJ<
1215     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
1216     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
1217!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1218!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1219!      si iflag_bergeron <> 2
1220!      A SUPPRIMER A TERME
1221              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
1222              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
1223              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
1224!>AJ
1225              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1226           ENDIF
1227        ENDDO  ! i = 1,klon
1228     ENDDO     ! n = 1,ninter
1229
1230     ! ----------------------------------------------------------------
1231     !
1232     IF (.NOT. ice_thermo) THEN
1233       DO i = 1, klon
1234         IF (rneb(i,k).GT.0.0) THEN
1235           d_ql(i,k) = zoliq(i)
1236           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1237                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1238         ENDIF
1239       ENDDO
1240     ELSE
1241!
1242!CR&JYG<
1243! On prend en compte l'effet Bergeron dans les flux de precipitation :
1244! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1245! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1246! et les precipitations est grossierement pris en compte en linearisant les equations
1247! et en approximant le processus de precipitation liquide par un processus a seuil.
1248! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1249! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1250! Le condensat precipitant solide est augmente.
1251! La vapeur d'eau est augmentee.
1252!
1253       IF ((iflag_bergeron .EQ. 2)) THEN
1254         DO i = 1, klon
1255           IF (rneb(i,k) .GT. 0.0) THEN
1256             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1257             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1258           if (fl_cor_ebil .GT. 0) then
1259             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1260             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1261!            Calcul de DT si toute les precips liquides congelent
1262             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1263!            On ne veut pas que T devienne superieur a la temp. de congelation.
1264!            donc que Delta > RTT-zt(i
1265             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1266             zt(i)      = zt(i)      + DeltaT
1267!            Eau vaporisee du fait de l'augmentation de T
1268             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1269!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1270             zq(i) = zq(i) +  Deltaq
1271!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1272             zcond(i)   = max( zcond(i)- Deltaq, 0. )
1273!            precip liquide qui congele ou qui s'evapore
1274             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1275             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1276!            bilan eau glacee
1277             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1278           else ! if (fl_cor_ebil .GT. 0)
1279!            ancien calcul
1280             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
1281             coef1 = RLMLT*zdqs(i)/RLVTT
1282             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1283             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1284             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
1285             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
1286             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
1287             zt(i)      = zt(i)      + DeltaT
1288           end if ! if (fl_cor_ebil .GT. 0)
1289           ENDIF  ! rneb(i,k) .GT. 0.0
1290         ENDDO
1291         DO i = 1, klon
1292           IF (rneb(i,k).GT.0.0) THEN
1293             d_ql(i,k) = (1-zfice(i))*zoliq(i)
1294             d_qi(i,k) = zfice(i)*zoliq(i)
1295             zrfl(i) = zrfl(i)+ zqprecl(i) &
1296                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1297             zifl(i) = zifl(i)+ zqpreci(i) &
1298                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1299           ENDIF                     
1300         ENDDO
1301!!
1302       ELSE  ! iflag_bergeron
1303!>CR&JYG
1304!!
1305       DO i = 1, klon
1306         IF (rneb(i,k).GT.0.0) THEN
1307!CR on prend en compte la phase glace
1308!JLD inutile car on ne passe jamais ici si .not.ice_thermo
1309!           if (.not.ice_thermo) then
1310!           d_ql(i,k) = zoliq(i)
1311!           d_qi(i,k) = 0.
1312!           else
1313           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1314           d_qi(i,k) = zfice(i)*zoliq(i)
1315!           endif
1316!AJ<
1317           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1318               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1319           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1320                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1321     !      zrfl(i) = zrfl(i)+  zpluie                         &
1322     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1323     !      zifl(i) = zifl(i)+  zice                    &
1324     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
1325
1326!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1327           IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
1328              zsolid = zrfl(i)
1329              zifl(i) = zifl(i)+zrfl(i)
1330              zrfl(i) = 0.
1331           if (fl_cor_ebil .GT. 0) then
1332              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1333                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
1334           else
1335              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1336                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
1337           end if
1338           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1339!RC   
1340
1341         ENDIF ! rneb(i,k).GT.0.0
1342       ENDDO
1343
1344       ENDIF  ! iflag_bergeron .EQ. 2
1345     ENDIF  ! .NOT. ice_thermo
1346
1347!CR: la fonte est faite au debut
1348!      IF (ice_thermo) THEN
1349!       DO i = 1, klon
1350!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1351!           zmelt = MIN(MAX(zmelt,0.),1.)
1352!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
1353!           zifl(i)=zifl(i)*(1.-zmelt)
1354!           print*,zt(i),'octavio1'
1355!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1356!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1357!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1358!       ENDDO
1359!     ENDIF
1360
1361       
1362     IF (.NOT. ice_thermo) THEN
1363       DO i = 1, klon
1364         IF (zt(i).LT.RTT) THEN
1365           psfl(i,k)=zrfl(i)
1366         ELSE
1367           prfl(i,k)=zrfl(i)
1368         ENDIF
1369       ENDDO
1370     ELSE
1371     ! JAM*************************************************
1372     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1373     ! *****************************************************
1374
1375       DO i = 1, klon
1376     !   IF (zt(i).LT.RTT) THEN
1377           psfl(i,k)=zifl(i)
1378     !   ELSE
1379           prfl(i,k)=zrfl(i)
1380     !   ENDIF
1381!>AJ
1382       ENDDO
1383     ENDIF
1384     ! ----------------------------------------------------------------
1385     ! Fin de formation des precipitations
1386     ! ----------------------------------------------------------------
1387     !
1388     ! Calculer les tendances de q et de t:
1389     !
1390     DO i = 1, klon
1391        d_q(i,k) = zq(i) - q(i,k)
1392        d_t(i,k) = zt(i) - t(i,k)
1393     ENDDO
1394     !
1395     !AA--------------- Calcul du lessivage stratiforme  -------------
1396
1397     DO i = 1,klon
1398        !
1399        if(zcond(i).gt.zoliq(i)+1.e-10) then
1400         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1401        else
1402         beta(i,k) = 0.
1403        endif
1404        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1405             * (paprs(i,k)-paprs(i,k+1))/RG
1406        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1407           !AA lessivage nucleation LMD5 dans la couche elle-meme
1408          IF (iflag_t_glace.EQ.0) THEN
1409           if (t(i,k) .GE. t_glace_min_old) THEN
1410              zalpha_tr = a_tr_sca(3)
1411           else
1412              zalpha_tr = a_tr_sca(4)
1413           endif
1414          ELSE ! of IF (iflag_t_glace.EQ.0)
1415           if (t(i,k) .GE. t_glace_min) THEN
1416              zalpha_tr = a_tr_sca(3)
1417           else
1418              zalpha_tr = a_tr_sca(4)
1419           endif
1420          ENDIF
1421           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1422           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1423           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1424           !
1425           ! nucleation avec un facteur -1 au lieu de -0.5
1426           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1427           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1428        ENDIF
1429        !
1430     ENDDO      ! boucle sur i
1431     !
1432     !AA Lessivage par impaction dans les couches en-dessous
1433     DO kk = k-1, 1, -1
1434        DO i = 1, klon
1435           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1436             IF (iflag_t_glace.EQ.0) THEN
1437              if (t(i,kk) .GE. t_glace_min_old) THEN
1438                 zalpha_tr = a_tr_sca(1)
1439              else
1440                 zalpha_tr = a_tr_sca(2)
1441              endif
1442             ELSE ! of IF (iflag_t_glace.EQ.0)
1443              if (t(i,kk) .GE. t_glace_min) THEN
1444                 zalpha_tr = a_tr_sca(1)
1445              else
1446                 zalpha_tr = a_tr_sca(2)
1447              endif
1448             ENDIF
1449              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1450              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1451              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1452           ENDIF
1453        ENDDO
1454     ENDDO
1455     !
1456     !AA===============================================================
1457     !                     FIN DE LA BOUCLE VERTICALE 
1458  end DO
1459  !
1460  !AA==================================================================
1461  !
1462  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1463  !
1464!CR: si la thermo de la glace est active, on calcule zifl directement
1465  IF (.NOT.ice_thermo) THEN
1466  DO i = 1, klon
1467     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1468!AJ<
1469!        snow(i) = zrfl(i)
1470        snow(i) = zrfl(i)+zifl(i)
1471!>AJ
1472        zlh_solid(i) = RLSTT-RLVTT
1473     ELSE
1474        rain(i) = zrfl(i)
1475        zlh_solid(i) = 0.
1476     ENDIF
1477  ENDDO
1478
1479  ELSE
1480     DO i = 1, klon
1481        snow(i) = zifl(i)
1482        rain(i) = zrfl(i)
1483     ENDDO
1484   
1485   ENDIF
1486  !
1487  ! For energy conservation : when snow is present, the solification
1488  ! latent heat is considered.
1489!CR: si thermo de la glace, neige deja prise en compte
1490  IF (.not.ice_thermo) THEN
1491  DO k = 1, klev
1492     DO i = 1, klon
1493        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
1494        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
1495        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1496        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
1497     END DO
1498  END DO
1499  ENDIF
1500  !
1501
1502  if (ncoreczq>0) then
1503     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1504  endif
1505
1506END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.