Index: LMDZ6/trunk/libf/phylmdiso/fisrtilp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/fisrtilp.F90	(revision 4653)
+++ 	(revision )
@@ -1,2742 +1,0 @@
-!
-! $Id: fisrtilp.F90 3493 2019-05-07 08:45:07Z idelkadi $
-!
-!
-SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
-     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
-     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
-     frac_impa, frac_nucl, beta,                        &
-     prfl, psfl, rhcl, zqta, fraca,                     &
-     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
-     iflag_ice_thermo &
-#ifdef ISO
-     &           ,xt,xtrain,xtsnow                      &
-     &           ,d_xt,d_xtl,d_xti,pxtrainfl,pxtsnowfl        &
-#endif      
-     &  )
-
-  !
-  USE dimphy
-  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
-  USE print_control_mod, ONLY: prt_level, lunout
-  USE lmdz_cloudth, only : cloudth, cloudth_v3, cloudth_v6
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
-  USE phys_local_var_mod, ONLY: rneblsvol
-  ! flag to include modifications to ensure energy conservation (if flag >0)
-  USE add_phys_tend_mod, only : fl_cor_ebil 
-  USE lscp_ini_mod, ONLY: iflag_t_glace,t_glace_min, t_glace_max, exposant_glace
-  USE lscp_ini_mod, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol 
-  USE lscp_ini_mod, ONLY: coef_eva, coef_eva_i, ffallv_lsc, ffallv_con
-  USE lscp_ini_mod, ONLY: cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
-  USE lscp_ini_mod, ONLY: reevap_ice, iflag_bergeron, iflag_fisrtilp_qsat, iflag_pdf 
-      use phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv
-      use phys_output_var_mod, ONLY : cloudth_sigmath,cloudth_sigmaenv
-
-
-
-
-#ifdef ISO
-  USE infotrac_phy, ONLY: ntraciso=>ntiso,niso,itZonIso
-  USE isotopes_mod
-!, ONLY: essai_convergence,bidouille_anti_divergence, &
-!        &       Rdefault,ridicule,pxtice,pxtmelt,iso_eau,iso_HDO, &
-!        &       iso_O18,ridicule_rain
-  USE isotopes_routines_mod, ONLY: iso_revap_fisrtilp, &
-        condiso_liq_ice_vectall
-#ifdef ISOVERIF
-  USE isotopes_verif_mod, ONLY: errmax,errmaxrel,deltalim_snow,deltalim, &
-        iso_verif_egalite_choix,iso_verif_aberrant_choix, &
-        iso_verif_positif,iso_verif_positif_choix,iso_verif_noNaN, &
-        iso_verif_aberrant,iso_verif_positif_choix_nostop, &
-        iso_verif_aberrant_encadre,iso_verif_egalite_choix_nostop, &
-        iso_verif_aberrant_nostop,iso_verif_o18_aberrant_nostop, &
-        deltaD,deltaO,iso_verif_egalite
-#endif
-#ifdef ISOTRAC
-    use isotrac_mod, only: option_traceurs,option_tmin,seuil_tag_tmin_ls, &
-&       izone_poubelle, option_cond,izone_cond,index_iso, index_zone,izone_oce,izone_ddft
-    use isotrac_routines_mod, only: iso_recolorise_condensation
-    use isotopes_routines_mod, only: condiso_liq_ice_vectall_trac
-#ifdef ISOVERIF
-  USE isotopes_verif_mod, ONLY: iso_verif_traceur_nostop,iso_verif_traceur
-  USE isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille
-#endif
-#endif
-#endif
-
-  IMPLICIT none
-  !======================================================================
-  ! Auteur(s): Z.X. Li (LMD/CNRS)
-  ! Date: le 20 mars 1995
-  ! Objet: condensation et precipitation stratiforme.
-  !        schema de nuage
-  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
-  !             et ilp (il pleut, L. Li)
-  ! Principales parties:
-  ! P0> Thermalisation des precipitations venant de la couche du dessus
-  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
-  ! P2> Formation du nuage (en k) 
-  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
-  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
-  ! les valeurs de T et Q initiales
-  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
-  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
-  ! P2.B> Nuage "tout ou rien"
-  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
-  ! P3> Formation de la precipitation (en k)
-  !======================================================================
-  ! JLD: 
-  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
-  ! * fl_cor_ebil doit etre > 0 ; 
-  !   fl_cor_ebil= 0 pour reproduire anciens bugs
-  !======================================================================
-  include "YOMCST.h"
-  !
-  ! Principaux inputs:
-  !
-  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
-  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
-  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
-  INTEGER,                         INTENT(IN)   :: iflag_cld_th
-  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
-  !
-  ! Inputs lies aux thermiques
-  ! 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zthl
-  !
-  !  Input/output
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
-  !
-  ! Principaux outputs:
-  !
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb ! fraction nuageuse
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
-  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
-  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
-  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
-  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl 
-
-#ifdef ISO
-      real double_to_real
-      double precision real_to_double
-      !integer niso
-      real xt(ntraciso,klon,klev)
-      real d_xt(ntraciso,klon,klev)
-      real d_xtl(ntraciso,klon,klev)
-      real d_xti(ntraciso,klon,klev)
-      real xtrain(ntraciso,klon)
-      real xtsnow(ntraciso,klon)
-      real pxtrainfl(ntraciso,klon,klev+1)
-      real pxtsnowfl(ntraciso,klon,klev+1)      
-      ! xt <-> q
-      ! xtrain et xtsonw <-> rain et snow
-      ! d_xt et d_xtl <-> d_q et d_ql
-      ! pxtrainfl et pxtsnowflux <-> prfl et psfl
-#endif
-
-
-  !AA
-  ! Coeffients de fraction lessivee : pour OFF-LINE
-  !
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_nucl
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_1nucl
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_impa
-  !
-  ! Fraction d'aerosols lessivee par impaction et par nucleation
-  ! POur ON-LINE
-  !
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
-  !AA
-  ! --------------------------------------------------------------------------------
-  !
-  ! Options du programme:
-  !
-  REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
-  !$OMP THREADPRIVATE(seuil_neb)
-
-
-  INTEGER ninter ! sous-intervals pour la precipitation
-  PARAMETER (ninter=5)
-  INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
-  !$OMP THREADPRIVATE(iflag_evap_prec)
-  !
-  LOGICAL cpartiel ! condensation partielle
-  PARAMETER (cpartiel=.TRUE.)
-  REAL t_coup
-  PARAMETER (t_coup=234.0)
-  REAL DDT0
-  PARAMETER (DDT0=.01)
-  REAL ztfondue
-  PARAMETER (ztfondue=278.15)
-  ! --------------------------------------------------------------------------------
-  !
-  ! Variables locales:
-  !
-  INTEGER i, k, n, kk
-  INTEGER,save::itap=0
-  !$OMP THREADPRIVATE(itap)
-
-  REAL qsl, qsi
-  real zct      ,zcl
-  INTEGER ncoreczq  
-  REAL ctot(klon,klev)
-  REAL ctot_vol(klon,klev)
-  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
-  REAL zdqsdT_raw(klon)
-  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
-
-  logical lognormale(klon)
-  logical ice_thermo
-  LOGICAL convergence(klon)
-  INTEGER n_i(klon), iter
-  REAL cste
-
-  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
-  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
-  real erf   
-  REAL qcloud(klon)
- 
-  REAL zrfl(klon), zrfln(klon), zqev, zqevt 
-  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
-  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq 
-  REAL zoliqp(klon), zoliqi(klon)
-  REAL zt(klon)
-! JBM (3/14) nexpo is replaced by exposant_glace
-! REAL nexpo ! exponentiel pour glace/eau
-! INTEGER, PARAMETER :: nexpo=6
-  INTEGER exposant_glace_old
-  REAL t_glace_min_old
-  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
-  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
-  REAL zmelt, zpluie, zice
-  REAL dzfice(klon)
-  REAL zsolid
-!!!!
-!  Variables pour Bergeron
-  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
-  REAL zqpreci(klon), zqprecl(klon)
-! Variable pour conservation enegie des precipitations
-  REAL zmqc(klon)
-
-#ifdef ISO
-      real zxtrfl(ntraciso,klon)
-      real zxtrfln(ntraciso,klon)
-      REAL zxtifl(ntraciso,klon), zxtifln(ntraciso,klon)  
-      real zxt(ntraciso,klon) ! equivalent local de xt
-      real zxtn(ntraciso,klon) ! isotpes nuageux
-      real zxtn_reduit(niso)
-      real zxtcond(ntraciso,klon)
-      real zxtoliq(ntraciso,klon)
-      real zq_ancien(klon)
-      REAL zxtpreci(ntraciso,klon), zxtprecl(ntraciso,klon)
-      integer ixt
-#ifdef ISOVERIF   
-!      integer iso_verif_aberrant_nostop
-!      integer iso_verif_egalite_choix_nostop
-!      integer iso_verif_positif_nostop
-!      integer iso_verif_positif_choix_nostop
-      integer trace_cas(klon)
-!      integer iso_verif_traceur_nostop
-#endif      
-      integer iso_verif_nonan_nostop
-      real zxtice(ntraciso,klon),zxtliq(ntraciso,klon)
-      real zrfl_ancien(klon)
-      real zqev_diag(klon)
-      real zxtrfl_ancien(ntraciso,klon)
-      real zxt_ancien(ntraciso,klon)
-#ifdef ISOTRAC
-        real zxtres(ntraciso) ! humidite spec nuageuse après condensation et precip
-        real zqcs(klon) ! humidite spec en ciel clair
-        real zxtcs(ntraciso,klon)        
-        real zcondn(klon) ! humidite spec de l'eau liquide dans le nuage
-        real zxtcondn(ntraciso,klon)
-        integer izone
-#ifdef ISOVERIF           
-        real zq_verif
-#endif        
-#endif 
-!      logical negation ! pour nier conditions logiques
-#endif
-
-  !
-  LOGICAL appel1er
-  SAVE appel1er
-  !$OMP THREADPRIVATE(appel1er)
-  !
-! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
-! iflag_oldbug_fisrtilp=1 ajoute le BUG
-  INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
-  !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
-  !---------------------------------------------------------------
-  !
-  !AA Variables traceurs:
-  !AA  Provisoire !!! Parametres alpha du lessivage
-  !AA  A priori on a 4 scavenging # possibles
-  !
-  REAL a_tr_sca(4)
-  save a_tr_sca
-  !$OMP THREADPRIVATE(a_tr_sca)
-  !
-  ! Variables intermediaires
-  !
-  REAL zalpha_tr
-  REAL zfrac_lessi
-  REAL zprec_cond(klon)
-  !AA
-! RomP >>> 15 nov 2012
-  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
-! RomP <<<
-  REAL zmair(klon), zcpair, zcpeau
-  !     Pour la conversion eau-neige
-  REAL zlh_solid(klon), zm_solid
-  !---------------------------------------------------------------
-  !
-  ! Fonctions en ligne:
-  !
-  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
-                     ! (Heymsfield & Donner, 1990)
-  REAL zzz
-
-  include "YOETHF.h"
-  include "FCTTRE.h"
-  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
-  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
-  !
-  DATA appel1er /.TRUE./
-  !ym
-!CR: pour iflag_ice_thermo=2, on active que la convection
-!  ice_thermo = iflag_ice_thermo .GE. 1
-
-  itap=itap+1
-  znebprecip(:)=0.
-
-  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
-  zdelq=0.0
-  ctot_vol(1:klon,1:klev)=0.0
-  rneblsvol(1:klon,1:klev)=0.0
-
-  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
-  IF (appel1er) THEN
-     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
-     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
-     CALL getin_p('seuil_neb',seuil_neb)
-     write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp 
-     !
-     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
-     WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
-     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
-     
-     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
-        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
-        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
-        !         CALL abort
-     ENDIF
-     appel1er = .FALSE.
-
-#ifdef ISO
-        ! cam verif 
-#ifdef ISOVERIF
-        write(*,*) 'ilp 310'
-          do k=1,klev
-           do i=1,klon             
-            if (iso_eau.gt.0) then        
-             call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &   
-     &                   'il pleut 205',errmax,errmaxrel)
-            endif !if (iso_eau.gt.0) then
-            if (iso_HDO.gt.0) then    
-             if (q(i,k).gt.ridicule) then
-               call iso_verif_aberrant(xt(iso_HDO,i,k)/q(i,k), &
-     &           'il pleut 210')
-             endif !if (zq(i).gt.ridicule) then
-            endif  !if (iso_HDO.gt.0) then
-           enddo !do k=1,klev
-          enddo !do i=1,klon
-        do k=1,klev
-           do i=1,klon
-             call iso_verif_positif(370.0-t(i,k),'il pleut 250')
-             call iso_verif_positif(t(i,k)-100.0,'il pleut 251')
-           enddo
-        enddo
-        write(*,*) 'ilp 331'
-#endif
-         ! end cam verif
-#endif 
-
-     !
-     !AA initialiation provisoire
-     a_tr_sca(1) = -0.5
-     a_tr_sca(2) = -0.5
-     a_tr_sca(3) = -0.5
-     a_tr_sca(4) = -0.5
-     !
-     !AA Initialisation a 1 des coefs des fractions lessivees 
-     !
-     !cdir collapse
-     DO k = 1, klev
-        DO i = 1, klon
-           pfrac_nucl(i,k)=1.
-           pfrac_1nucl(i,k)=1.
-           pfrac_impa(i,k)=1.
-           beta(i,k)=0.  !RomP initialisation
-        ENDDO
-     ENDDO
-
-  ENDIF          !  test sur appel1er
-  !
-  !MAf Initialisation a 0 de zoliq
-  !      DO i = 1, klon
-  !         zoliq(i)=0.
-  !      ENDDO 
-  ! Determiner les nuages froids par leur temperature
-  !  nexpo regle la raideur de la transition eau liquide / eau glace.
-  !
-!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
-  IF (iflag_t_glace.EQ.0) THEN
-!   ztglace = RTT - 15.0
-    t_glace_min_old = RTT - 15.0
-    !AJ<
-    IF (ice_thermo) THEN
-!     nexpo = 2
-      exposant_glace_old = 2
-    ELSE
-!     nexpo = 6
-      exposant_glace_old = 6
-    ENDIF
-    
-  ENDIF
-  
-!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
-!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
-!>AJ
-  !cc      nexpo = 1
-  !
-  ! Initialiser les sorties:
-  !
-  !cdir collapse
-  DO k = 1, klev+1
-     DO i = 1, klon
-        prfl(i,k) = 0.0
-        psfl(i,k) = 0.0
-#ifdef ISO                 
-          do ixt=1,ntraciso
-            pxtrainfl(ixt,i,k)=0.0
-            pxtsnowfl(ixt,i,k)=0.0
-          enddo !do ixt=1,niso
-#endif
-     ENDDO
-  ENDDO
-
-  !cdir collapse
-  DO k = 1, klev
-     DO i = 1, klon
-        d_t(i,k) = 0.0
-        d_q(i,k) = 0.0
-        d_ql(i,k) = 0.0
-        d_qi(i,k) = 0.0
-#ifdef ISO
-          do ixt=1,ntraciso
-            d_xt(ixt,i,k) = 0.0
-            d_xtl(ixt,i,k) = 0.0
-            d_xti(ixt,i,k) = 0.0
-          enddo !do ixt=1,niso
-#endif 
-        rneb(i,k) = 0.0
-        radliq(i,k) = 0.0
-        frac_nucl(i,k) = 1. 
-        frac_impa(i,k) = 1. 
-     ENDDO
-  ENDDO
-  DO i = 1, klon
-     rain(i) = 0.0
-     snow(i) = 0.0
-     zoliq(i)=0.
-#ifdef ISO
-         do ixt=1,ntraciso
-            xtrain(ixt,i) = 0.0
-            xtsnow(ixt,i) = 0.0
-            zxtoliq(ixt,i)=0.
-         enddo !do ixt=1,niso
-#endif
-     !     ENDDO
-     !
-     ! Initialiser le flux de precipitation a zero
-     !
-     !     DO i = 1, klon
-     zrfl(i) = 0.0
-     zifl(i) = 0.0
-     zrfln(i) = 0.0 ! C Risi pour plus de securite
-     zifln(i) = 0.0 ! C Risi pour plus de securite
-#ifdef ISO
-          do ixt=1,ntraciso
-            zxtrfl(ixt,i)=0.0
-            zxtifl(ixt,i)=0.0
-            zxtrfln(ixt,i)=0.0
-            zxtifln(ixt,i)=0.0
-          enddo
-#endif
-     zneb(i) = seuil_neb
-  ENDDO
-  !
-  !
-  !AA Pour plus de securite 
-
-  zalpha_tr   = 0.
-  zfrac_lessi = 0.
-
-  !AA==================================================================
-  !
-  ncoreczq=0
-  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
-  !
-  DO k = klev, 1, -1
-     !
-     !
-#ifdef ISO      
-      ! cam debug!      
-       ! verif en debut de boucle
-#ifdef ISOVERIF
-       !write(*,*) 'ilp 478: boucle en k: k=',k         
-       do i=1,klon
-        if (iso_eau.gt.0) then
-        call iso_verif_egalite_choix(zxtrfl(iso_eau,i), &
-     &           zrfl(i), &
-     &           'il pleut 310',errmax,errmaxrel)
-        endif !if (iso_eau.gt.0) then
-        if (iso_HDO.gt.0) then  
-             call iso_verif_aberrant_choix(zxtrfl(iso_HDO,i) &
-     &           ,zrfl(i),ridicule_rain,deltalim_snow,'il pleut 316')
-        endif  !if (iso_HDO.gt.0) then  
-#ifdef ISOTRAC   
-        call iso_verif_traceur(zxtrfl(1,i), &
-     &          'il pleut 365')   
-        call iso_verif_traceur_pbidouille(zxtrfl(1,i), &
-     &          'il pleut 388')
-#endif         
-       enddo !do i=1,klon
-#endif
-#ifdef ISOVERIF   
-        !write(*,*) 'ilp tmp 412'      
-       do i=1,klon
-         do ixt=1,ntraciso
-           call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 406')
-           call iso_verif_noNaN(zxtoliq(ixt,i),'ilp 407')
-         enddo !do ixt=1,ntraciso
-       enddo !do i=1,klon
-#endif
-      ! end verif en debut de boucle
-      ! end cam debug
-#endif
-     !AA===============================================================
-     !
-     ! Initialisation temperature et vapeur
-     DO i = 1, klon
-        zt(i)=t(i,k)
-        zq(i)=q(i,k)
-#ifdef ISOVERIF
-        if (q(i,k).lt.0.0) then
-               write(*,*) 'ilp 400: avant clmain 31janv: i,k,q=', &
-     &           i,k,q(i,k)
-               CALL abort_physic('ilp 484', 'on stop', 1)
-        endif 
-#endif
-#ifdef ISO
-           zq(i)=max(0.0,zq(i))
-#ifdef ISOTRAC
-#ifdef ISOVERIF         
-          call iso_verif_traceur(xt(1,i,k),'ilp 404')           
-          call iso_verif_traceur_pbidouille(xt(1,i,k),'ilp 407')
-#endif          
-#endif          
-         do ixt=1,ntraciso
-           zxt(ixt,i)=xt(ixt,i,k)
-           zxt(ixt,i)=max(0.0,zxt(ixt,i))
-         enddo !do ixt=1,niso  
-
-         ! cam verif
-#ifdef ISOVERIF
-         !write(*,*) 'ilp tmp 562'
-             call iso_verif_positif(370.0-zt(i),'il pleut 425')
-             call iso_verif_positif(zt(i)-100.0,'il pleut 426')
-          if (iso_eau.gt.0) then        
-            call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
-     &                   'il pleut 289',errmax,errmaxrel)
-          endif !if (iso_eau.gt.0) then
-          if (iso_HDO.gt.0) then    
-             if (zq(i).gt.ridicule) then
-               call iso_verif_aberrant(zxt(iso_HDO,i)/zq(i), &
-     &           'il pleut 337')
-             endif !if (zq(i).gt.ridicule) then
-          endif  !if (iso_HDO.gt.0) then  
-#ifdef ISOTRAC
-          call iso_verif_traceur(zxt(1,i),'ilp 431')
-#endif           
-#endif
-         ! end cam verif
-#endif  
-
-     ENDDO
-     !
-     ! ----------------------------------------------------------------
-     ! P0> Thermalisation des precipitations venant de la couche du dessus
-     ! ----------------------------------------------------------------
-     ! Calculer la varition de temp. de l'air du a la chaleur sensible
-     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
-     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
-     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
-     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
-     ! de l'enthalpie  de la couche avec la temperature
-     ! Variables calculees ou modifiees:
-     !   -  zt: temperature de la cocuhe
-     !   - zmqc: masse de precip qui doit etre thermalisee
-     !
-     IF(k.LE.klevm1) THEN         
-        DO i = 1, klon
-           !IM
-           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
-           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
-           zcpair=RCPD*(1.0+RVTMP2*zq(i))
-           zcpeau=RCPD*RVTMP2
-         if (fl_cor_ebil .GT. 0) then
-           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm 
-           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la 
-           ! derniere couche
-           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
-           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
-           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
-                 / (zcpair + zmqc(i)*zcpeau)
-         else ! si on maintient les anciennes erreurs
-           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
-                + zmair(i)*zcpair*zt(i) ) &
-                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
-         end if
-        ENDDO
-     ELSE  ! IF(k.LE.klevm1) 
-        DO i = 1, klon 
-           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 
-           zmqc(i) = 0.
-        ENDDO
-     ENDIF ! end IF(k.LE.klevm1)
-     !
-     ! ----------------------------------------------------------------
-     ! P1> Calcul de l'evaporation de la precipitation
-     ! ----------------------------------------------------------------
-     ! On evapore une partie des precipitations venant de la maille du dessus.
-     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
-     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
-     ! Variables calculees ou modifiees:
-     !   - zrfl et zifl: flux de precip liquide et glace
-     !   - zq, zt: humidite et temperature de la cocuhe
-     !   - zmqc: masse de precip qui doit etre thermalisee
-     !
-#ifdef ISO
-      DO i = 1, klon
-         zq_ancien(i)=zq(i)
-         zqev_diag(i)=0.0
-         zrfl_ancien(i)=zrfl(i)
-         do ixt=1,ntraciso
-           zxtrfl_ancien(ixt,i)=zxtrfl(ixt,i)
-           zxt_ancien(ixt,i)=zxt(ixt,i)
-         enddo 
-      enddo   ! DO i = 1, klon  
-#ifdef ISOVERIF
-        !write(*,*) 'ilp tmp 616'
-      if (iso_eau.gt.0) then
-       DO i = 1, klon           
-          call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
-     &           'il pleut 426a',errmax,errmaxrel) 
-          call iso_verif_egalite_choix(zxtrfl_ancien(iso_eau,i), &
-     &            zrfl_ancien(i), &
-     &           'il pleut 426b',errmax,errmaxrel)
-        enddo
-      endif
-#ifdef ISOTRAC
-       do i=1,klon      
-            call iso_verif_traceur(zxt_ancien(1,i),'ilp 532')
-            call iso_verif_traceur(zxtrfl_ancien(1,i), &
-     &           'ilp 533')  
-       enddo !do i=1,klon  
-#endif      
-#endif      
-#endif
-
-
-     IF (iflag_evap_prec>=1) THEN
-        DO i = 1, klon
-!          S'il y a des precipitations
-           IF (zrfl(i)+zifl(i).GT.0.) THEN
-              ! Calcul du qsat
-              IF (thermcep) THEN
-                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
-                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-                 zqs(i)=MIN(0.5,zqs(i))
-                 zcor=1./(1.-RETV*zqs(i))
-                 zqs(i)=zqs(i)*zcor
-              ELSE
-                 IF (zt(i) .LT. t_coup) THEN
-                    zqs(i) = qsats(zt(i)) / pplay(i,k)
-                 ELSE
-                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
-                 ENDIF
-              ENDIF
-           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
-        ENDDO
-!AJ<
-
-       IF (.NOT. ice_thermo) THEN
-        DO i = 1, klon
-!          S'il y a des precipitations
-           IF (zrfl(i)+zifl(i).GT.0.) THEN                
-                ! Evap max pour ne pas saturer la fraction sous le nuage
-                ! Evap max jusqu'à atteindre la saturation dans la partie
-                ! de la maille qui est sous le nuage de la couche du dessus
-                !!! On ne tient compte de cette fraction que sous une seule
-                !!! couche sous le nuage
-                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
-             ! Ajout de la prise en compte des precip a thermiser
-             ! avec petite reecriture 
-             if  (fl_cor_ebil .GT. 0) then ! nouveau
-                ! Calcul de l'evaporation du flux de precip herite
-                !   d'au-dessus
-                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
-                     * zmair(i)/pplay(i,k)*zt(i)*RD
-                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
-                ! Seuil pour ne pas saturer la fraction sous le nuage
-                zqev = MIN (zqev, zqevt)
-                ! Nouveau flux de precip
-                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
-                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
-                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
-                  zrfln(i)=0.
-                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
-                END IF
-                ! Nouvelle vapeur 
-                zq(i) = zq(i) + zqev
-                zmqc(i) = zmqc(i)-zqev
-                ! Nouvelle temperature (chaleur latente)
-                zt(i) = zt(i) - zqev &
-                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-!!JLD debut de partie a supprimer a terme
-            else ! if  (fl_cor_ebil .GT. 0)
-                ! Calcul de l'evaporation du flux de precip herite
-                !   d'au-dessus
-                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
-                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
-                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
-                ! Seuil pour ne pas saturer la fraction sous le nuage
-                zqev = MIN (zqev, zqevt)
-                ! Nouveau flux de precip
-                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
-                     /RG/dtime
-                ! Aucun flux liquide pour T < t_coup
-                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-                ! Nouvelle vapeur 
-                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
-                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-                ! Nouvelle temperature (chaleur latente)
-                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
-                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
-                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-              end if ! if  (fl_cor_ebil .GT. 0)
-
-!!JLD fin de partie a supprimer a terme
-                zrfl(i) = zrfln(i)
-                zifl(i) = 0.
-#ifdef ISO
-         zqev_diag(i)=zqev
-#endif
-
-           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
-
-
-        ENDDO
-#ifdef ISO
-
-      ! * traitement de l'évaporation
-      call iso_revap_fisrtilp(klon,klev,k, &
-     &            zrfl_ancien,zrfl,zrfln,zt,zxt_ancien, &
-     &            zxtrfl,zxtrfl_ancien,zxtrfln,zxt, &
-     &            paprs,dtime, &
-     &            zqs,zq_ancien,zqev_diag,zq)  
-#endif
-
-!
-       ELSE ! (.NOT. ice_thermo) 
-
-#ifdef ISO
-        CALL abort_physic('ilp 664', 'isotopes pas encore prevus dans ice_thermo', 1)        
-#endif
-!      ================================
-!      Avec thermodynamique de la glace
-!      ================================
-        DO i = 1, klon
-!AJ<
-!        S'il y a des precipitations
-         IF (zrfl(i)+zifl(i).GT.0.) THEN
-
-         IF (iflag_evap_prec==1) THEN
-            znebprecip(i)=zneb(i)
-         ELSE
-            znebprecip(i)=MAX(zneb(i),znebprecip(i))
-         ENDIF
-     
-        ! Evap max pour ne pas saturer la fraction sous le nuage
-         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
-
-         !JAM
-         ! On differencie qsat pour l'eau et la glace
-         ! Si zdelta=1. --> glace
-         ! Si zdelta=0. --> eau liquide
-        
-         ! Calcul du qsat par rapport a l'eau liquide
-         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
-         qsl= MIN(0.5,qsl)
-         zcor= 1./(1.-RETV*qsl)
-         qsl= qsl*zcor
-         
-         ! Calcul de l'evaporation du flux de precip venant du dessus
-         ! Formulation en racine du flux de precip
-         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
-         IF (iflag_evap_prec==3) THEN
-         zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
-              *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
-              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         ELSE
-         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
-              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         ENDIF
-
-
-         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
-              *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
-         
-         ! Calcul du qsat par rapport a la glace
-         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
-         qsi= MIN(0.5,qsi)
-         zcor= 1./(1.-RETV*qsi)
-         qsi= qsi*zcor
-
-         ! Calcul de la sublimation du flux de precip solide herite
-         !   d'au-dessus
-         IF (iflag_evap_prec==3) THEN
-         zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
-              *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
-              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         ELSE
-         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
-              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         ENDIF
-         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
-              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
-
-        !JAM
-        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas 
-        ! la fraction de la couche sous le nuage sinon on repartit zqev0 
-        ! en conservant la proportion liquide / glace
-     
-         IF (zqevt+zqevti.GT.zqev0) THEN
-            zqev=zqev0*zqevt/(zqevt+zqevti)
-            zqevi=zqev0*zqevti/(zqevt+zqevti)
-         ELSE
-!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
-!       liquides et solides meme si on ne sature pas la couche. 
-!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
-!            zqev=zqevt
-!            zqevi=zqevti
-             IF (zqevt+zqevti.GT.0.) THEN
-                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
-                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
-             ELSE
-             zqev=0.
-             zqevi=0.
-             ENDIF
-         ENDIF
-         ! Nouveaux flux de precip liquide et solide 
-         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
-                                 /RG/dtime)
-         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
-                                 /RG/dtime)
-
-         ! Mise a jour de la vapeur, temperature et flux de precip
-         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
-         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
-                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
-                  + (zifln(i)-zifl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
-                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-       else ! sans correction thermalisation des precips
-         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
-                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
-                  + (zifln(i)-zifl(i)) &
-                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
-                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
-       end if
-        ! Nouvelles vaeleurs des precips liquides et solides 
-         zrfl(i) = zrfln(i)
-         zifl(i) = zifln(i)
-!        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
-
-!CR ATTENTION: deplacement de la fonte de la glace
-!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
-!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
-!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
-           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
-           ! fraction de la precip solide qui est fondue
-           zmelt = MIN(MAX(zmelt,0.),1.)
-           ! Fusion de la glace
-           zrfl(i)=zrfl(i)+zmelt*zifl(i)
-           if (fl_cor_ebil .LE. 0) then
-             ! the following line should not be here. Indeed, if zifl is modified
-             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
-             ! and therefore the change in temperature computed below is wrong
-             zifl(i)=zifl(i)*(1.-zmelt)
-           end if
-           ! Chaleur latente de fusion
-        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
-           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-        else ! sans correction thermalisation des precips
-           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
-        end if
-           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
-             zifl(i)=zifl(i)*(1.-zmelt)
-           end if
-
-           ELSE
-              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
-              ! sous nuageuse utilisee pour la pluie.
-              znebprecip(i)=0.
-           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
-        ENDDO
-
-      ENDIF ! (.NOT. ice_thermo) 
-     
-     ! ----------------------------------------------------------------
-     ! Fin evaporation de la precipitation
-     ! ----------------------------------------------------------------
-     ENDIF ! (iflag_evap_prec>=1)
-     !
-     ! Calculer Qs et L/Cp*dQs/dT:
-#ifdef ISOVERIF
-!        write(*,*) 'ilp tmp 849'
-           do i=1,klon
-             call iso_verif_positif(370.0-zt(i),'il pleut 938')
-             call iso_verif_positif(zt(i)-100.0,'il pleut 939')
-             do ixt=1,ntraciso
-              call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 614')
-             enddo   
-             if (zq(i).gt.ridicule) then
-               if (iso_HDO.gt.0) then
-                  call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
-     &                   'il pleut 856')                 
-                  if (iso_O18.gt.0) then                
-                   if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
-     &               zxt(iso_O18,i)/zq(i),'il pleut 859').eq.1) then
-                     write(*,*) 'i,k,zq=',i,k,zq(i)
-                     write(*,*) 'zqev_diag(i)=',zqev_diag(i)
-                     write(*,*) 'zrfl_ancien(i)=',zrfl_ancien(i)
-                     write(*,*) 'zrfln(i)=',zrfln(i)
-                     write(*,*) 'zq_ancien(i)=',zq_ancien(i)
-                     write(*,*) 'deltaD(zrfl_ancien)=',deltaD(zxtrfl_ancien(iso_HDO,i))
-                     write(*,*) 'deltaO(zrfl_ancien)=',deltaO(zxtrfl_ancien(iso_O18,i))
-                     write(*,*) 'deltaD(zq_ancien)=',deltaD(zxt_ancien(iso_HDO,i))
-                     write(*,*) 'deltaO(zq_ancien)=',deltaO(zxt_ancien(iso_O18,i))
-                     write(*,*) 'zt=',zt(i)
-                     stop
-                endif !if (iso_verif_O18_aberrant_nostop
-               endif !if (iso_O18.gt.0) then
-              endif !if (iso_HDO.gt.0) then
-             endif !if (zq(i).gt.ridicule) then                
-           enddo !do i=1,klon
-#ifdef ISOTRAC
-!      write(*,*) 'ilp tmp 633'
-           do i=1,klon
-              call iso_verif_traceur(zxtrfl(1,i),'il pleut 632')
-           enddo
-!      write(*,*) 'ilp tmp 637,thermcep'
-#endif           
-#endif    
-#ifdef ISOVERIF
-!        write(*,*) 'ilp tmp 888'
-           do i=1,klon
-            do ixt=1,ntraciso
-             call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 624')
-            enddo !do ixt=1,ntraciso
-           enddo ! do i=1,klon
-!        write(*,*) 'ilp tmp 638: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
-#endif
-     !
-     IF (thermcep) THEN
-        DO i = 1, klon
-           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
-           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-       if  (fl_cor_ebil .GT. 0) then ! nouveau
-           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-       else    
-           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
-       endif
-           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-           zqs(i) = MIN(0.5,zqs(i))
-           zcor = 1./(1.-RETV*zqs(i))
-           zqs(i) = zqs(i)*zcor
-           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
-           zdqsdT_raw(i) = zdqs(i)*  &
-         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
-        ENDDO
-     ELSE
-        DO i = 1, klon
-           IF (zt(i).LT.t_coup) THEN
-              zqs(i) = qsats(zt(i))/pplay(i,k)
-              zdqs(i) = dqsats(zt(i),zqs(i))
-           ELSE
-              zqs(i) = qsatl(zt(i))/pplay(i,k)
-              zdqs(i) = dqsatl(zt(i),zqs(i))
-           ENDIF
-        ENDDO
-     ENDIF
-     !
-     ! Determiner la condensation partielle et calculer la quantite
-     ! de l'eau condensee:
-     !
-!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
-!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
-!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
-!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
-!         stop
-!       endif
-
-     ! ----------------------------------------------------------------
-     ! P2> Formation du nuage
-     ! ----------------------------------------------------------------
-     ! Variables calculees:
-     !   rneb  : fraction nuageuse
-     !   zcond : eau condensee moyenne dans la maille.
-     !   rhcl: humidite relative ciel-clair
-     !   zt : temperature de la maille
-     ! ----------------------------------------------------------------
-     !
-     IF (cpartiel) THEN
-        ! -------------------------
-        ! P2.A> Nuage fractionnaire
-        ! -------------------------
-        !
-        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
-        !   nuageuse a partir des PDF de Sandrine Bony.
-        !   rneb  : fraction nuageuse
-        !   zqn   : eau totale dans le nuage
-        !   zcond : eau condensee moyenne dans la maille.
-        !  on prend en compte le réchauffement qui diminue la partie
-        ! condensee
-        !
-        !   Version avec les raqts
-
-        ! ----------------------------------------------------------------
-        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
-        ! ----------------------------------------------------------------
-
-#ifdef ISO          
-#ifdef ISOVERIF
-        do i=1,klon
-          if (iso_eau.gt.0) then    
-            if (iso_verif_egalite_choix_nostop(zxt(iso_eau,i),zq(i), &
-     &            'il pleut 608',errmax,errmaxrel).eq.1) then
-              write(*,*) 'i=',i
-              stop
-            endif
-          endif !if (iso_eau.gt.0) then
-#ifdef ISOTRAC
-          call iso_verif_traceur(zxt(1,i),'il pleut 1106')
-#endif          
-        enddo !do i=1,klon        
-#endif
-#endif 
-
-        if (iflag_pdf.eq.0) then
-
-           ! version creneau de (Li, 1998)
-           do i=1,klon
-              zdelq = min(ratqs(i,k),0.99) * zq(i)
-              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
-              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
-           enddo
-
-        else !  if (iflag_pdf.eq.0)
-           ! ----------------------------------------------------------------
-           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
-           ! les valeurs de T et Q initiales
-           ! ----------------------------------------------------------------
-           do i=1,klon
-              if(zq(i).lt.1.e-15) then
-                 ncoreczq=ncoreczq+1
-#ifdef ISO
-                zq_ancien(i)=zq(i)
-#endif
-                 zq(i)=1.e-15
-#ifdef ISO
-                ! on conserve le même rapport pour les isotopes
-                if (zq_ancien(i).gt.0.0) then
-                   do ixt=1,ntraciso
-                     zxt(ixt,i)=zxt(ixt,i)/zq_ancien(i)*zq(i)
-                     zxt(ixt,i)=max(0.0,zxt(ixt,i))
-                   enddo   
-                else !if (zq_ancien(i).gt.0.0) then
-                    do ixt=1,ntraciso
-                     zxt(ixt,i)=0.0
-                   enddo 
-                endif !if (zq_ancien(i).gt.0.0) then
-
-                ! verif
-#ifdef ISOVERIF
-                  if (iso_eau.gt.0) then
-                    call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
-     &                 'il pleut 654',errmax,errmaxrel)
-                  endif !if (iso_eau.gt.0) then
-#ifdef ISOTRAC
-             call iso_verif_traceur(zxt(1,i),'il pleut 1207')
-#endif                    
-#endif
-                ! end verif
-#endif 
-              endif
-           enddo
-
-           if (iflag_cld_th>=5) then
-
-              if (iflag_cloudth_vert<=2) then
-               call cloudth(klon,klev,k,ztv, &
-                   zq,zqta,fraca, &
-                   qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
-                   ratqs,zqs,t, &
-                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
-
-              elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
-               call cloudth_v3(klon,klev,k,ztv, &
-                   zq,zqta,fraca, &
-                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
-                   ratqs,zqs,t, &
-                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
-
-              !----------------------------------
-              !Version these Jean Jouhaud, Decembre 2018 
-              !----------------------------------              
-              elseif (iflag_cloudth_vert==6) then
-               call cloudth_v6(klon,klev,k,ztv, &
-                   zq,zqta,fraca, &
-                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
-                   ratqs,zqs,t, &
-                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
-
-
-              endif
-              do i=1,klon
-                 rneb(i,k)=ctot(i,k)
-                 rneblsvol(i,k)=ctot_vol(i,k)
-                 zqn(i)=qcloud(i)
-              enddo
-
-           endif
-
-           if (iflag_cld_th <= 4) then
-              lognormale = .true.
-           elseif (iflag_cld_th >= 6) then
-              ! lognormale en l'absence des thermiques
-              lognormale = fraca(:,k) < 1e-10
-           else
-              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
-              ! bi-gaussienne
-              lognormale = .false.
-           end if
-
-!CR: variation de qsat avec T en presence de glace ou non
-!initialisations
-           do i=1,klon
-              DT(i) = 0.
-              n_i(i)=0
-              Tbef(i)=zt(i)
-              qlbef(i)=0.
-           enddo
-
-        ! ----------------------------------------------------------------
-        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
-        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
-        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
-        ! qui en dependent. Ce changement de temperature est provisoire, et
-        ! la valeur definitive sera calcule plus tard.
-        ! Variables calculees: 
-        !   rneb : nebulosite
-        !   zcond: eau condensee en moyenne dans la maille
-        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
-        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
-        ! T sur qsat
-        ! ----------------------------------------------------------------
-
-!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
-           if (iflag_fisrtilp_qsat.ge.0) then
-             ! Iteration pour condensation avec variation de qsat(T) 
-             ! -----------------------------------------------------
-             do iter=1,iflag_fisrtilp_qsat+1
-                
-               do i=1,klon
-!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
-                 ! !! convergence = .true. tant que l'on n'a pas converge !!
-                 ! ------------------------------
-                 convergence(i)=abs(DT(i)).gt.DDT0
-                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
-                 ! si on n'a pas converge
-                 !
-                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
-                 ! ---------------------------------------------------------------
-                 ! Variables calculees: 
-                 ! rneb : nebulosite
-                 ! zqn : eau condensee, dans le nuage (in cloud water content)
-                 ! zcond: eau condensee en moyenne dans la maille
-                 ! rhcl: humidite relative ciel-clair
-                 !
-                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
-                 if (.not.ice_thermo) then   
-                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
-                 else
-                    if (iflag_t_glace.eq.0) then
-                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
-                    else if (iflag_t_glace.ge.1) then
-                       if (iflag_oldbug_fisrtilp.EQ.0) then
-                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
-                       else
-                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
-                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
-                       endif
-                    endif
-                 endif
-                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
-                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-               if (fl_cor_ebil .GT. 0) then
-                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-               else
-                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
-               end if
-                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
-                 zqs(i) = MIN(0.5,zqs(i))
-                 zcor = 1./(1.-RETV*zqs(i))
-                 zqs(i) = zqs(i)*zcor
-                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
-                 zpdf_sig(i)=ratqs(i,k)*zq(i)
-                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
-                 zpdf_delta(i)=log(zq(i)/zqs(i))
-                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
-                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
-                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
-                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
-                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
-                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
-                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
-                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
-             
-                 if (zpdf_e1(i).lt.1.e-10) then
-                    rneb(i,k)=0.
-                    zqn(i)=zqs(i)
-                 else
-                    rneb(i,k)=0.5*zpdf_e1(i)
-                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
-                 endif
-
-                 ! If vertical heterogeneity, change fraction by volume as well
-                 if (iflag_cloudth_vert>=3) then
-                   ctot_vol(i,k)=rneb(i,k)
-                   rneblsvol(i,k)=ctot_vol(i,k)
-                 endif
-
-                 endif !convergence
-
-               enddo ! boucle en i
-
-                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT 
-                 !         due a la condensation.
-                 ! ---------------------------------------------------------------
-                 ! Variables calculees: 
-                 ! DT : variation de temperature due a la condensation 
-
-                 if (.not. ice_thermo) then
-                 ! --------------------------
-                 do i=1,klon
-                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
-
-                 qlbef(i)=max(0.,zqn(i)-zqs(i))
-               if (fl_cor_ebil .GT. 0) then
-                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
-               else
-                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
-               end if
-                 denom=1.+rneb(i,k)*zdqs(i)
-                 DT(i)=num/denom
-                 n_i(i)=n_i(i)+1
-                 endif
-                 enddo
-
-                 else ! if (.not. ice_thermo)
-                 ! --------------------------
-                 if (iflag_t_glace.ge.1) then
-                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
-                 endif
-
-                 do i=1,klon
-                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
-                 
-                 if (iflag_t_glace.eq.0) then
-                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
-                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-                    zfice(i) = zfice(i)**exposant_glace_old
-                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
-          &                     / (t_glace_min_old - RTT)
-                 endif
-                 
-                 if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
-                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
-          &                    / (t_glace_min - t_glace_max)
-                 endif
-               
-                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
-                    dzfice(i)=0.
-                 endif
-
-                 if (zfice(i).lt.1) then
-                    cste=RLVTT
-                 else
-                    cste=RLSTT
-                 endif
-
-                 qlbef(i)=max(0.,zqn(i)-zqs(i))
-               if (fl_cor_ebil .GT. 0) then
-                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
-           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
-                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
-                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
-           &               *qlbef(i)*dzfice(i)
-               else
-                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
-           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
-                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
-                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
-               end if
-                 DT(i)=num/denom
-                 n_i(i)=n_i(i)+1
-
-                 endif ! fin convergence
-                 enddo ! fin boucle i
-
-                 endif !ice_thermo
-
-             enddo ! iter=1,iflag_fisrtilp_qsat+1
-             ! Fin d'iteration pour condensation avec variation de qsat(T) 
-             ! -----------------------------------------------------------
-           endif !  if (iflag_fisrtilp_qsat.ge.0)
-     ! ----------------------------------------------------------------
-     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
-     ! ----------------------------------------------------------------
-
-        endif ! iflag_pdf
-
-!        if (iflag_fisrtilp_qsat.eq.-1) then 
-!------------------------------------------
-!CR: ATTENTION option fausse mais a existe:
-! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
-! activer les lignes suivantes: 
-       IF (1.eq.0) THEN 
-       DO i=1,klon
-           IF (rneb(i,k) .LE. 0.0) THEN
-              zqn(i) = 0.0
-              rneb(i,k) = 0.0
-              zcond(i) = 0.0
-              rhcl(i,k)=zq(i)/zqs(i)
-           ELSE IF (rneb(i,k) .GE. 1.0) THEN
-              zqn(i) = zq(i)
-              rneb(i,k) = 1.0                 
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
-              rhcl(i,k)=1.0
-           ELSE
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
-              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-           ENDIF
-       ENDDO
-       ENDIF
-!------------------------------------------
-
-!        ELSE
-        ! ----------------------------------------------------------------
-        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
-        ! Variables calculees: 
-        !   rneb : nebulosite
-        !   zcond: eau condensee en moyenne dans la maille
-        !   zq : eau vapeur dans la maille
-        !   zt : temperature de la maille
-        !   rhcl: humidite relative ciel-clair
-        ! ----------------------------------------------------------------
-        !
-        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb) 
-        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
-        ! et de l'humidite relative ciel-clair (rhcl)
-        DO i=1,klon
-           IF (rneb(i,k) .LE. 0.0) THEN
-              zqn(i) = 0.0
-              rneb(i,k) = 0.0
-              zcond(i) = 0.0
-              rhcl(i,k)=zq(i)/zqs(i)
-           ELSE IF (rneb(i,k) .GE. 1.0) THEN
-              zqn(i) = zq(i)
-              rneb(i,k) = 1.0
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
-              rhcl(i,k)=1.0
-           ELSE
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-           ENDIF
-        ENDDO
-        ! If vertical heterogeneity, change fraction by volume as well
-        if (iflag_cloudth_vert>=3) then
-          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
-          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
-        endif
-
-!        ENDIF
-
-     ELSE ! de IF (cpartiel) 
-        ! -------------------------
-        ! P2.B> Nuage "tout ou rien"
-        ! -------------------------
-        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
-        DO i = 1, klon
-           IF (zq(i).GT.zqs(i)) THEN
-              rneb(i,k) = 1.0
-           ELSE
-              rneb(i,k) = 0.0
-           ENDIF
-           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
-        ENDDO
-     ENDIF ! de IF (cpartiel) 
-
-#ifdef ISO
-        ! on calcule la compo iso de l'eau dans le nuage et en ciel
-        ! clair: zxtn/zqn et zxtcs/zqcs
-
-              ! Hypothèsesimplificatrice mais peut-être fausse: 
-              ! eau nuageuse a même composition que eau de la maille 
-
-              ! ce calcul n'eat pas sensé dépendre du type de PDF utilisé     
-         do i=1,klon  
-#ifdef ISOTRAC
-            if (option_tmin.eq.1) then
-                ! c'est mieux si zqn>=0.0
-            if (essai_convergence) then
-            else
-                zqn(i)=max(0.0,zqn(i))
-            endif    
-
-            if (rneb(i,k).lt.1.0) then
-                zqcs(i)=(zq(i)-zqn(i)*rneb(i,k))/(1.0-rneb(i,k))
-#ifdef ISOVERIF
-!                write(*,*) 'rneb,zq,zqn,zqcs=',rneb(i,k),zq(i), &
-!     &                zqn(i),zqcs(i)
-                call iso_verif_positif(zqcs(i),'ipleut 853')
-#endif                
-            else
-                zqcs(i)=0.0
-            endif
-            endif ! if (option_tmin.eq.1) then
-#endif                  
-#ifdef ISOVERIF
-              if (iso_eau.gt.0) then  
-                 call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
-     &                   'il pleut 748',errmax,errmaxrel) 
-              endif !if (iso_eau.gt.0) then
-              do ixt=1,ntraciso
-                call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 857')
-              enddo 
-#ifdef ISOTRAC
-              call iso_verif_traceur_pbidouille(zxt(1,i), &
-     &           'il pleut 874')
-#endif              
-#endif
-              if ((bidouille_anti_divergence).and. &
-     &           (iso_eau.gt.0)) then
-                zxt(iso_eau,i)= max(zq(i),0.0)  
-!                if (zxt(iso_eau,i).gt.0.0) then
-!                    do ixt=1,niso
-!                        zxt(ixt,i)=zxt(ixt,i)*zq(i)/zxt(iso_eau,i)
-!                    enddo
-!                else !if (zxt(iso_eau,i).gt.0.0) then
-!                    do ixt=1,niso
-!                        zxt(ixt,i)=0.0
-!                    enddo    
-!                endif !if (zxt(iso_eau,i).gt.0.0) then
-              endif !if ((bidouille_anti_divergence).and. 
-
-              if (zq(i).gt.1e-15) then   
-                ! 2 fev 2010: on modifie seuil de zq  
-                ! avant, c'était 0. Vérifier que ça n'a pas de
-                ! conséquences trop graves pour les isotopes.
-                do ixt=1,ntraciso             
-                   zxtn(ixt,i)=zqn(i)*(zxt(ixt,i)/zq(i))  
-                enddo !do ixt=1,niso         
-#ifdef ISOTRAC
-               if (option_tmin.eq.1) then
-                 do ixt=1,ntraciso
-                   zxtcs(ixt,i)=(zqcs(i))*(zxt(ixt,i)/zq(i))
-                 enddo !do ixt=1,ntraciso
-               endif ! if (option_tmin.eq.1) then
-#endif                
-              else !if (zq(i).gt.0.0) then 
-                  ! problème: zq~0 et sans compo, mais zqn >0
-                  ! on fait ce qu'on peut dans ce cas pathologique
-                do ixt=1,ntraciso           
-                   zxtn(ixt,i)=0.0              
-                enddo !do ixt=1,niso 
-#ifdef ISOTRAC
-                if (option_tmin.eq.1) then
-                 do ixt=1,ntraciso
-                   zxtcs(ixt,i)=0.0
-                 enddo !do ixt=1,niso    
-                endif ! if (option_tmin.eq.1) then
-#endif                     
-                if ((bidouille_anti_divergence) &
-     &                   .and.(iso_eau.gt.0)) then
-!                   write(*,*) 'ilp tmp 848: warning: bidouille'
-                   zxtn(iso_eau,i)=zqn(i)
-#ifdef ISOTRAC
-                   zxtn(itZonIso(izone_poubelle,iso_eau),i)=zqn(i)  
-                   if (option_tmin.eq.1) then                   
-                     zxtcs(iso_eau,i)=zqcs(i)
-                   endif ! if (option_tmin.eq.1) then
-!                   write(*,*) 'zxtn(:,i)=',zxtn(:,i)
-#endif                   
-                endif !if ((bidouille_anti_divergence)
-              endif  !if (zq(i).gt.0.0) then 
-
-#ifdef ISOVERIF
-             do ixt=1,niso     
-                  call iso_verif_positif(zxtn(ixt,i),'il pleut 1062')  
-             enddo !do ixt=1,niso
-             if (iso_eau.gt.0) then
-                   call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
-     &                   'il pleut 932',errmax,errmaxrel)
-             endif !if ((iso_eau.gt.0)) then
-#ifdef ISOTRAC
-!             write(*,*) 'zxt(:,i)=',zxt(:,i)   
-!             write(*,*) 'zq,zqn=',zq(i) ,zqn(i)  
-             call iso_verif_traceur_pbidouille(zxtn(1,i), &
-     &           'il pleut 1277')
-             if (option_tmin.eq.1) then
-              call iso_verif_positif(zxtcs(ixt,i),'il pleut 895')  
-              do ixt=1,niso     
-               call iso_verif_positif(zxtcs(ixt,i),'il pleut 897')  
-              enddo !do ixt=1,niso
-              if (iso_eau.gt.0) then
-               call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
-     &              zqcs(i),'il pleut 733',errmax,errmaxrel)
-              endif !if ((iso_eau.gt.0)) then
-             endif ! if (option_tmin.eq.1) then
-#endif 
-#endif                   
-!                 write(*,*) 'il pleut 695: calcul de zxtn, i=',i
-!                 write(*,*) 'zxtn(4,i)=',zxtn(4,i)
-!                 write(*,*) 'zqn(i)=',zqn(i)
-              do ixt=1,ntraciso    
-               zxtn(ixt,i)=max(0.0,zxtn(ixt,i))          
-              enddo !do ixt=1,niso
-#ifdef ISOTRAC
-              if (option_tmin.eq.1) then
-              do ixt=1,ntraciso
-               zxtcs(ixt,i)=max(0.0,zxtcs(ixt,i))
-              enddo
-              endif ! if (option_tmin.eq.1) then
-#endif                            
-
-#ifdef ISOVERIF
-               if (iso_eau.gt.0) then
-                  call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
-     &                   'il pleut 746',errmax,errmaxrel) 
-               endif !if (iso_eau.gt.0) then
-               if (iso_HDO.gt.0) then
-                if (zqn(i).gt.ridicule) then
-                   call iso_verif_aberrant(zxtn(iso_HDO,i)/zqn(i), &
-     &                   'il pleut 977') 
-                endif  !if (zqn(i).gt.ridicule) then
-               endif  !if (iso_HDO.gt.0) then  
-#ifdef ISOTRAC               
-             call iso_verif_traceur_pbidouille  &
-     &                  (zxtn(1,i),'il pleut 1300')
-             if (option_tmin.eq.1) then
-               if (iso_eau.gt.0) then   
-                   call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
-     &                   zqcs(i), &
-     &                   'il pleut 971',errmax,errmaxrel)
-               endif !if (iso_eau.gt.0) then
-               call iso_verif_traceur_pbidouille &
-     &                  (zxtcs(1,i),'il pleut 1300')
-             
-             endif ! if (option_tmin.eq.1) then
-#endif                
-#endif   
-         enddo ! do i=1,klon
-         
-!        write(*,*) 'ilp tmp 1311: zoliq,zxtoliq(1,363)=',zoliq(363),zxtoliq(1,363)
-#endif    
-     !
-     ! Mise a jour vapeur d'eau
-     ! -------------------------
-     DO i = 1, klon
-        zq(i) = zq(i) - zcond(i)
-        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
-     ENDDO
-!AJ<
-     ! ------------------------------------
-     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
-     ! -------------------------------------
-     ! Variable calcule:
-     !   zt : temperature de la maille
-     !
-     IF (.NOT. ice_thermo) THEN
-        if (iflag_fisrtilp_qsat.lt.1) then
-           DO i = 1, klon
-             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-           ENDDO
-        else if (iflag_fisrtilp_qsat.gt.0) then
-           DO i= 1, klon
-    if (fl_cor_ebil .GT. 0) then
-             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
-    else
-             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
-    end if
-           ENDDO
-        endif
-     ELSE
-         if (iflag_t_glace.ge.1) then
-            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
-         endif
-         if (iflag_fisrtilp_qsat.lt.1) then
-           DO i = 1, klon
-! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
-!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                     t_glace_max, exposant_glace)
-              if (iflag_t_glace.eq.0) then
-                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
-                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-                    zfice(i) = zfice(i)**exposant_glace_old
-              endif
-              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
-                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
-           ENDDO
-         else
-           DO i=1, klon
-! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
-!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                     t_glace_max, exposant_glace)
-              if (iflag_t_glace.eq.0) then
-                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
-                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-                    zfice(i) = zfice(i)**exposant_glace_old
-              endif
-        if (fl_cor_ebil .GT. 0) then
-              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 
-           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
-                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
-        else 
-              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
-                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
-        end if
-           ENDDO
-         endif
-!         print*,zt(i),zrfl(i),zifl(i),'temp1'
-       ENDIF
-!>AJ
-
-      ! isotopes lors de la condensation
-#ifdef ISO          
-#ifdef ISOVERIF
-!        write(*,*) 'ilp tmp 1823'
-      do i=1,klon
-        ! verif
-        if (iso_eau.gt.0) then
-          call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
-     &           'il pleut 810',errmax,errmaxrel) 
-        endif !if (iso_eau.gt.0) then 
-#ifdef ISOTRAC
-        call iso_verif_traceur(zxtn(1,i),'il pleut 1408')
-!        write(*,*) 'fisrtilp 1421: appel condiso_liq_ice_vectall'
-#endif
-       enddo !do i=1,klon 
-!       write(*,*) 'il pleut 961: avant condensation'
-       ! end verif
-#endif
-      
-      do i=1,klon  
-!        call calcul_zfice(zt(i),zfice(i)) ! inliner pour optimisation:
-        zfice(i)=1.0-(zt(i)-pxtice)/(pxtmelt-pxtice)
-        zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-       enddo  
-        
-        ! calcul de la composition du condensat
-        ! hypothèse simplificatrice (mais à garder en mémoire!): l'eau
-        ! nuageuse a la même composition que l'eau totale de la maille:
-        ! c'est à dire zxtn=(zxt/zq)*zqn
-        
-        call condiso_liq_ice_vectall(zxtn,zqn,zcond,zt,zfice, &
-     &            zxtice,zxtliq,klon)
-#ifdef ISOTRAC
-!#ifdef ISOVERIF
-!        write(*,*) 'il pleut 980: condensation pour les traceurs'
-!#endif        
-        call condiso_liq_ice_vectall_trac(zxtn,zqn,zcond, &
-     &            zt,zfice,zxtice,zxtliq,klon)
-#ifdef ISOVERIF
-!        write(*,*) 'ilp tmp 1859'
-        do i=1,klon             
-             call iso_verif_traceur(zxtn(1,i),'il pleut 1428a')
-             call iso_verif_traceur(zxtliq(1,i),'il pleut 1428b')
-             call iso_verif_traceur(zxtice(1,i),'il pleut 1428c')
-        enddo
-#endif
-#endif
-        do i=1,klon
-        do ixt=1,ntraciso                 
-          zxtcond(ixt,i)=zxtice(ixt,i)+zxtliq(ixt,i)
-          zxtcond(ixt,i)=max(zxtcond(ixt,i),0.0)
-          zxtcond(ixt,i)=min(zxtcond(ixt,i),zxt(ixt,i))
-        enddo !do ixt=1,niso 
-        enddo !do i=1,klon
-
-       do i=1,klon   
-       enddo  !do i=1,klon       
-        
-!#ifdef ISOVERIF
-#ifdef ISOVERIF
-        do i=1,klon
-          do ixt=1,ntraciso  
-            call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 638')
-          enddo
-        enddo !do i=1,klon 
-#endif       
-#ifdef ISOVERIF       
-        if (iso_eau.gt.0) then
-          do i=1,klon      
-              call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
-     &                   'il pleut 645',errmax,errmaxrel)
-          enddo ! do i=1,klon     
-        endif
-        if (iso_HDO.gt.0) then
-           do i=1,klon  
-                  call iso_verif_aberrant_choix(zxtcond(iso_hdo,i), &
-     &           zcond(i),ridicule_rain,deltalim_snow,'il pleut 1276')
-           enddo ! do i=1,klon  
-        endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
-#ifdef ISOTRAC
-        do i=1,klon
-             if (iso_verif_traceur_nostop(zxtcond(1,i),'il pleut 1455') &
-     &                   .eq.1) then
-               write(*,*) 'zxtn(:,i)=',zxtn(:,i)
-               write(*,*) 'zxt(:,i)=',zxt(:,i)
-               write(*,*) 'zxtliq(:,i)=',zxtliq(:,i)
-               write(*,*) 'zxtice(:,i)=',zxtice(:,i)
-               stop
-             endif
-        enddo     
-#endif         
-#endif
-
-       ! retranchage à zxt 
-      do i=1,klon 
-       do ixt=1,ntraciso
-          zxt(ixt,i)=zxt(ixt,i)-zxtcond(ixt,i)
-       enddo !do ixt=1,niso   
-      enddo !do i=1,klon 
-
-#ifdef ISOVERIF
-      do i=1,klon
-         do ixt=1,niso
-          call iso_verif_positif(zxt(ixt,i),'il pleut 1088')
-         enddo !do ixt=1,niso
-#ifdef ISOTRAC  
-         ! on sépare pour mieux débugger 8 fev 2010
-         do ixt=1,ntraciso
-          call iso_verif_positif(zxt(ixt,i),'il pleut 1093')
-         enddo !do ixt=1,niso
-#endif      
-      enddo   !do i=1,klon
-#endif
-
-#ifdef ISOTRAC        
-        if (option_tmin.ge.1) then
-        do i=1,klon    
-        ! colorier la vapeur résiduelle selon température de
-        ! condensation, et le condensat en un tag spécifique
-        ! Attention: zqn,zxtn ne servent qu'au calcul de zcond,zxtcond.
-        ! C'est en fait zxt qui se fait retranché
-        ! correction le 31 mars 2010: c'est à qn qu'il faut retrancher
-        ! le condensat, car la condensation LS est un processus sous
-        ! maille -> ne tagguer que le résidus de qn.
-          if ((zcond(i).gt.0.0).and.(zq(i).gt.0.0)) then 
-
-            if (rneb(i,k).gt.0.0) then
-              zcondn(i)=zcond(i)/rneb(i,k)
-              do ixt=1,ntraciso
-                zxtcondn(ixt,i)=zxtcond(ixt,i)/rneb(i,k)
-              enddo !do ixt=1,ntraciso
-            else !if (rneb.eq.0.0)
-              zcondn(i)=0.0
-              do ixt=1,ntraciso
-                zxtcondn(ixt,i)=0.0
-              enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-              call iso_verif_egalite(zcond(i),0.0,'ilp 1225')
-#endif              
-            endif !if (rneb.eq.0.0)
-            
-#ifdef ISOVERIF
-          call iso_verif_egalite_choix(zxtcondn(iso_eau,i),zcondn(i), &
-     &                   'il pleut 1244',errmax*10,errmaxrel*10)
-          call iso_verif_traceur(zxtcondn(1,i),'cv3_routines 2302')
-          zq_verif=rneb(i,k)*zqn(i)+(1.0-rneb(i,k))*zqcs(i)-zcond(i)
-          call iso_verif_egalite_choix(zq(i),zq_verif, &
-     &                   'il pleut 1235',errmax,errmaxrel)
-          do ixt=1,niso
-            zq_verif=rneb(i,k)*zxtn(ixt,i) &
-     &           +(1.0-rneb(i,k))*zxtcs(ixt,i)-zxtcond(ixt,i)
-            call iso_verif_egalite_choix(zxt(ixt,i),zq_verif, &
-     &                   'il pleut 1235',errmax,errmaxrel)
-          enddo
-#endif            
-          ! bidouille
-          if (bidouille_anti_divergence) then
-                zxtcondn(iso_eau,i)=zcondn(i)
-#ifdef ISOVERIF                
-                call iso_verif_traceur_pbidouille(zxtcondn(1,i), &
-     &          'il pleut 1261')
-#endif                
-          endif
-
-            if (option_traceurs.eq.17) then      
-            call iso_recolorise_condensation(zqn(i),zcondn(i), &
-     &           zxtn(1,i),zxtcondn(1,i),zt(i),1.0,zxtres, &
-     &           seuil_tag_tmin_ls)
-            else !if (option_traceurs.eq.17) then
-                ! cas du tag 18 ou 19
-                call iso_recolorise_condensation(zqn(i),zcondn(i), &
-     &           zxtn(1,i),zxtcondn(1,i),zqs(i),1.0,zxtres, &
-     &           seuil_tag_tmin_ls)
-            endif !if (option_traceurs.eq.17) then
-            do ixt=1+niso,ntraciso
-               zxt(ixt,i)=rneb(i,k)*zxtres(ixt) &
-     &           +(1.0-rneb(i,k))*zxtcs(ixt,i)
-            enddo         
-#ifdef ISOVERIF
-         call iso_verif_traceur(zxt(1,i),'ilp 1025')
-         do izone=izone_oce,izone_ddft
-          if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
-           ! avec condensation, on ne peut que perdre des traceurs à
-           ! part le tag résuel et le condensat
-           if (iso_verif_positif_choix_nostop( &
-     &           zxt_ancien(itZonIso(izone,iso_eau),i) &
-     &          -zxt(itZonIso(izone,iso_eau),i),1e-8,'ilp 1270') &
-     &          .eq.1) then
-            write(*,*) 'i,izone,rneb=',i,izone,rneb(i,k)
-            write(*,*) 'zxt(:,i)=',zxt(:,i)
-            write(*,*) 'zxt_ancien(:,i)=',zxt_ancien(:,i)
-            write(*,*) 'zxtcs(:,i)=',zxtcs(:,i)
-            write(*,*) 'zxtres(:)=',zxtres(:)
-            write(*,*) 'zxtcond(:,i)=',zxtcond(:,i)
-            write(*,*) 'zxtn(:,i)=',zxtn(:,i)
-            write(*,*) 'zxtcondn(:,i)=',zxtcondn(:,i)
-            stop
-           endif !if (iso_verif_positif_nostop(
-          endif !if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
-         enddo !do izone=izone_oce,izone_revap
-#endif                     
-          endif !if (cond.gt.0.0) then
-        enddo !do i=1,klon
-
-#ifdef ISOVERIF
-        ! vérifier qu'on recolorise bien le condensat
-        call iso_verif_positif(float(option_cond-1),'ilp 1310')
-#endif        
-        endif ! if (option_tmin.ge.1) then
-
-        if (option_cond.eq.1) then
-        ! colorier le condensat en un tag spécifique
-        do i=1,klon        
-          do ixt=1+niso,ntraciso
-            if (index_zone(ixt).eq.izone_cond) then
-              zxtcond(ixt,i)=zxtcond(index_iso(ixt),i)
-            else
-              zxtcond(ixt,i)=0.0
-            endif
-          enddo !do ixt=1,ntraciso
-        enddo !do i=1,klon
-#ifdef ISOVERIF
-        do i=1,klon    
-          call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
-     &                   'il pleut 1139',errmax,errmaxrel)
-          call iso_verif_traceur(zxtcond(1,i),'cv3_routines 2302')
-        enddo !do i=1,klon        
-#endif         
-      endif ! if (option_cond.ge.1) then
-#endif        
-     
-      
-#ifdef ISOVERIF
-!        write(*,*) 'ilp tmp 2066'
-      do i=1,klon
-         do ixt=1,niso
-          call iso_verif_positif(zxt(ixt,i),'il pleut 778')
-         enddo !do ixt=1,niso
-#ifdef ISOTRAC  
-         ! on sépare pour mieux débugger 8 fev 2010
-         do ixt=1,ntraciso
-          call iso_verif_positif(zxt(ixt,i),'il pleut 1135')
-          call iso_verif_traceur(zxt(1,i),'il pleut 1485')
-         enddo !do ixt=1,niso
-#endif
-         if (iso_eau.gt.0) then
-            call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
-     &                   'il pleut 644',errmax,errmaxrel)
-         endif !if (iso_eau.gt.0) then
-         if (iso_HDO.gt.0) then  
-            if (zq(i).gt.ridicule) then
-             call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
-     &                   'il pleut 648')                 
-              if (iso_O18.gt.0) then                
-                if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
-     &           zxt(iso_O18,i)/zq(i),'il pleut 2059').eq.1) then
-                   write(*,*) 'i,k,zq=',i,k,zq(i)
-                   stop
-                endif !if (iso_verif_O18_aberrant_nostop
-              endif !if (iso_O18.gt.0) then
-            endif !if (zq(i).gt.ridicule) then
-         endif  ! if (iso_HDO.gt.0) then 
-        enddo  ! i=1,klon
-#ifdef ISOTRAC
-        do i=1,klon
-             call iso_verif_traceur(zxt(1,i),'il pleut 1485')
-        enddo     
-#endif         
-#endif
-#ifdef ISOVERIF   
-       do i=1,klon   
-           do ixt=1,ntraciso
-             call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 1560')
-             call iso_verif_noNAN(zxt(ixt,i),'il pleut 1561')
-           enddo   ! ixt=1,niso
-        enddo  !do i=1,klon         
-#endif   
-#endif    
-     ! ----------------------------------------------------------------
-     ! P3> Formation des precipitations
-     ! ----------------------------------------------------------------
-     !
-     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
-     !
-
-     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
-     DO i = 1, klon
-        IF (rneb(i,k).GT.0.0) THEN
-           zoliq(i) = zcond(i)
-           zrho(i) = pplay(i,k) / zt(i) / RD
-           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
-        ENDIF
-     ENDDO
-!AJ<
-     IF (.NOT. ice_thermo) THEN
-       IF (iflag_t_glace.EQ.0) THEN
-         DO i = 1, klon
-            IF (rneb(i,k).GT.0.0) THEN
-               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
-               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-               zfice(i) = zfice(i)**exposant_glace_old
-!              zfice(i) = zfice(i)**nexpo
-         !!      zfice(i)=0.
-            ENDIF
-         ENDDO
-       ELSE ! of IF (iflag_t_glace.EQ.0)
-         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
-!         DO i = 1, klon
-!            IF (rneb(i,k).GT.0.0) THEN
-! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
-!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                     t_glace_max, exposant_glace)
-!            ENDIF
-!         ENDDO
-       ENDIF
-     ENDIF
-
-     ! Calcul de radliq (eau condensee pour le rayonnement)
-     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
-     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
-     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
-     ! transmise au rayonnement;
-     ! ----------------------------------------------------------------
-     DO i = 1, klon
-        IF (rneb(i,k).GT.0.0) THEN
-           zneb(i) = MAX(rneb(i,k), seuil_neb)
-     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i))  
-!           print*,zt(i),'fractionglace'
-!>AJ
-           radliq(i,k) = zoliq(i)/REAL(ninter+1)
-        ENDIF
-     ENDDO
-#ifdef ISO
-      ! initilsation des quantités totales d'eau liquide et solide qui
-      ! tombent
-      do i=1,klon
-         if (rneb(i,k).gt.0.0) then
-            do ixt=1,ntraciso
-               zxtoliq(ixt,i)=zxtcond(ixt,i)
-            enddo
-         endif ! if (rneb(i,k).gt.0.0)
-       enddo
-          
-#ifdef ISOVERIF
-!       write(*,*) 'ilp tmp 1626: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
-       do i=1,klon   
-         do ixt=1,ntraciso
-           if (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
-     &           'il pleut 1628').eq.1) then
-             write(*,*) 'i,k,ixt=',i,k,ixt
-             write(*,*) 'rneb(i,k)=',rneb(i,k)
-             write(*,*) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
-             stop
-           endif
-         enddo   ! ixt=1,niso
-        enddo  !do i=1,klon  
-#endif                
-#ifdef ISOVERIF
-       !i=519
-       !write(*,*) 'k,i,zoliq(i)=',k,i,zoliq(i)
-       do i=1,klon   
-           ! cam verif
-           if (iso_eau.gt.0) then
-            call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
-     &            'il pleut 678',errmax,errmaxrel)
-           endif !if (iso_eau.gt.0) then
-           if (zoliq(i).gt.ridicule) then
-            if (iso_HDO.gt.0) then  
-                ! Camille 9 mars 2023: on est moins stricte avec le condensat
-                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
-     &           ridicule_rain,deltalim_snow, 'il pleut 895a')
-              if (iso_O18.gt.0) then
-                if (iso_verif_O18_aberrant_nostop(zxtoliq(iso_HDO,i)/zoliq(i), &
-     &           zxtoliq(iso_O18,i)/zoliq(i),'il pleut 895b').eq.1) then
-                   write(*,*) 'i,k,zoliq,zfice=',i,k,zoliq(i),zfice(i)
-                   !stop
-                endif
-              endif ! if (iso_HDO.gt.0) then
-            endif  !if (iso_HDO.gt.0) then
-           endif !if (zoliq(i).gt.ridicule) then
-          ! end cam verif
-         enddo ! do i=1,klon
-#endif
-#endif
-
-     !
-     DO n = 1, ninter
-        DO i = 1, klon
-           IF (rneb(i,k).GT.0.0) THEN
-              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-              ! Initialization of zpluie and zice:
-              zpluie=0
-              zice=0
-              IF (zneb(i).EQ.seuil_neb) THEN
-                 ztot = 0.0
-              ELSE
-                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
-                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
-                 if (ptconv(i,k)) then
-                    zcl   =cld_lc_con
-                    zct   =1./cld_tau_con
-                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
-                         *fallvc(zrhol(i)) * zfice(i)
-                 else
-                    zcl   =cld_lc_lsc
-                    zct   =1./cld_tau_lsc
-                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
-                         *fallvs(zrhol(i)) * zfice(i)
-                 endif
-
-                 ! si l'heterogeneite verticale est active, on utilise 
-                 ! la fraction volumique "vraie" plutot que la fraction
-                 ! surfacique modifiee, qui est plus grande et reduit
-                 ! sinon l'eau in-cloud de facon artificielle
-                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
-                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
-                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
-                 else
-                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
-                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
-                 endif
-!AJ<
-                 IF (.NOT. ice_thermo) THEN
-                   ztot    = zchau + zfroi
-                 ELSE
-                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
-                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
-                   ztot    = zpluie    + zice
-                 ENDIF
-!>AJ
-                 ztot    = MAX(ztot   ,0.0)
-              ENDIF
-              ztot    = MIN(ztot,zoliq(i))
-!AJ<
-     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
-     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
-!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
-!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
-!      si iflag_bergeron <> 2
-!      A SUPPRIMER A TERME
-              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
-              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
-              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
-!>AJ
-              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
-           ENDIF
-        ENDDO  ! i = 1,klon
-     ENDDO     ! n = 1,ninter
-
-     ! ----------------------------------------------------------------
-     !
-#ifdef ISO
-      do i=1,klon
-       IF (rneb(i,k).GT.0.0) THEN
-        if (zcond(i).gt.ridicule**2) then
-            ! le 21 dec: on change .gt.0 en .gt.ridicule**2 pour éviter valeur
-            ! ridiculement petites
-           do ixt=1,ntraciso
-             zxtoliq(ixt,i)=zoliq(i)*(zxtcond(ixt,i)/zcond(i))
-           enddo !do ixt=1,niso           
-        else !if (zcond(i).gt.0.0) then            
-           if (zoliq(i).lt.ridicule) then
-               do ixt=1,ntraciso
-                zxtoliq(ixt,i)=0.0
-               enddo
-           else !if (zoliq(i).lt.ridicule) then
-               ! modif 28 oct 2008
-              do ixt=1,niso
-                zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
-              enddo !do ixt=1,niso
-#ifdef ISOTRAC
-              do ixt=niso+1,ntraciso
-                if (index_zone(ixt).eq.izone_poubelle) then
-                    zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
-                else
-                    zxtoliq(ixt,i)=0.0
-                endif
-              enddo  
-#endif    
-#ifdef ISOVERIF                   
-              write(*,*) 'il pleut 1412: zcond(i)=',zcond(i)
-              write(*,*) 'zoliq(i)=',zoliq(i)
-              stop
-#endif                   
-            endif   !if (zoliq(i).lt.ridicule) then         
-        endif !if (zcond(i).gt.0.0) then
-       endif !IF (rneb(i,k).GT.0.0) THEN
-      enddo ! do i=1,klon 
-           
-         ! cam verif
-#ifdef ISOVERIF
-      do i=1,klon  
-        IF (rneb(i,k).GT.0.0) THEN 
-           do ixt=1,niso
-             call iso_verif_noNAN(zxtoliq(ixt,i),'il pleut 718')
-           ! end cam verif
-           enddo ! do ixt=1,niso
-          ! cam verif
-        endif !IF (rneb(i,k).GT.0.0) THEN 
-      enddo !do i=1,klon  
-#endif
-#ifdef ISOVERIF
-      do i=1,klon  
-        IF (rneb(i,k).GT.0.0) THEN
-          if (iso_eau.gt.0) then
-            call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
-     &           'il pleut 749',errmax,errmaxrel)
-          endif !if (iso_eau.gt.0) then
-          if (iso_HDO.gt.0) then
-                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
-     &           ridicule_rain, deltalim_snow, 'il pleut 963')
-          endif !if (iso_HDO.gt.0) then
-        ! end cam verif
-#ifdef ISOTRAC
-!             write(*,*) 'ilp 1201 tmp: i=',i
-!             write(*,*) 'zoliq(i)=',zoliq(i)
-!             write(*,*) 'zxtoliq(i)=',zxtoliq(iso_eau:ntraciso:3,i)
-!             write(*,*) 'zcond(i)=',zcond(i)
-!             write(*,*) 'zxtcond(i)=',zxtcond(iso_eau:ntraciso:3,i)
-             call iso_verif_traceur(zxtoliq(1,i),'il pleut 1634')
-#endif        
-         endif !IF (rneb(i,k).GT.0.0) THEN
-       enddo ! do i=1,klon
-#endif        
-#endif      
-     IF (.NOT. ice_thermo) THEN
-       DO i = 1, klon
-         IF (rneb(i,k).GT.0.0) THEN
-           d_ql(i,k) = zoliq(i)
-           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
-                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-
-#ifdef ISO   
-         do ixt=1,ntraciso
-            d_xtl(ixt,i,k)=zxtoliq(ixt,i)
-         enddo ! do ixt=1,niso   
-            ! Rq: 28 oct 2008: si zoliq a la même compo que zcond, alors
-            ! on évite pb numériques si on prend directement zrfl avec
-            ! la compo de zcond?
-!          do ixt=1,niso  
-!            zxtrfl(ixt,i)=zxtrfl(ixt,i) &
-!     &           +max(zxtcond(ixt,i)-zxtoliq(ixt,i),0.0) &
-!     &           * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-!          enddo  
-            if (zcond(i).gt.ridicule**2) then
-                ! modif le 19 dec 2012 pour éviter zcond pathologiquement petits
-                do ixt=1,ntraciso
-                  zxtrfl(ixt,i)=zxtrfl(ixt,i) &
-     &                  +(zxtcond(ixt,i)/zcond(i)) &
-     &                  * MAX(zcond(i)-zoliq(i),0.0) &
-     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-                enddo !do ixt=1,niso
-            else !if (zcond(i).gt.0.0) then
-                if (zcond(i)-zoliq(i).gt.0.0) then
-                    do ixt=1,niso
-                        zxtrfl(ixt,i)=zxtrfl(ixt,i) &
-     &                  +Rdefault(ixt) &
-     &                  * MAX(zcond(i)-zoliq(i),0.0) &
-     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-                    enddo !do ixt=1,niso  
-#ifdef ISOVERIF                   
-                    write(*,*) 'il pleut 1023: zcond(i)=',zcond(i)
-                    write(*,*) 'zcond(i),zoliq(i)=',zcond(i),zoliq(i)
-                    stop
-#endif        
-#ifdef ISOTRAC
-                   do ixt=niso+1,ntraciso
-                      if (index_zone(ixt).eq.izone_poubelle) then
-                        zxtrfl(ixt,i)=zxtrfl(ixt,i) &
-     &                  +Rdefault(ixt) &
-     &                  * MAX(zcond(i)-zoliq(i),0.0) &
-     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-                      endif
-                    enddo !do ixt=1,niso
-#endif                 
-                else !if (zcond(i)-zoliq(i).gt.0.0) then
-                    ! zrfl non modifié
-                endif !if (zcond(i)-zoliq(i).gt.0.0) then
-            endif !if (zcond(i).gt.0.0) then
-
-         ! cam verif         
-!#ifdef ISOVERIF
-#ifdef ISOVERIF
-            do ixt=1,ntraciso
-             if ((iso_verif_noNAN_nostop(zxtrfl(ixt,i), &
-     &                   'il pleut 772').eq.1).or. &
-     &                   (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
-     &                   'il pleut 772b').eq.1)) then
-               write(*,*) 'zcond(i)=',zcond(i)
-               write(*,*) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
-               write(*,*) 'zoliq(i)=',zoliq(i)
-               write(*,*) 'Rdefault(ixt)=',Rdefault(ixt)
-               write(*,*) 'paprs(i,k),paprs(i,k+1)=',paprs(i,k),paprs(i,k+1)
-               stop
-             endif
-            enddo
-#endif
-#ifdef ISOVERIF
-            if (iso_eau.gt.0) then
-              call iso_verif_egalite_choix(d_xtl(iso_eau,i,k), &
-     &           d_ql(i,k),'il pleut 771',errmax,errmaxrel)
-              if (iso_verif_egalite_choix_nostop( &
-     &           zxtrfl(iso_eau,i),zrfl(i), &
-     &           'il pleut 771b',errmax,errmaxrel).eq.1)then
-                 write(*,*) 'trace_cas(i)=',trace_cas(i)
-                 write(*,*) 'zxtcond(iso_eau,i)',zxtcond(iso_eau,i)
-                 write(*,*) 'zcond(i)',zcond(i)
-                 write(*,*) 'zxtoliq(iso_eau,i)',zxtoliq(iso_eau,i)
-                 write(*,*) 'zoliq(i)',zoliq(i)
-                 stop
-              endif
-            endif !if (iso_eau.gt.0) then
-            if (iso_HDO.gt.0) then
-                IF ((t(i,1)) .LT. RTT) THEN
-                call iso_verif_aberrant_choix(zxtrfl(iso_hdo,i), &
-     &           zrfl(i),ridicule_rain,deltalim_snow,'il pleut 1458')
-                endif
-            endif !if (iso_HDO.gt.0) then
-#ifdef ISOTRAC
-             call iso_verif_traceur(zxtrfl(1,i),'il pleut 1729')
-#endif            
-#endif
-         ! bidouille
-         if ((bidouille_anti_divergence).and.(iso_eau.gt.0)) then 
-               zxtrfl(iso_eau,i)= zrfl(i)   
-!             if (zxtrfl(iso_eau,i).gt.0.0) then
-!             do ixt=1,niso
-!               zxtrfl(ixt,i)=zxtrfl(ixt,i)*zrfl(i)/zxtrfl(iso_eau,i)
-!             enddo
-!           endif !if (zxtrfl(iso_eau,i).gt.0.0) then
-         endif ! if (bidouille_anti_divergence) then
-         ! end bidouille
-         ! end cam verif
-#endif
-         ENDIF !IF (rneb(i,k).GT.0.0) THEN
-       ENDDO
-     ELSE
-
-#ifdef ISO
-        CALL abort_physic('ilp 2347', 'isotopes pas encore prevus ici', 1)
-#endif
-!
-!CR&JYG<
-! On prend en compte l'effet Bergeron dans les flux de precipitation :
-! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
-! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
-! et les precipitations est grossierement pris en compte en linearisant les equations 
-! et en approximant le processus de precipitation liquide par un processus a seuil.
-! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
-! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
-! Le condensat precipitant solide est augmente.
-! La vapeur d'eau est augmentee.
-!
-       IF ((iflag_bergeron .EQ. 2)) THEN
-         DO i = 1, klon
-           IF (rneb(i,k) .GT. 0.0) THEN
-             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
-             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
-           if (fl_cor_ebil .GT. 0) then
-             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
-             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
-!            Calcul de DT si toute les precips liquides congelent
-             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
-!            On ne veut pas que T devienne superieur a la temp. de congelation.
-!            donc que Delta > RTT-zt(i
-             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
-             zt(i)      = zt(i)      + DeltaT
-!            Eau vaporisee du fait de l'augmentation de T
-             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
-!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
-             zq(i) = zq(i) +  Deltaq
-!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
-             zcond(i)   = max( zcond(i)- Deltaq, 0. )
-!            precip liquide qui congele ou qui s'evapore
-             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
-             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
-!            bilan eau glacee
-             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
-           else ! if (fl_cor_ebil .GT. 0)
-!            ancien calcul
-             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
-             coef1 = RLMLT*zdqs(i)/RLVTT
-             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
-             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
-             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
-             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
-             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
-             zt(i)      = zt(i)      + DeltaT
-           end if ! if (fl_cor_ebil .GT. 0)
-           ENDIF  ! rneb(i,k) .GT. 0.0
-         ENDDO
-         DO i = 1, klon
-           IF (rneb(i,k).GT.0.0) THEN
-             d_ql(i,k) = (1-zfice(i))*zoliq(i)
-             d_qi(i,k) = zfice(i)*zoliq(i)
-             zrfl(i) = zrfl(i)+ zqprecl(i) &
-                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-             zifl(i) = zifl(i)+ zqpreci(i) &
-                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
-           ENDIF                      
-         ENDDO
-!!
-       ELSE  ! iflag_bergeron
-!>CR&JYG
-!!
-       DO i = 1, klon
-         IF (rneb(i,k).GT.0.0) THEN
-!CR on prend en compte la phase glace
-!JLD inutile car on ne passe jamais ici si .not.ice_thermo
-!           if (.not.ice_thermo) then
-!           d_ql(i,k) = zoliq(i)
-!           d_qi(i,k) = 0.
-!           else
-           d_ql(i,k) = (1-zfice(i))*zoliq(i)
-           d_qi(i,k) = zfice(i)*zoliq(i)
-!           endif
-!AJ<
-           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
-               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
-                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
-     !      zrfl(i) = zrfl(i)+  zpluie                         &
-     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-     !      zifl(i) = zifl(i)+  zice                    &
-     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                    
-
-!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
-           IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
-              zsolid = zrfl(i)
-              zifl(i) = zifl(i)+zrfl(i)
-              zrfl(i) = 0.
-           if (fl_cor_ebil .GT. 0) then
-              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-           else
-              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
-           end if
-           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
-!RC   
-
-         ENDIF ! rneb(i,k).GT.0.0
-       ENDDO
-
-       ENDIF  ! iflag_bergeron .EQ. 2
-     ENDIF  ! .NOT. ice_thermo
-
-!CR: la fonte est faite au debut
-!      IF (ice_thermo) THEN
-!       DO i = 1, klon
-!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
-!           zmelt = MIN(MAX(zmelt,0.),1.)
-!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
-!           zifl(i)=zifl(i)*(1.-zmelt)
-!           print*,zt(i),'octavio1'
-!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
-!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
-!       ENDDO
-!     ENDIF
-
-       
-     IF (.NOT. ice_thermo) THEN
-       DO i = 1, klon
-         IF (zt(i).LT.RTT) THEN
-           psfl(i,k)=zrfl(i)
-#ifdef ISO
-         do ixt=1,ntraciso
-           pxtsnowfl(ixt,i,k)=zxtrfl(ixt,i)
-         enddo
-#endif
-         ELSE
-           prfl(i,k)=zrfl(i)
-#ifdef ISO
-         do ixt=1,ntraciso
-           pxtrainfl(ixt,i,k)=zxtrfl(ixt,i)
-         enddo         
-#endif
-         ENDIF
-       ENDDO
-     ELSE ! if ice_thermo
-     ! JAM*************************************************
-     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
-     ! *****************************************************
-
-#ifdef ISO
-     CALL abort_physic('ilp 2465', 'isos pas prevus ici', 1)   
-#endif
-       DO i = 1, klon
-     !   IF (zt(i).LT.RTT) THEN
-           psfl(i,k)=zifl(i) 
-     !   ELSE
-           prfl(i,k)=zrfl(i)
-     !   ENDIF
-!>AJ
-       ENDDO
-     ENDIF
-
-#ifdef ISO
-#ifdef ISOVERIF
-      DO i = 1, klon
-        if (iso_eau.gt.0) then
-          call iso_verif_egalite_choix(pxtrainfl(iso_eau,i,k),prfl(i,k), & 
-     &           'il pleut 804',errmax,errmaxrel)
-          call iso_verif_egalite_choix(pxtsnowfl(iso_eau,i,k),psfl(i,k), &  
-     &           'il pleut 805',errmax,errmaxrel) 
-          call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &  
-     &           'il pleut 1612',errmax,errmaxrel) 
-          if (iso_verif_positif_choix_nostop(pxtrainfl(iso_eau,i,k), &
-     &           ridicule_rain,'ilp 1547').eq.1) then
-            write(*,*) 'prfl(i,k)=',prfl(i,k)
-            stop
-          endif  
-        endif !if (iso_eau.gt.0) then
-        if (iso_HDO.gt.0) then
-          if (prfl(i,k).gt.ridicule_rain) then
-            call iso_verif_aberrant(pxtrainfl(iso_HDO,i,k)/prfl(i,k), &
-     &            'il pleut 1020')
-          endif !if (prfl(i,k).gt.ridicule_rain) then           
-          call iso_verif_aberrant_choix(pxtsnowfl(iso_HDO,i,k),psfl(i,k), &
-     &           ridicule_rain, deltalim_snow, 'il pleut 1024')
-          if (zq(i).gt.ridicule) then
-              call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
-     &           'il pleut 1204')
-          endif !if (zq(i).gt.ridicule) then          
-         endif !if (iso_HDO.gt.0) then
-#ifdef ISOTRAC
-          call iso_verif_traceur(zxt(1,i),'il pleut 1792')
-          call iso_verif_traceur(pxtsnowfl(1,i,k),'il pleut 1793')
-          call iso_verif_traceur(pxtrainfl(1,i,k),'il pleut 1794')
-#endif         
-       ENDDO !DO i = 1, klon  
-#endif
-#endif 
-     ! ----------------------------------------------------------------
-     ! Fin de formation des precipitations
-     ! ----------------------------------------------------------------
-     !
-     ! Calculer les tendances de q et de t:
-     !
-     DO i = 1, klon
-        d_q(i,k) = zq(i) - q(i,k)
-#ifdef ISO
-         do ixt=1,ntraciso
-             d_xt(ixt,i,k)=zxt(ixt,i)-xt(ixt,i,k)   
-         enddo
-         ! cam verif
-#ifdef ISOVERIF
-           if (iso_eau.gt.0) then
-             if (iso_verif_egalite_choix_nostop(d_xt(iso_eau,i,k), &
-     &            d_q(i,k),'il pleut 838',errmax,errmaxrel).eq.1) then
-                write(*,*) 'i,k=',i,k
-                write(*,*) 'zxt,zq=',zxt(iso_eau,i),zq(i)
-                write(*,*) 'xt,q=',xt(iso_eau,i,k),q(i,k)
-                stop
-             endif
-           endif !if (iso_eau.gt.0) then
-           if (1.eq.0) then
-           if (iso_HDO.gt.0) then
-            if (abs(d_q(i,k)).gt.max(errmaxrel*q(i,k),errmax)) then
-              if (iso_verif_aberrant_nostop( &
-     &           d_xt(iso_HDO,i,k)/d_q(i,k),'il pleut 1229').eq.1) then
-!                write(*,*) 'd_q(i,k)=',d_q(i,k)
-!                write(*,*) 'q(i,k)=',q(i,k)
-!                write(*,*) 'zq(i)=',zq(i)
-!                write(*,*) 'deltaD_prec=',deltaD(xt(iso_HDO,i,k)/q(i,k))
-!                write(*,*) 'deltaD=',deltaD(zxt(iso_HDO,i)/zq(i))
-!                stop       
-              endif !   if (iso_verif_aberrant_nostop(
-            endif !if (zq(i).gt.ridicule) then          
-           endif !if (iso_HDO.gt.0) then
-          endif !if (1.eq.0) then
-#endif
-         ! end cam verif
-#endif
-        d_t(i,k) = zt(i) - t(i,k)
-     ENDDO
-     !
-     !AA--------------- Calcul du lessivage stratiforme  -------------
-
-     DO i = 1,klon
-        !
-        if(zcond(i).gt.zoliq(i)+1.e-10) then
-         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
-        else
-         beta(i,k) = 0.
-        endif
-        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
-             * (paprs(i,k)-paprs(i,k+1))/RG
-        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-           !AA lessivage nucleation LMD5 dans la couche elle-meme
-          IF (iflag_t_glace.EQ.0) THEN
-           if (t(i,k) .GE. t_glace_min_old) THEN
-              zalpha_tr = a_tr_sca(3)
-           else
-              zalpha_tr = a_tr_sca(4)
-           endif
-          ELSE ! of IF (iflag_t_glace.EQ.0)
-           if (t(i,k) .GE. t_glace_min) THEN
-              zalpha_tr = a_tr_sca(3)
-           else
-              zalpha_tr = a_tr_sca(4)
-           endif
-          ENDIF
-           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
-           !
-           ! nucleation avec un facteur -1 au lieu de -0.5
-           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
-           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-        ENDIF
-        !
-     ENDDO      ! boucle sur i
-     !
-     !AA Lessivage par impaction dans les couches en-dessous
-     DO kk = k-1, 1, -1
-        DO i = 1, klon
-           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-             IF (iflag_t_glace.EQ.0) THEN
-              if (t(i,kk) .GE. t_glace_min_old) THEN
-                 zalpha_tr = a_tr_sca(1)
-              else
-                 zalpha_tr = a_tr_sca(2)
-              endif
-             ELSE ! of IF (iflag_t_glace.EQ.0)
-              if (t(i,kk) .GE. t_glace_min) THEN
-                 zalpha_tr = a_tr_sca(1)
-              else
-                 zalpha_tr = a_tr_sca(2)
-              endif
-             ENDIF
-              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
-              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
-           ENDIF
-        ENDDO
-     ENDDO !DO kk = k-1, 1, -1
-
-#ifdef ISO      
-      ! verif en fin de boucle
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        do i=1,klon       
-        call iso_verif_egalite_choix(zxtrfl(iso_eau,i) &
-     &      ,zrfl(i), 'il pleut 1089',errmax,errmaxrel)
-        call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
-     &           'il pleut 1205',errmax,errmaxrel) 
-#ifdef ISOTRAC
-        call iso_verif_traceur(zxtrfl(1,i),'il pleut 1894')
-        call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1895')
-#endif        
-        enddo !do i=1,klon
-      endif !  if (iso_eau.gt.0) then        
-#endif
-#endif      
-      ! end verif en fin de boucle
-
-     !
-     !AA===============================================================
-     !                     FIN DE LA BOUCLE VERTICALE  
-  end DO
-  !
-  !AA==================================================================
-  !
-  ! Pluie ou neige au sol selon la temperature de la 1ere couche
-  !
-!CR: si la thermo de la glace est active, on calcule zifl directement
-  IF (.NOT.ice_thermo) THEN
-  DO i = 1, klon
-     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
-!AJ<
-!        snow(i) = zrfl(i)
-        snow(i) = zrfl(i)+zifl(i)
-!>AJ
-        zlh_solid(i) = RLSTT-RLVTT
-#ifdef ISO
-        do ixt=1,ntraciso
-          xtsnow(ixt,i) = zxtrfl(ixt,i)
-        enddo !do ixt=1,ntraciso
-#endif
-     ELSE
-        rain(i) = zrfl(i)
-        zlh_solid(i) = 0.
-#ifdef ISO
-        do ixt=1,ntraciso
-          xtrain(ixt,i) = zxtrfl(ixt,i)
-        enddo !do ixt=1,ntraciso
-#endif
-     ENDIF
-#ifdef ISO
-      ! cam verif 
-#ifdef ISOVERIF
-       do ixt=1,ntraciso
-        call iso_verif_noNAN(xtrain(ixt,i),'il pleut 927')
-        call iso_verif_noNAN(xtsnow(ixt,i),'il pleut 927')
-       enddo ! do ixt=1,niso
-#endif         
-#ifdef ISOVERIF
-       if (iso_eau.gt.0) then
-        call iso_verif_egalite_choix(xtrain(iso_eau,i),rain(i), &
-     &           'il pleut 926a',errmax,errmaxrel)  
-        if (snow(i)*dtime.gt.ridicule) then
-           call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
-     &            'il pleut 938a',errmax,errmaxrel)
-        endif !if (snow(i)*dtime.gt.ridicule) then
-       endif !if (iso_eau.gt.0) then
-       if (iso_HDO.gt.0) then
-           call iso_verif_aberrant_choix(xtrain(iso_hdo,i),rain(i), &
-     &           ridicule_rain,deltalim,'il pleut 926b')  
-           call iso_verif_aberrant_choix(xtsnow(iso_hdo,i),snow(i), &
-     &            ridicule_rain,deltalim_snow,'il pleut 938b')
-       endif !if (iso_eau.gt.0) then
-#ifdef ISOTRAC
-       call iso_verif_traceur(xtrain(1,i),'il pleut 1952')
-       call iso_verif_traceur(xtsnow(1,i),'il pleut 1953')
-#endif       
-#endif
-       ! end cam verif   
-#endif
-  ENDDO
-
-  ELSE
-#ifdef ISO
-        CALL abort_physic('ilp 2708', 'isos pas prevus ici', 1)
-#endif
-     DO i = 1, klon
-        snow(i) = zifl(i)
-        rain(i) = zrfl(i)
-     ENDDO
-   
-   ENDIF
-  !
-  ! For energy conservation : when snow is present, the solification
-  ! latent heat is considered.
-!CR: si thermo de la glace, neige deja prise en compte
-  IF (.not.ice_thermo) THEN
-  DO k = 1, klev
-     DO i = 1, klon
-        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
-        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
-        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
-        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
-     END DO
-  END DO
-  ENDIF
-  !
-
-  if (ncoreczq>0) then
-     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
-  endif
-
-
-#ifdef ISO
-      ! verif en sortie de il pleut 
-#ifdef ISOVERIF
-       do k=1,klev
-        do i=1,klon
-         if (iso_eau.gt.0) then
-           call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
-     &            'il pleut 1289',errmax,errmaxrel)
-         endif !if (iso_eau.gt.0) then
-#ifdef ISOTRAC
-         call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1982')
-#endif
-        enddo !do i=1,klon
-       enddo !do k=1,klev
-#endif
-      ! end verif
-#endif
-
-END SUBROUTINE fisrtilp
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_lscp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_lscp.F90	(revision 4664)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_lscp.F90	(revision 4664)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_lscp.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_ini.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_ini.F90	(revision 4664)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_ini.F90	(revision 4664)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_lscp_ini.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_old.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_old.F90	(revision 4664)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_old.F90	(revision 4664)
@@ -0,0 +1,2745 @@
+!
+! $Id: fisrtilp.F90 3493 2019-05-07 08:45:07Z idelkadi $
+!
+!
+MODULE lmdz_lscp_old
+CONTAINS
+SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
+     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
+     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
+     frac_impa, frac_nucl, beta,                        &
+     prfl, psfl, rhcl, zqta, fraca,                     &
+     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
+     iflag_ice_thermo &
+#ifdef ISO
+     &           ,xt,xtrain,xtsnow                      &
+     &           ,d_xt,d_xtl,d_xti,pxtrainfl,pxtsnowfl        &
+#endif      
+     &  )
+
+  !
+  USE dimphy
+  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
+  USE print_control_mod, ONLY: prt_level, lunout
+  USE lmdz_cloudth, only : cloudth, cloudth_v3, cloudth_v6
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
+  USE phys_local_var_mod, ONLY: rneblsvol
+  ! flag to include modifications to ensure energy conservation (if flag >0)
+  USE add_phys_tend_mod, only : fl_cor_ebil 
+  USE lmdz_lscp_ini, ONLY: iflag_t_glace,t_glace_min, t_glace_max, exposant_glace
+  USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol 
+  USE lmdz_lscp_ini, ONLY: coef_eva, coef_eva_i, ffallv_lsc, ffallv_con
+  USE lmdz_lscp_ini, ONLY: cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
+  USE lmdz_lscp_ini, ONLY: reevap_ice, iflag_bergeron, iflag_fisrtilp_qsat, iflag_pdf 
+      use phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv
+      use phys_output_var_mod, ONLY : cloudth_sigmath,cloudth_sigmaenv
+
+
+
+
+#ifdef ISO
+  USE infotrac_phy, ONLY: ntraciso=>ntiso,niso,itZonIso
+  USE isotopes_mod
+!, ONLY: essai_convergence,bidouille_anti_divergence, &
+!        &       Rdefault,ridicule,pxtice,pxtmelt,iso_eau,iso_HDO, &
+!        &       iso_O18,ridicule_rain
+  USE isotopes_routines_mod, ONLY: iso_revap_fisrtilp, &
+        condiso_liq_ice_vectall
+#ifdef ISOVERIF
+  USE isotopes_verif_mod, ONLY: errmax,errmaxrel,deltalim_snow,deltalim, &
+        iso_verif_egalite_choix,iso_verif_aberrant_choix, &
+        iso_verif_positif,iso_verif_positif_choix,iso_verif_noNaN, &
+        iso_verif_aberrant,iso_verif_positif_choix_nostop, &
+        iso_verif_aberrant_encadre,iso_verif_egalite_choix_nostop, &
+        iso_verif_aberrant_nostop,iso_verif_o18_aberrant_nostop, &
+        deltaD,deltaO,iso_verif_egalite
+#endif
+#ifdef ISOTRAC
+    use isotrac_mod, only: option_traceurs,option_tmin,seuil_tag_tmin_ls, &
+&       izone_poubelle, option_cond,izone_cond,index_iso, index_zone,izone_oce,izone_ddft
+    use isotrac_routines_mod, only: iso_recolorise_condensation
+    use isotopes_routines_mod, only: condiso_liq_ice_vectall_trac
+#ifdef ISOVERIF
+  USE isotopes_verif_mod, ONLY: iso_verif_traceur_nostop,iso_verif_traceur
+  USE isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille
+#endif
+#endif
+#endif
+
+  IMPLICIT none
+  !======================================================================
+  ! Auteur(s): Z.X. Li (LMD/CNRS)
+  ! Date: le 20 mars 1995
+  ! Objet: condensation et precipitation stratiforme.
+  !        schema de nuage
+  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
+  !             et ilp (il pleut, L. Li)
+  ! Principales parties:
+  ! P0> Thermalisation des precipitations venant de la couche du dessus
+  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
+  ! P2> Formation du nuage (en k) 
+  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
+  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
+  ! les valeurs de T et Q initiales
+  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
+  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
+  ! P2.B> Nuage "tout ou rien"
+  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
+  ! P3> Formation de la precipitation (en k)
+  !======================================================================
+  ! JLD: 
+  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
+  ! * fl_cor_ebil doit etre > 0 ; 
+  !   fl_cor_ebil= 0 pour reproduire anciens bugs
+  !======================================================================
+  include "YOMCST.h"
+  !
+  ! Principaux inputs:
+  !
+  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
+  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
+  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
+  INTEGER,                         INTENT(IN)   :: iflag_cld_th
+  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
+  !
+  ! Inputs lies aux thermiques
+  ! 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zthl
+  !
+  !  Input/output
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
+  !
+  ! Principaux outputs:
+  !
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb ! fraction nuageuse
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
+  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
+  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl 
+
+#ifdef ISO
+      real double_to_real
+      double precision real_to_double
+      !integer niso
+      real xt(ntraciso,klon,klev)
+      real d_xt(ntraciso,klon,klev)
+      real d_xtl(ntraciso,klon,klev)
+      real d_xti(ntraciso,klon,klev)
+      real xtrain(ntraciso,klon)
+      real xtsnow(ntraciso,klon)
+      real pxtrainfl(ntraciso,klon,klev+1)
+      real pxtsnowfl(ntraciso,klon,klev+1)      
+      ! xt <-> q
+      ! xtrain et xtsonw <-> rain et snow
+      ! d_xt et d_xtl <-> d_q et d_ql
+      ! pxtrainfl et pxtsnowflux <-> prfl et psfl
+#endif
+
+
+  !AA
+  ! Coeffients de fraction lessivee : pour OFF-LINE
+  !
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_nucl
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_1nucl
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_impa
+  !
+  ! Fraction d'aerosols lessivee par impaction et par nucleation
+  ! POur ON-LINE
+  !
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
+  !AA
+  ! --------------------------------------------------------------------------------
+  !
+  ! Options du programme:
+  !
+  REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
+  !$OMP THREADPRIVATE(seuil_neb)
+
+
+  INTEGER ninter ! sous-intervals pour la precipitation
+  PARAMETER (ninter=5)
+  INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie
+  !$OMP THREADPRIVATE(iflag_evap_prec)
+  !
+  LOGICAL cpartiel ! condensation partielle
+  PARAMETER (cpartiel=.TRUE.)
+  REAL t_coup
+  PARAMETER (t_coup=234.0)
+  REAL DDT0
+  PARAMETER (DDT0=.01)
+  REAL ztfondue
+  PARAMETER (ztfondue=278.15)
+  ! --------------------------------------------------------------------------------
+  !
+  ! Variables locales:
+  !
+  INTEGER i, k, n, kk
+  INTEGER,save::itap=0
+  !$OMP THREADPRIVATE(itap)
+
+  REAL qsl, qsi
+  real zct      ,zcl
+  INTEGER ncoreczq  
+  REAL ctot(klon,klev)
+  REAL ctot_vol(klon,klev)
+  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
+  REAL zdqsdT_raw(klon)
+  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
+
+  logical lognormale(klon)
+  logical ice_thermo
+  LOGICAL convergence(klon)
+  INTEGER n_i(klon), iter
+  REAL cste
+
+  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
+  real erf   
+  REAL qcloud(klon)
+ 
+  REAL zrfl(klon), zrfln(klon), zqev, zqevt 
+  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
+  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq 
+  REAL zoliqp(klon), zoliqi(klon)
+  REAL zt(klon)
+! JBM (3/14) nexpo is replaced by exposant_glace
+! REAL nexpo ! exponentiel pour glace/eau
+! INTEGER, PARAMETER :: nexpo=6
+  INTEGER exposant_glace_old
+  REAL t_glace_min_old
+  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
+  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
+  REAL zmelt, zpluie, zice
+  REAL dzfice(klon)
+  REAL zsolid
+!!!!
+!  Variables pour Bergeron
+  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
+  REAL zqpreci(klon), zqprecl(klon)
+! Variable pour conservation enegie des precipitations
+  REAL zmqc(klon)
+
+#ifdef ISO
+      real zxtrfl(ntraciso,klon)
+      real zxtrfln(ntraciso,klon)
+      REAL zxtifl(ntraciso,klon), zxtifln(ntraciso,klon)  
+      real zxt(ntraciso,klon) ! equivalent local de xt
+      real zxtn(ntraciso,klon) ! isotpes nuageux
+      real zxtn_reduit(niso)
+      real zxtcond(ntraciso,klon)
+      real zxtoliq(ntraciso,klon)
+      real zq_ancien(klon)
+      REAL zxtpreci(ntraciso,klon), zxtprecl(ntraciso,klon)
+      integer ixt
+#ifdef ISOVERIF   
+!      integer iso_verif_aberrant_nostop
+!      integer iso_verif_egalite_choix_nostop
+!      integer iso_verif_positif_nostop
+!      integer iso_verif_positif_choix_nostop
+      integer trace_cas(klon)
+!      integer iso_verif_traceur_nostop
+#endif      
+      integer iso_verif_nonan_nostop
+      real zxtice(ntraciso,klon),zxtliq(ntraciso,klon)
+      real zrfl_ancien(klon)
+      real zqev_diag(klon)
+      real zxtrfl_ancien(ntraciso,klon)
+      real zxt_ancien(ntraciso,klon)
+#ifdef ISOTRAC
+        real zxtres(ntraciso) ! humidite spec nuageuse après condensation et precip
+        real zqcs(klon) ! humidite spec en ciel clair
+        real zxtcs(ntraciso,klon)        
+        real zcondn(klon) ! humidite spec de l'eau liquide dans le nuage
+        real zxtcondn(ntraciso,klon)
+        integer izone
+#ifdef ISOVERIF           
+        real zq_verif
+#endif        
+#endif 
+!      logical negation ! pour nier conditions logiques
+#endif
+
+  !
+  LOGICAL appel1er
+  SAVE appel1er
+  !$OMP THREADPRIVATE(appel1er)
+  !
+! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
+! iflag_oldbug_fisrtilp=1 ajoute le BUG
+  INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug
+  !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp)
+  !---------------------------------------------------------------
+  !
+  !AA Variables traceurs:
+  !AA  Provisoire !!! Parametres alpha du lessivage
+  !AA  A priori on a 4 scavenging # possibles
+  !
+  REAL a_tr_sca(4)
+  save a_tr_sca
+  !$OMP THREADPRIVATE(a_tr_sca)
+  !
+  ! Variables intermediaires
+  !
+  REAL zalpha_tr
+  REAL zfrac_lessi
+  REAL zprec_cond(klon)
+  !AA
+! RomP >>> 15 nov 2012
+  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
+! RomP <<<
+  REAL zmair(klon), zcpair, zcpeau
+  !     Pour la conversion eau-neige
+  REAL zlh_solid(klon), zm_solid
+  !---------------------------------------------------------------
+  !
+  ! Fonctions en ligne:
+  !
+  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
+                     ! (Heymsfield & Donner, 1990)
+  REAL zzz
+
+  include "YOETHF.h"
+  include "FCTTRE.h"
+  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
+  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
+  !
+  DATA appel1er /.TRUE./
+  !ym
+!CR: pour iflag_ice_thermo=2, on active que la convection
+!  ice_thermo = iflag_ice_thermo .GE. 1
+
+  itap=itap+1
+  znebprecip(:)=0.
+
+  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
+  zdelq=0.0
+  ctot_vol(1:klon,1:klev)=0.0
+  rneblsvol(1:klon,1:klev)=0.0
+
+  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
+  IF (appel1er) THEN
+     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
+     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
+     CALL getin_p('seuil_neb',seuil_neb)
+     write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp 
+     !
+     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
+     WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
+     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
+     
+     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
+        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
+        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
+        !         CALL abort
+     ENDIF
+     appel1er = .FALSE.
+
+#ifdef ISO
+        ! cam verif 
+#ifdef ISOVERIF
+        write(*,*) 'ilp 310'
+          do k=1,klev
+           do i=1,klon             
+            if (iso_eau.gt.0) then        
+             call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &   
+     &                   'il pleut 205',errmax,errmaxrel)
+            endif !if (iso_eau.gt.0) then
+            if (iso_HDO.gt.0) then    
+             if (q(i,k).gt.ridicule) then
+               call iso_verif_aberrant(xt(iso_HDO,i,k)/q(i,k), &
+     &           'il pleut 210')
+             endif !if (zq(i).gt.ridicule) then
+            endif  !if (iso_HDO.gt.0) then
+           enddo !do k=1,klev
+          enddo !do i=1,klon
+        do k=1,klev
+           do i=1,klon
+             call iso_verif_positif(370.0-t(i,k),'il pleut 250')
+             call iso_verif_positif(t(i,k)-100.0,'il pleut 251')
+           enddo
+        enddo
+        write(*,*) 'ilp 331'
+#endif
+         ! end cam verif
+#endif 
+
+     !
+     !AA initialiation provisoire
+     a_tr_sca(1) = -0.5
+     a_tr_sca(2) = -0.5
+     a_tr_sca(3) = -0.5
+     a_tr_sca(4) = -0.5
+     !
+     !AA Initialisation a 1 des coefs des fractions lessivees 
+     !
+     !cdir collapse
+     DO k = 1, klev
+        DO i = 1, klon
+           pfrac_nucl(i,k)=1.
+           pfrac_1nucl(i,k)=1.
+           pfrac_impa(i,k)=1.
+           beta(i,k)=0.  !RomP initialisation
+        ENDDO
+     ENDDO
+
+  ENDIF          !  test sur appel1er
+  !
+  !MAf Initialisation a 0 de zoliq
+  !      DO i = 1, klon
+  !         zoliq(i)=0.
+  !      ENDDO 
+  ! Determiner les nuages froids par leur temperature
+  !  nexpo regle la raideur de la transition eau liquide / eau glace.
+  !
+!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
+  IF (iflag_t_glace.EQ.0) THEN
+!   ztglace = RTT - 15.0
+    t_glace_min_old = RTT - 15.0
+    !AJ<
+    IF (ice_thermo) THEN
+!     nexpo = 2
+      exposant_glace_old = 2
+    ELSE
+!     nexpo = 6
+      exposant_glace_old = 6
+    ENDIF
+    
+  ENDIF
+  
+!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
+!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
+!>AJ
+  !cc      nexpo = 1
+  !
+  ! Initialiser les sorties:
+  !
+  !cdir collapse
+  DO k = 1, klev+1
+     DO i = 1, klon
+        prfl(i,k) = 0.0
+        psfl(i,k) = 0.0
+#ifdef ISO                 
+          do ixt=1,ntraciso
+            pxtrainfl(ixt,i,k)=0.0
+            pxtsnowfl(ixt,i,k)=0.0
+          enddo !do ixt=1,niso
+#endif
+     ENDDO
+  ENDDO
+
+  !cdir collapse
+  DO k = 1, klev
+     DO i = 1, klon
+        d_t(i,k) = 0.0
+        d_q(i,k) = 0.0
+        d_ql(i,k) = 0.0
+        d_qi(i,k) = 0.0
+#ifdef ISO
+          do ixt=1,ntraciso
+            d_xt(ixt,i,k) = 0.0
+            d_xtl(ixt,i,k) = 0.0
+            d_xti(ixt,i,k) = 0.0
+          enddo !do ixt=1,niso
+#endif 
+        rneb(i,k) = 0.0
+        radliq(i,k) = 0.0
+        frac_nucl(i,k) = 1. 
+        frac_impa(i,k) = 1. 
+     ENDDO
+  ENDDO
+  DO i = 1, klon
+     rain(i) = 0.0
+     snow(i) = 0.0
+     zoliq(i)=0.
+#ifdef ISO
+         do ixt=1,ntraciso
+            xtrain(ixt,i) = 0.0
+            xtsnow(ixt,i) = 0.0
+            zxtoliq(ixt,i)=0.
+         enddo !do ixt=1,niso
+#endif
+     !     ENDDO
+     !
+     ! Initialiser le flux de precipitation a zero
+     !
+     !     DO i = 1, klon
+     zrfl(i) = 0.0
+     zifl(i) = 0.0
+     zrfln(i) = 0.0 ! C Risi pour plus de securite
+     zifln(i) = 0.0 ! C Risi pour plus de securite
+#ifdef ISO
+          do ixt=1,ntraciso
+            zxtrfl(ixt,i)=0.0
+            zxtifl(ixt,i)=0.0
+            zxtrfln(ixt,i)=0.0
+            zxtifln(ixt,i)=0.0
+          enddo
+#endif
+     zneb(i) = seuil_neb
+  ENDDO
+  !
+  !
+  !AA Pour plus de securite 
+
+  zalpha_tr   = 0.
+  zfrac_lessi = 0.
+
+  !AA==================================================================
+  !
+  ncoreczq=0
+  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
+  !
+  DO k = klev, 1, -1
+     !
+     !
+#ifdef ISO      
+      ! cam debug!      
+       ! verif en debut de boucle
+#ifdef ISOVERIF
+       !write(*,*) 'ilp 478: boucle en k: k=',k         
+       do i=1,klon
+        if (iso_eau.gt.0) then
+        call iso_verif_egalite_choix(zxtrfl(iso_eau,i), &
+     &           zrfl(i), &
+     &           'il pleut 310',errmax,errmaxrel)
+        endif !if (iso_eau.gt.0) then
+        if (iso_HDO.gt.0) then  
+             call iso_verif_aberrant_choix(zxtrfl(iso_HDO,i) &
+     &           ,zrfl(i),ridicule_rain,deltalim_snow,'il pleut 316')
+        endif  !if (iso_HDO.gt.0) then  
+#ifdef ISOTRAC   
+        call iso_verif_traceur(zxtrfl(1,i), &
+     &          'il pleut 365')   
+        call iso_verif_traceur_pbidouille(zxtrfl(1,i), &
+     &          'il pleut 388')
+#endif         
+       enddo !do i=1,klon
+#endif
+#ifdef ISOVERIF   
+        !write(*,*) 'ilp tmp 412'      
+       do i=1,klon
+         do ixt=1,ntraciso
+           call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 406')
+           call iso_verif_noNaN(zxtoliq(ixt,i),'ilp 407')
+         enddo !do ixt=1,ntraciso
+       enddo !do i=1,klon
+#endif
+      ! end verif en debut de boucle
+      ! end cam debug
+#endif
+     !AA===============================================================
+     !
+     ! Initialisation temperature et vapeur
+     DO i = 1, klon
+        zt(i)=t(i,k)
+        zq(i)=q(i,k)
+#ifdef ISOVERIF
+        if (q(i,k).lt.0.0) then
+               write(*,*) 'ilp 400: avant clmain 31janv: i,k,q=', &
+     &           i,k,q(i,k)
+               CALL abort_physic('ilp 484', 'on stop', 1)
+        endif 
+#endif
+#ifdef ISO
+           zq(i)=max(0.0,zq(i))
+#ifdef ISOTRAC
+#ifdef ISOVERIF         
+          call iso_verif_traceur(xt(1,i,k),'ilp 404')           
+          call iso_verif_traceur_pbidouille(xt(1,i,k),'ilp 407')
+#endif          
+#endif          
+         do ixt=1,ntraciso
+           zxt(ixt,i)=xt(ixt,i,k)
+           zxt(ixt,i)=max(0.0,zxt(ixt,i))
+         enddo !do ixt=1,niso  
+
+         ! cam verif
+#ifdef ISOVERIF
+         !write(*,*) 'ilp tmp 562'
+             call iso_verif_positif(370.0-zt(i),'il pleut 425')
+             call iso_verif_positif(zt(i)-100.0,'il pleut 426')
+          if (iso_eau.gt.0) then        
+            call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
+     &                   'il pleut 289',errmax,errmaxrel)
+          endif !if (iso_eau.gt.0) then
+          if (iso_HDO.gt.0) then    
+             if (zq(i).gt.ridicule) then
+               call iso_verif_aberrant(zxt(iso_HDO,i)/zq(i), &
+     &           'il pleut 337')
+             endif !if (zq(i).gt.ridicule) then
+          endif  !if (iso_HDO.gt.0) then  
+#ifdef ISOTRAC
+          call iso_verif_traceur(zxt(1,i),'ilp 431')
+#endif           
+#endif
+         ! end cam verif
+#endif  
+
+     ENDDO
+     !
+     ! ----------------------------------------------------------------
+     ! P0> Thermalisation des precipitations venant de la couche du dessus
+     ! ----------------------------------------------------------------
+     ! Calculer la varition de temp. de l'air du a la chaleur sensible
+     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
+     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
+     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
+     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
+     ! de l'enthalpie  de la couche avec la temperature
+     ! Variables calculees ou modifiees:
+     !   -  zt: temperature de la cocuhe
+     !   - zmqc: masse de precip qui doit etre thermalisee
+     !
+     IF(k.LE.klevm1) THEN         
+        DO i = 1, klon
+           !IM
+           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
+           zcpair=RCPD*(1.0+RVTMP2*zq(i))
+           zcpeau=RCPD*RVTMP2
+         if (fl_cor_ebil .GT. 0) then
+           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm 
+           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la 
+           ! derniere couche
+           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
+           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
+           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+                 / (zcpair + zmqc(i)*zcpeau)
+         else ! si on maintient les anciennes erreurs
+           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
+                + zmair(i)*zcpair*zt(i) ) &
+                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
+         end if
+        ENDDO
+     ELSE  ! IF(k.LE.klevm1) 
+        DO i = 1, klon 
+           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 
+           zmqc(i) = 0.
+        ENDDO
+     ENDIF ! end IF(k.LE.klevm1)
+     !
+     ! ----------------------------------------------------------------
+     ! P1> Calcul de l'evaporation de la precipitation
+     ! ----------------------------------------------------------------
+     ! On evapore une partie des precipitations venant de la maille du dessus.
+     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
+     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
+     ! Variables calculees ou modifiees:
+     !   - zrfl et zifl: flux de precip liquide et glace
+     !   - zq, zt: humidite et temperature de la cocuhe
+     !   - zmqc: masse de precip qui doit etre thermalisee
+     !
+#ifdef ISO
+      DO i = 1, klon
+         zq_ancien(i)=zq(i)
+         zqev_diag(i)=0.0
+         zrfl_ancien(i)=zrfl(i)
+         do ixt=1,ntraciso
+           zxtrfl_ancien(ixt,i)=zxtrfl(ixt,i)
+           zxt_ancien(ixt,i)=zxt(ixt,i)
+         enddo 
+      enddo   ! DO i = 1, klon  
+#ifdef ISOVERIF
+        !write(*,*) 'ilp tmp 616'
+      if (iso_eau.gt.0) then
+       DO i = 1, klon           
+          call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
+     &           'il pleut 426a',errmax,errmaxrel) 
+          call iso_verif_egalite_choix(zxtrfl_ancien(iso_eau,i), &
+     &            zrfl_ancien(i), &
+     &           'il pleut 426b',errmax,errmaxrel)
+        enddo
+      endif
+#ifdef ISOTRAC
+       do i=1,klon      
+            call iso_verif_traceur(zxt_ancien(1,i),'ilp 532')
+            call iso_verif_traceur(zxtrfl_ancien(1,i), &
+     &           'ilp 533')  
+       enddo !do i=1,klon  
+#endif      
+#endif      
+#endif
+
+
+     IF (iflag_evap_prec>=1) THEN
+        DO i = 1, klon
+!          S'il y a des precipitations
+           IF (zrfl(i)+zifl(i).GT.0.) THEN
+              ! Calcul du qsat
+              IF (thermcep) THEN
+                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
+                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+                 zqs(i)=MIN(0.5,zqs(i))
+                 zcor=1./(1.-RETV*zqs(i))
+                 zqs(i)=zqs(i)*zcor
+              ELSE
+                 IF (zt(i) .LT. t_coup) THEN
+                    zqs(i) = qsats(zt(i)) / pplay(i,k)
+                 ELSE
+                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
+                 ENDIF
+              ENDIF
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+        ENDDO
+!AJ<
+
+       IF (.NOT. ice_thermo) THEN
+        DO i = 1, klon
+!          S'il y a des precipitations
+           IF (zrfl(i)+zifl(i).GT.0.) THEN                
+                ! Evap max pour ne pas saturer la fraction sous le nuage
+                ! Evap max jusqu'à atteindre la saturation dans la partie
+                ! de la maille qui est sous le nuage de la couche du dessus
+                !!! On ne tient compte de cette fraction que sous une seule
+                !!! couche sous le nuage
+                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+             ! Ajout de la prise en compte des precip a thermiser
+             ! avec petite reecriture 
+             if  (fl_cor_ebil .GT. 0) then ! nouveau
+                ! Calcul de l'evaporation du flux de precip herite
+                !   d'au-dessus
+                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
+                     * zmair(i)/pplay(i,k)*zt(i)*RD
+                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
+                ! Seuil pour ne pas saturer la fraction sous le nuage
+                zqev = MIN (zqev, zqevt)
+                ! Nouveau flux de precip
+                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
+                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
+                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
+                  zrfln(i)=0.
+                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
+                END IF
+                ! Nouvelle vapeur 
+                zq(i) = zq(i) + zqev
+                zmqc(i) = zmqc(i)-zqev
+                ! Nouvelle temperature (chaleur latente)
+                zt(i) = zt(i) - zqev &
+                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+!!JLD debut de partie a supprimer a terme
+            else ! if  (fl_cor_ebil .GT. 0)
+                ! Calcul de l'evaporation du flux de precip herite
+                !   d'au-dessus
+                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
+                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
+                ! Seuil pour ne pas saturer la fraction sous le nuage
+                zqev = MIN (zqev, zqevt)
+                ! Nouveau flux de precip
+                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
+                     /RG/dtime
+                ! Aucun flux liquide pour T < t_coup
+                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
+                ! Nouvelle vapeur 
+                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
+                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                ! Nouvelle temperature (chaleur latente)
+                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+              end if ! if  (fl_cor_ebil .GT. 0)
+
+!!JLD fin de partie a supprimer a terme
+                zrfl(i) = zrfln(i)
+                zifl(i) = 0.
+#ifdef ISO
+         zqev_diag(i)=zqev
+#endif
+
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+
+
+        ENDDO
+#ifdef ISO
+
+      ! * traitement de l'évaporation
+      call iso_revap_fisrtilp(klon,klev,k, &
+     &            zrfl_ancien,zrfl,zrfln,zt,zxt_ancien, &
+     &            zxtrfl,zxtrfl_ancien,zxtrfln,zxt, &
+     &            paprs,dtime, &
+     &            zqs,zq_ancien,zqev_diag,zq)  
+#endif
+
+!
+       ELSE ! (.NOT. ice_thermo) 
+
+#ifdef ISO
+        CALL abort_physic('ilp 664', 'isotopes pas encore prevus dans ice_thermo', 1)        
+#endif
+!      ================================
+!      Avec thermodynamique de la glace
+!      ================================
+        DO i = 1, klon
+!AJ<
+!        S'il y a des precipitations
+         IF (zrfl(i)+zifl(i).GT.0.) THEN
+
+         IF (iflag_evap_prec==1) THEN
+            znebprecip(i)=zneb(i)
+         ELSE
+            znebprecip(i)=MAX(zneb(i),znebprecip(i))
+         ENDIF
+     
+        ! Evap max pour ne pas saturer la fraction sous le nuage
+         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
+
+         !JAM
+         ! On differencie qsat pour l'eau et la glace
+         ! Si zdelta=1. --> glace
+         ! Si zdelta=0. --> eau liquide
+        
+         ! Calcul du qsat par rapport a l'eau liquide
+         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
+         qsl= MIN(0.5,qsl)
+         zcor= 1./(1.-RETV*qsl)
+         qsl= qsl*zcor
+         
+         ! Calcul de l'evaporation du flux de precip venant du dessus
+         ! Formulation en racine du flux de precip
+         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
+         IF (iflag_evap_prec==3) THEN
+         zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
+              *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+         ELSE
+         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+         ENDIF
+
+
+         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+              *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
+         
+         ! Calcul du qsat par rapport a la glace
+         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
+         qsi= MIN(0.5,qsi)
+         zcor= 1./(1.-RETV*qsi)
+         qsi= qsi*zcor
+
+         ! Calcul de la sublimation du flux de precip solide herite
+         !   d'au-dessus
+         IF (iflag_evap_prec==3) THEN
+         zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
+              *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+         ELSE
+         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+         ENDIF
+         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
+              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
+
+        !JAM
+        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas 
+        ! la fraction de la couche sous le nuage sinon on repartit zqev0 
+        ! en conservant la proportion liquide / glace
+     
+         IF (zqevt+zqevti.GT.zqev0) THEN
+            zqev=zqev0*zqevt/(zqevt+zqevti)
+            zqevi=zqev0*zqevti/(zqevt+zqevti)
+         ELSE
+!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
+!       liquides et solides meme si on ne sature pas la couche. 
+!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
+!            zqev=zqevt
+!            zqevi=zqevti
+             IF (zqevt+zqevti.GT.0.) THEN
+                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
+                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
+             ELSE
+             zqev=0.
+             zqevi=0.
+             ENDIF
+         ENDIF
+         ! Nouveaux flux de precip liquide et solide 
+         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
+                                 /RG/dtime)
+         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
+                                 /RG/dtime)
+
+         ! Mise a jour de la vapeur, temperature et flux de precip
+         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
+         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
+                  + (zifln(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+       else ! sans correction thermalisation des precips
+         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
+                  + (zifln(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
+       end if
+        ! Nouvelles vaeleurs des precips liquides et solides 
+         zrfl(i) = zrfln(i)
+         zifl(i) = zifln(i)
+!        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
+
+!CR ATTENTION: deplacement de la fonte de la glace
+!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
+!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
+!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
+           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
+           ! fraction de la precip solide qui est fondue
+           zmelt = MIN(MAX(zmelt,0.),1.)
+           ! Fusion de la glace
+           zrfl(i)=zrfl(i)+zmelt*zifl(i)
+           if (fl_cor_ebil .LE. 0) then
+             ! the following line should not be here. Indeed, if zifl is modified
+             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
+             ! and therefore the change in temperature computed below is wrong
+             zifl(i)=zifl(i)*(1.-zmelt)
+           end if
+           ! Chaleur latente de fusion
+        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
+           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+        else ! sans correction thermalisation des precips
+           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
+        end if
+           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
+             zifl(i)=zifl(i)*(1.-zmelt)
+           end if
+
+           ELSE
+              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
+              ! sous nuageuse utilisee pour la pluie.
+              znebprecip(i)=0.
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+        ENDDO
+
+      ENDIF ! (.NOT. ice_thermo) 
+     
+     ! ----------------------------------------------------------------
+     ! Fin evaporation de la precipitation
+     ! ----------------------------------------------------------------
+     ENDIF ! (iflag_evap_prec>=1)
+     !
+     ! Calculer Qs et L/Cp*dQs/dT:
+#ifdef ISOVERIF
+!        write(*,*) 'ilp tmp 849'
+           do i=1,klon
+             call iso_verif_positif(370.0-zt(i),'il pleut 938')
+             call iso_verif_positif(zt(i)-100.0,'il pleut 939')
+             do ixt=1,ntraciso
+              call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 614')
+             enddo   
+             if (zq(i).gt.ridicule) then
+               if (iso_HDO.gt.0) then
+                  call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
+     &                   'il pleut 856')                 
+                  if (iso_O18.gt.0) then                
+                   if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
+     &               zxt(iso_O18,i)/zq(i),'il pleut 859').eq.1) then
+                     write(*,*) 'i,k,zq=',i,k,zq(i)
+                     write(*,*) 'zqev_diag(i)=',zqev_diag(i)
+                     write(*,*) 'zrfl_ancien(i)=',zrfl_ancien(i)
+                     write(*,*) 'zrfln(i)=',zrfln(i)
+                     write(*,*) 'zq_ancien(i)=',zq_ancien(i)
+                     write(*,*) 'deltaD(zrfl_ancien)=',deltaD(zxtrfl_ancien(iso_HDO,i))
+                     write(*,*) 'deltaO(zrfl_ancien)=',deltaO(zxtrfl_ancien(iso_O18,i))
+                     write(*,*) 'deltaD(zq_ancien)=',deltaD(zxt_ancien(iso_HDO,i))
+                     write(*,*) 'deltaO(zq_ancien)=',deltaO(zxt_ancien(iso_O18,i))
+                     write(*,*) 'zt=',zt(i)
+                     stop
+                endif !if (iso_verif_O18_aberrant_nostop
+               endif !if (iso_O18.gt.0) then
+              endif !if (iso_HDO.gt.0) then
+             endif !if (zq(i).gt.ridicule) then                
+           enddo !do i=1,klon
+#ifdef ISOTRAC
+!      write(*,*) 'ilp tmp 633'
+           do i=1,klon
+              call iso_verif_traceur(zxtrfl(1,i),'il pleut 632')
+           enddo
+!      write(*,*) 'ilp tmp 637,thermcep'
+#endif           
+#endif    
+#ifdef ISOVERIF
+!        write(*,*) 'ilp tmp 888'
+           do i=1,klon
+            do ixt=1,ntraciso
+             call iso_verif_noNaN(zxtrfl(ixt,i),'ilp 624')
+            enddo !do ixt=1,ntraciso
+           enddo ! do i=1,klon
+!        write(*,*) 'ilp tmp 638: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
+#endif
+     !
+     IF (thermcep) THEN
+        DO i = 1, klon
+           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
+           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+       if  (fl_cor_ebil .GT. 0) then ! nouveau
+           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+       else    
+           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+       endif
+           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+           zqs(i) = MIN(0.5,zqs(i))
+           zcor = 1./(1.-RETV*zqs(i))
+           zqs(i) = zqs(i)*zcor
+           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
+           zdqsdT_raw(i) = zdqs(i)*  &
+         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
+        ENDDO
+     ELSE
+        DO i = 1, klon
+           IF (zt(i).LT.t_coup) THEN
+              zqs(i) = qsats(zt(i))/pplay(i,k)
+              zdqs(i) = dqsats(zt(i),zqs(i))
+           ELSE
+              zqs(i) = qsatl(zt(i))/pplay(i,k)
+              zdqs(i) = dqsatl(zt(i),zqs(i))
+           ENDIF
+        ENDDO
+     ENDIF
+     !
+     ! Determiner la condensation partielle et calculer la quantite
+     ! de l'eau condensee:
+     !
+!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
+!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
+!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
+!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
+!         stop
+!       endif
+
+     ! ----------------------------------------------------------------
+     ! P2> Formation du nuage
+     ! ----------------------------------------------------------------
+     ! Variables calculees:
+     !   rneb  : fraction nuageuse
+     !   zcond : eau condensee moyenne dans la maille.
+     !   rhcl: humidite relative ciel-clair
+     !   zt : temperature de la maille
+     ! ----------------------------------------------------------------
+     !
+     IF (cpartiel) THEN
+        ! -------------------------
+        ! P2.A> Nuage fractionnaire
+        ! -------------------------
+        !
+        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
+        !   nuageuse a partir des PDF de Sandrine Bony.
+        !   rneb  : fraction nuageuse
+        !   zqn   : eau totale dans le nuage
+        !   zcond : eau condensee moyenne dans la maille.
+        !  on prend en compte le réchauffement qui diminue la partie
+        ! condensee
+        !
+        !   Version avec les raqts
+
+        ! ----------------------------------------------------------------
+        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
+        ! ----------------------------------------------------------------
+
+#ifdef ISO          
+#ifdef ISOVERIF
+        do i=1,klon
+          if (iso_eau.gt.0) then    
+            if (iso_verif_egalite_choix_nostop(zxt(iso_eau,i),zq(i), &
+     &            'il pleut 608',errmax,errmaxrel).eq.1) then
+              write(*,*) 'i=',i
+              stop
+            endif
+          endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+          call iso_verif_traceur(zxt(1,i),'il pleut 1106')
+#endif          
+        enddo !do i=1,klon        
+#endif
+#endif 
+
+        if (iflag_pdf.eq.0) then
+
+           ! version creneau de (Li, 1998)
+           do i=1,klon
+              zdelq = min(ratqs(i,k),0.99) * zq(i)
+              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
+              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
+           enddo
+
+        else !  if (iflag_pdf.eq.0)
+           ! ----------------------------------------------------------------
+           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
+           ! les valeurs de T et Q initiales
+           ! ----------------------------------------------------------------
+           do i=1,klon
+              if(zq(i).lt.1.e-15) then
+                 ncoreczq=ncoreczq+1
+#ifdef ISO
+                zq_ancien(i)=zq(i)
+#endif
+                 zq(i)=1.e-15
+#ifdef ISO
+                ! on conserve le même rapport pour les isotopes
+                if (zq_ancien(i).gt.0.0) then
+                   do ixt=1,ntraciso
+                     zxt(ixt,i)=zxt(ixt,i)/zq_ancien(i)*zq(i)
+                     zxt(ixt,i)=max(0.0,zxt(ixt,i))
+                   enddo   
+                else !if (zq_ancien(i).gt.0.0) then
+                    do ixt=1,ntraciso
+                     zxt(ixt,i)=0.0
+                   enddo 
+                endif !if (zq_ancien(i).gt.0.0) then
+
+                ! verif
+#ifdef ISOVERIF
+                  if (iso_eau.gt.0) then
+                    call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
+     &                 'il pleut 654',errmax,errmaxrel)
+                  endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+             call iso_verif_traceur(zxt(1,i),'il pleut 1207')
+#endif                    
+#endif
+                ! end verif
+#endif 
+              endif
+           enddo
+
+           if (iflag_cld_th>=5) then
+
+              if (iflag_cloudth_vert<=2) then
+               call cloudth(klon,klev,k,ztv, &
+                   zq,zqta,fraca, &
+                   qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
+                   ratqs,zqs,t, &
+                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+              elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
+               call cloudth_v3(klon,klev,k,ztv, &
+                   zq,zqta,fraca, &
+                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+                   ratqs,zqs,t, &
+                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+              !----------------------------------
+              !Version these Jean Jouhaud, Decembre 2018 
+              !----------------------------------              
+              elseif (iflag_cloudth_vert==6) then
+               call cloudth_v6(klon,klev,k,ztv, &
+                   zq,zqta,fraca, &
+                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
+                   ratqs,zqs,t, &
+                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+
+
+              endif
+              do i=1,klon
+                 rneb(i,k)=ctot(i,k)
+                 rneblsvol(i,k)=ctot_vol(i,k)
+                 zqn(i)=qcloud(i)
+              enddo
+
+           endif
+
+           if (iflag_cld_th <= 4) then
+              lognormale = .true.
+           elseif (iflag_cld_th >= 6) then
+              ! lognormale en l'absence des thermiques
+              lognormale = fraca(:,k) < 1e-10
+           else
+              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
+              ! bi-gaussienne
+              lognormale = .false.
+           end if
+
+!CR: variation de qsat avec T en presence de glace ou non
+!initialisations
+           do i=1,klon
+              DT(i) = 0.
+              n_i(i)=0
+              Tbef(i)=zt(i)
+              qlbef(i)=0.
+           enddo
+
+        ! ----------------------------------------------------------------
+        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
+        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
+        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
+        ! qui en dependent. Ce changement de temperature est provisoire, et
+        ! la valeur definitive sera calcule plus tard.
+        ! Variables calculees: 
+        !   rneb : nebulosite
+        !   zcond: eau condensee en moyenne dans la maille
+        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
+        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
+        ! T sur qsat
+        ! ----------------------------------------------------------------
+
+!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
+           if (iflag_fisrtilp_qsat.ge.0) then
+             ! Iteration pour condensation avec variation de qsat(T) 
+             ! -----------------------------------------------------
+             do iter=1,iflag_fisrtilp_qsat+1
+                
+               do i=1,klon
+!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
+                 ! !! convergence = .true. tant que l'on n'a pas converge !!
+                 ! ------------------------------
+                 convergence(i)=abs(DT(i)).gt.DDT0
+                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
+                 ! si on n'a pas converge
+                 !
+                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
+                 ! ---------------------------------------------------------------
+                 ! Variables calculees: 
+                 ! rneb : nebulosite
+                 ! zqn : eau condensee, dans le nuage (in cloud water content)
+                 ! zcond: eau condensee en moyenne dans la maille
+                 ! rhcl: humidite relative ciel-clair
+                 !
+                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
+                 if (.not.ice_thermo) then   
+                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
+                 else
+                    if (iflag_t_glace.eq.0) then
+                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
+                    else if (iflag_t_glace.ge.1) then
+                       if (iflag_oldbug_fisrtilp.EQ.0) then
+                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
+                       else
+                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
+                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
+                       endif
+                    endif
+                 endif
+                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
+                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+               if (fl_cor_ebil .GT. 0) then
+                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+               else
+                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+               end if
+                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
+                 zqs(i) = MIN(0.5,zqs(i))
+                 zcor = 1./(1.-RETV*zqs(i))
+                 zqs(i) = zqs(i)*zcor
+                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
+                 zpdf_sig(i)=ratqs(i,k)*zq(i)
+                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
+                 zpdf_delta(i)=log(zq(i)/zqs(i))
+                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
+                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
+                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
+                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
+                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
+                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
+                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
+                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
+             
+                 if (zpdf_e1(i).lt.1.e-10) then
+                    rneb(i,k)=0.
+                    zqn(i)=zqs(i)
+                 else
+                    rneb(i,k)=0.5*zpdf_e1(i)
+                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                 endif
+
+                 ! If vertical heterogeneity, change fraction by volume as well
+                 if (iflag_cloudth_vert>=3) then
+                   ctot_vol(i,k)=rneb(i,k)
+                   rneblsvol(i,k)=ctot_vol(i,k)
+                 endif
+
+                 endif !convergence
+
+               enddo ! boucle en i
+
+                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT 
+                 !         due a la condensation.
+                 ! ---------------------------------------------------------------
+                 ! Variables calculees: 
+                 ! DT : variation de temperature due a la condensation 
+
+                 if (.not. ice_thermo) then
+                 ! --------------------------
+                 do i=1,klon
+                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
+
+                 qlbef(i)=max(0.,zqn(i)-zqs(i))
+               if (fl_cor_ebil .GT. 0) then
+                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
+               else
+                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+               end if
+                 denom=1.+rneb(i,k)*zdqs(i)
+                 DT(i)=num/denom
+                 n_i(i)=n_i(i)+1
+                 endif
+                 enddo
+
+                 else ! if (.not. ice_thermo)
+                 ! --------------------------
+                 if (iflag_t_glace.ge.1) then
+                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
+                 endif
+
+                 do i=1,klon
+                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
+                 
+                 if (iflag_t_glace.eq.0) then
+                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
+                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+                    zfice(i) = zfice(i)**exposant_glace_old
+                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
+          &                     / (t_glace_min_old - RTT)
+                 endif
+                 
+                 if (iflag_t_glace.ge.1.and.zfice(i)>0.) then
+                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
+          &                    / (t_glace_min - t_glace_max)
+                 endif
+               
+                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
+                    dzfice(i)=0.
+                 endif
+
+                 if (zfice(i).lt.1) then
+                    cste=RLVTT
+                 else
+                    cste=RLSTT
+                 endif
+
+                 qlbef(i)=max(0.,zqn(i)-zqs(i))
+               if (fl_cor_ebil .GT. 0) then
+                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
+           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
+                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
+                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
+           &               *qlbef(i)*dzfice(i)
+               else
+                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
+           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
+                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
+               end if
+                 DT(i)=num/denom
+                 n_i(i)=n_i(i)+1
+
+                 endif ! fin convergence
+                 enddo ! fin boucle i
+
+                 endif !ice_thermo
+
+             enddo ! iter=1,iflag_fisrtilp_qsat+1
+             ! Fin d'iteration pour condensation avec variation de qsat(T) 
+             ! -----------------------------------------------------------
+           endif !  if (iflag_fisrtilp_qsat.ge.0)
+     ! ----------------------------------------------------------------
+     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
+     ! ----------------------------------------------------------------
+
+        endif ! iflag_pdf
+
+!        if (iflag_fisrtilp_qsat.eq.-1) then 
+!------------------------------------------
+!CR: ATTENTION option fausse mais a existe:
+! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
+! activer les lignes suivantes: 
+       IF (1.eq.0) THEN 
+       DO i=1,klon
+           IF (rneb(i,k) .LE. 0.0) THEN
+              zqn(i) = 0.0
+              rneb(i,k) = 0.0
+              zcond(i) = 0.0
+              rhcl(i,k)=zq(i)/zqs(i)
+           ELSE IF (rneb(i,k) .GE. 1.0) THEN
+              zqn(i) = zq(i)
+              rneb(i,k) = 1.0                 
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
+              rhcl(i,k)=1.0
+           ELSE
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
+              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+           ENDIF
+       ENDDO
+       ENDIF
+!------------------------------------------
+
+!        ELSE
+        ! ----------------------------------------------------------------
+        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
+        ! Variables calculees: 
+        !   rneb : nebulosite
+        !   zcond: eau condensee en moyenne dans la maille
+        !   zq : eau vapeur dans la maille
+        !   zt : temperature de la maille
+        !   rhcl: humidite relative ciel-clair
+        ! ----------------------------------------------------------------
+        !
+        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb) 
+        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
+        ! et de l'humidite relative ciel-clair (rhcl)
+        DO i=1,klon
+           IF (rneb(i,k) .LE. 0.0) THEN
+              zqn(i) = 0.0
+              rneb(i,k) = 0.0
+              zcond(i) = 0.0
+              rhcl(i,k)=zq(i)/zqs(i)
+           ELSE IF (rneb(i,k) .GE. 1.0) THEN
+              zqn(i) = zq(i)
+              rneb(i,k) = 1.0
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+              rhcl(i,k)=1.0
+           ELSE
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+           ENDIF
+        ENDDO
+        ! If vertical heterogeneity, change fraction by volume as well
+        if (iflag_cloudth_vert>=3) then
+          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
+          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
+        endif
+
+!        ENDIF
+
+     ELSE ! de IF (cpartiel) 
+        ! -------------------------
+        ! P2.B> Nuage "tout ou rien"
+        ! -------------------------
+        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
+        DO i = 1, klon
+           IF (zq(i).GT.zqs(i)) THEN
+              rneb(i,k) = 1.0
+           ELSE
+              rneb(i,k) = 0.0
+           ENDIF
+           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
+        ENDDO
+     ENDIF ! de IF (cpartiel) 
+
+#ifdef ISO
+        ! on calcule la compo iso de l'eau dans le nuage et en ciel
+        ! clair: zxtn/zqn et zxtcs/zqcs
+
+              ! Hypothèsesimplificatrice mais peut-être fausse: 
+              ! eau nuageuse a même composition que eau de la maille 
+
+              ! ce calcul n'eat pas sensé dépendre du type de PDF utilisé     
+         do i=1,klon  
+#ifdef ISOTRAC
+            if (option_tmin.eq.1) then
+                ! c'est mieux si zqn>=0.0
+            if (essai_convergence) then
+            else
+                zqn(i)=max(0.0,zqn(i))
+            endif    
+
+            if (rneb(i,k).lt.1.0) then
+                zqcs(i)=(zq(i)-zqn(i)*rneb(i,k))/(1.0-rneb(i,k))
+#ifdef ISOVERIF
+!                write(*,*) 'rneb,zq,zqn,zqcs=',rneb(i,k),zq(i), &
+!     &                zqn(i),zqcs(i)
+                call iso_verif_positif(zqcs(i),'ipleut 853')
+#endif                
+            else
+                zqcs(i)=0.0
+            endif
+            endif ! if (option_tmin.eq.1) then
+#endif                  
+#ifdef ISOVERIF
+              if (iso_eau.gt.0) then  
+                 call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
+     &                   'il pleut 748',errmax,errmaxrel) 
+              endif !if (iso_eau.gt.0) then
+              do ixt=1,ntraciso
+                call iso_verif_positif_choix(zxt(ixt,i),0.0,'ilp 857')
+              enddo 
+#ifdef ISOTRAC
+              call iso_verif_traceur_pbidouille(zxt(1,i), &
+     &           'il pleut 874')
+#endif              
+#endif
+              if ((bidouille_anti_divergence).and. &
+     &           (iso_eau.gt.0)) then
+                zxt(iso_eau,i)= max(zq(i),0.0)  
+!                if (zxt(iso_eau,i).gt.0.0) then
+!                    do ixt=1,niso
+!                        zxt(ixt,i)=zxt(ixt,i)*zq(i)/zxt(iso_eau,i)
+!                    enddo
+!                else !if (zxt(iso_eau,i).gt.0.0) then
+!                    do ixt=1,niso
+!                        zxt(ixt,i)=0.0
+!                    enddo    
+!                endif !if (zxt(iso_eau,i).gt.0.0) then
+              endif !if ((bidouille_anti_divergence).and. 
+
+              if (zq(i).gt.1e-15) then   
+                ! 2 fev 2010: on modifie seuil de zq  
+                ! avant, c'était 0. Vérifier que ça n'a pas de
+                ! conséquences trop graves pour les isotopes.
+                do ixt=1,ntraciso             
+                   zxtn(ixt,i)=zqn(i)*(zxt(ixt,i)/zq(i))  
+                enddo !do ixt=1,niso         
+#ifdef ISOTRAC
+               if (option_tmin.eq.1) then
+                 do ixt=1,ntraciso
+                   zxtcs(ixt,i)=(zqcs(i))*(zxt(ixt,i)/zq(i))
+                 enddo !do ixt=1,ntraciso
+               endif ! if (option_tmin.eq.1) then
+#endif                
+              else !if (zq(i).gt.0.0) then 
+                  ! problème: zq~0 et sans compo, mais zqn >0
+                  ! on fait ce qu'on peut dans ce cas pathologique
+                do ixt=1,ntraciso           
+                   zxtn(ixt,i)=0.0              
+                enddo !do ixt=1,niso 
+#ifdef ISOTRAC
+                if (option_tmin.eq.1) then
+                 do ixt=1,ntraciso
+                   zxtcs(ixt,i)=0.0
+                 enddo !do ixt=1,niso    
+                endif ! if (option_tmin.eq.1) then
+#endif                     
+                if ((bidouille_anti_divergence) &
+     &                   .and.(iso_eau.gt.0)) then
+!                   write(*,*) 'ilp tmp 848: warning: bidouille'
+                   zxtn(iso_eau,i)=zqn(i)
+#ifdef ISOTRAC
+                   zxtn(itZonIso(izone_poubelle,iso_eau),i)=zqn(i)  
+                   if (option_tmin.eq.1) then                   
+                     zxtcs(iso_eau,i)=zqcs(i)
+                   endif ! if (option_tmin.eq.1) then
+!                   write(*,*) 'zxtn(:,i)=',zxtn(:,i)
+#endif                   
+                endif !if ((bidouille_anti_divergence)
+              endif  !if (zq(i).gt.0.0) then 
+
+#ifdef ISOVERIF
+             do ixt=1,niso     
+                  call iso_verif_positif(zxtn(ixt,i),'il pleut 1062')  
+             enddo !do ixt=1,niso
+             if (iso_eau.gt.0) then
+                   call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
+     &                   'il pleut 932',errmax,errmaxrel)
+             endif !if ((iso_eau.gt.0)) then
+#ifdef ISOTRAC
+!             write(*,*) 'zxt(:,i)=',zxt(:,i)   
+!             write(*,*) 'zq,zqn=',zq(i) ,zqn(i)  
+             call iso_verif_traceur_pbidouille(zxtn(1,i), &
+     &           'il pleut 1277')
+             if (option_tmin.eq.1) then
+              call iso_verif_positif(zxtcs(ixt,i),'il pleut 895')  
+              do ixt=1,niso     
+               call iso_verif_positif(zxtcs(ixt,i),'il pleut 897')  
+              enddo !do ixt=1,niso
+              if (iso_eau.gt.0) then
+               call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
+     &              zqcs(i),'il pleut 733',errmax,errmaxrel)
+              endif !if ((iso_eau.gt.0)) then
+             endif ! if (option_tmin.eq.1) then
+#endif 
+#endif                   
+!                 write(*,*) 'il pleut 695: calcul de zxtn, i=',i
+!                 write(*,*) 'zxtn(4,i)=',zxtn(4,i)
+!                 write(*,*) 'zqn(i)=',zqn(i)
+              do ixt=1,ntraciso    
+               zxtn(ixt,i)=max(0.0,zxtn(ixt,i))          
+              enddo !do ixt=1,niso
+#ifdef ISOTRAC
+              if (option_tmin.eq.1) then
+              do ixt=1,ntraciso
+               zxtcs(ixt,i)=max(0.0,zxtcs(ixt,i))
+              enddo
+              endif ! if (option_tmin.eq.1) then
+#endif                            
+
+#ifdef ISOVERIF
+               if (iso_eau.gt.0) then
+                  call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
+     &                   'il pleut 746',errmax,errmaxrel) 
+               endif !if (iso_eau.gt.0) then
+               if (iso_HDO.gt.0) then
+                if (zqn(i).gt.ridicule) then
+                   call iso_verif_aberrant(zxtn(iso_HDO,i)/zqn(i), &
+     &                   'il pleut 977') 
+                endif  !if (zqn(i).gt.ridicule) then
+               endif  !if (iso_HDO.gt.0) then  
+#ifdef ISOTRAC               
+             call iso_verif_traceur_pbidouille  &
+     &                  (zxtn(1,i),'il pleut 1300')
+             if (option_tmin.eq.1) then
+               if (iso_eau.gt.0) then   
+                   call iso_verif_egalite_choix(zxtcs(iso_eau,i), &
+     &                   zqcs(i), &
+     &                   'il pleut 971',errmax,errmaxrel)
+               endif !if (iso_eau.gt.0) then
+               call iso_verif_traceur_pbidouille &
+     &                  (zxtcs(1,i),'il pleut 1300')
+             
+             endif ! if (option_tmin.eq.1) then
+#endif                
+#endif   
+         enddo ! do i=1,klon
+         
+!        write(*,*) 'ilp tmp 1311: zoliq,zxtoliq(1,363)=',zoliq(363),zxtoliq(1,363)
+#endif    
+     !
+     ! Mise a jour vapeur d'eau
+     ! -------------------------
+     DO i = 1, klon
+        zq(i) = zq(i) - zcond(i)
+        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
+     ENDDO
+!AJ<
+     ! ------------------------------------
+     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
+     ! -------------------------------------
+     ! Variable calcule:
+     !   zt : temperature de la maille
+     !
+     IF (.NOT. ice_thermo) THEN
+        if (iflag_fisrtilp_qsat.lt.1) then
+           DO i = 1, klon
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+           ENDDO
+        else if (iflag_fisrtilp_qsat.gt.0) then
+           DO i= 1, klon
+    if (fl_cor_ebil .GT. 0) then
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+    else
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+    end if
+           ENDDO
+        endif
+     ELSE
+         if (iflag_t_glace.ge.1) then
+            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
+         endif
+         if (iflag_fisrtilp_qsat.lt.1) then
+           DO i = 1, klon
+! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
+!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
+!                                     t_glace_max, exposant_glace)
+              if (iflag_t_glace.eq.0) then
+                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
+                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+                    zfice(i) = zfice(i)**exposant_glace_old
+              endif
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
+                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
+           ENDDO
+         else
+           DO i=1, klon
+! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
+!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
+!                                     t_glace_max, exposant_glace)
+              if (iflag_t_glace.eq.0) then
+                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
+                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+                    zfice(i) = zfice(i)**exposant_glace_old
+              endif
+        if (fl_cor_ebil .GT. 0) then
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 
+           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
+                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+        else 
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
+                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+        end if
+           ENDDO
+         endif
+!         print*,zt(i),zrfl(i),zifl(i),'temp1'
+       ENDIF
+!>AJ
+
+      ! isotopes lors de la condensation
+#ifdef ISO          
+#ifdef ISOVERIF
+!        write(*,*) 'ilp tmp 1823'
+      do i=1,klon
+        ! verif
+        if (iso_eau.gt.0) then
+          call iso_verif_egalite_choix(zxtn(iso_eau,i),zqn(i), &
+     &           'il pleut 810',errmax,errmaxrel) 
+        endif !if (iso_eau.gt.0) then 
+#ifdef ISOTRAC
+        call iso_verif_traceur(zxtn(1,i),'il pleut 1408')
+!        write(*,*) 'fisrtilp 1421: appel condiso_liq_ice_vectall'
+#endif
+       enddo !do i=1,klon 
+!       write(*,*) 'il pleut 961: avant condensation'
+       ! end verif
+#endif
+      
+      do i=1,klon  
+!        call calcul_zfice(zt(i),zfice(i)) ! inliner pour optimisation:
+        zfice(i)=1.0-(zt(i)-pxtice)/(pxtmelt-pxtice)
+        zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+       enddo  
+        
+        ! calcul de la composition du condensat
+        ! hypothèse simplificatrice (mais à garder en mémoire!): l'eau
+        ! nuageuse a la même composition que l'eau totale de la maille:
+        ! c'est à dire zxtn=(zxt/zq)*zqn
+        
+        call condiso_liq_ice_vectall(zxtn,zqn,zcond,zt,zfice, &
+     &            zxtice,zxtliq,klon)
+#ifdef ISOTRAC
+!#ifdef ISOVERIF
+!        write(*,*) 'il pleut 980: condensation pour les traceurs'
+!#endif        
+        call condiso_liq_ice_vectall_trac(zxtn,zqn,zcond, &
+     &            zt,zfice,zxtice,zxtliq,klon)
+#ifdef ISOVERIF
+!        write(*,*) 'ilp tmp 1859'
+        do i=1,klon             
+             call iso_verif_traceur(zxtn(1,i),'il pleut 1428a')
+             call iso_verif_traceur(zxtliq(1,i),'il pleut 1428b')
+             call iso_verif_traceur(zxtice(1,i),'il pleut 1428c')
+        enddo
+#endif
+#endif
+        do i=1,klon
+        do ixt=1,ntraciso                 
+          zxtcond(ixt,i)=zxtice(ixt,i)+zxtliq(ixt,i)
+          zxtcond(ixt,i)=max(zxtcond(ixt,i),0.0)
+          zxtcond(ixt,i)=min(zxtcond(ixt,i),zxt(ixt,i))
+        enddo !do ixt=1,niso 
+        enddo !do i=1,klon
+
+       do i=1,klon   
+       enddo  !do i=1,klon       
+        
+!#ifdef ISOVERIF
+#ifdef ISOVERIF
+        do i=1,klon
+          do ixt=1,ntraciso  
+            call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 638')
+          enddo
+        enddo !do i=1,klon 
+#endif       
+#ifdef ISOVERIF       
+        if (iso_eau.gt.0) then
+          do i=1,klon      
+              call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
+     &                   'il pleut 645',errmax,errmaxrel)
+          enddo ! do i=1,klon     
+        endif
+        if (iso_HDO.gt.0) then
+           do i=1,klon  
+                  call iso_verif_aberrant_choix(zxtcond(iso_hdo,i), &
+     &           zcond(i),ridicule_rain,deltalim_snow,'il pleut 1276')
+           enddo ! do i=1,klon  
+        endif !if ((iso_eau.gt.0).and.(ixt.eq.iso_eau)) then
+#ifdef ISOTRAC
+        do i=1,klon
+             if (iso_verif_traceur_nostop(zxtcond(1,i),'il pleut 1455') &
+     &                   .eq.1) then
+               write(*,*) 'zxtn(:,i)=',zxtn(:,i)
+               write(*,*) 'zxt(:,i)=',zxt(:,i)
+               write(*,*) 'zxtliq(:,i)=',zxtliq(:,i)
+               write(*,*) 'zxtice(:,i)=',zxtice(:,i)
+               stop
+             endif
+        enddo     
+#endif         
+#endif
+
+       ! retranchage à zxt 
+      do i=1,klon 
+       do ixt=1,ntraciso
+          zxt(ixt,i)=zxt(ixt,i)-zxtcond(ixt,i)
+       enddo !do ixt=1,niso   
+      enddo !do i=1,klon 
+
+#ifdef ISOVERIF
+      do i=1,klon
+         do ixt=1,niso
+          call iso_verif_positif(zxt(ixt,i),'il pleut 1088')
+         enddo !do ixt=1,niso
+#ifdef ISOTRAC  
+         ! on sépare pour mieux débugger 8 fev 2010
+         do ixt=1,ntraciso
+          call iso_verif_positif(zxt(ixt,i),'il pleut 1093')
+         enddo !do ixt=1,niso
+#endif      
+      enddo   !do i=1,klon
+#endif
+
+#ifdef ISOTRAC        
+        if (option_tmin.ge.1) then
+        do i=1,klon    
+        ! colorier la vapeur résiduelle selon température de
+        ! condensation, et le condensat en un tag spécifique
+        ! Attention: zqn,zxtn ne servent qu'au calcul de zcond,zxtcond.
+        ! C'est en fait zxt qui se fait retranché
+        ! correction le 31 mars 2010: c'est à qn qu'il faut retrancher
+        ! le condensat, car la condensation LS est un processus sous
+        ! maille -> ne tagguer que le résidus de qn.
+          if ((zcond(i).gt.0.0).and.(zq(i).gt.0.0)) then 
+
+            if (rneb(i,k).gt.0.0) then
+              zcondn(i)=zcond(i)/rneb(i,k)
+              do ixt=1,ntraciso
+                zxtcondn(ixt,i)=zxtcond(ixt,i)/rneb(i,k)
+              enddo !do ixt=1,ntraciso
+            else !if (rneb.eq.0.0)
+              zcondn(i)=0.0
+              do ixt=1,ntraciso
+                zxtcondn(ixt,i)=0.0
+              enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+              call iso_verif_egalite(zcond(i),0.0,'ilp 1225')
+#endif              
+            endif !if (rneb.eq.0.0)
+            
+#ifdef ISOVERIF
+          call iso_verif_egalite_choix(zxtcondn(iso_eau,i),zcondn(i), &
+     &                   'il pleut 1244',errmax*10,errmaxrel*10)
+          call iso_verif_traceur(zxtcondn(1,i),'cv3_routines 2302')
+          zq_verif=rneb(i,k)*zqn(i)+(1.0-rneb(i,k))*zqcs(i)-zcond(i)
+          call iso_verif_egalite_choix(zq(i),zq_verif, &
+     &                   'il pleut 1235',errmax,errmaxrel)
+          do ixt=1,niso
+            zq_verif=rneb(i,k)*zxtn(ixt,i) &
+     &           +(1.0-rneb(i,k))*zxtcs(ixt,i)-zxtcond(ixt,i)
+            call iso_verif_egalite_choix(zxt(ixt,i),zq_verif, &
+     &                   'il pleut 1235',errmax,errmaxrel)
+          enddo
+#endif            
+          ! bidouille
+          if (bidouille_anti_divergence) then
+                zxtcondn(iso_eau,i)=zcondn(i)
+#ifdef ISOVERIF                
+                call iso_verif_traceur_pbidouille(zxtcondn(1,i), &
+     &          'il pleut 1261')
+#endif                
+          endif
+
+            if (option_traceurs.eq.17) then      
+            call iso_recolorise_condensation(zqn(i),zcondn(i), &
+     &           zxtn(1,i),zxtcondn(1,i),zt(i),1.0,zxtres, &
+     &           seuil_tag_tmin_ls)
+            else !if (option_traceurs.eq.17) then
+                ! cas du tag 18 ou 19
+                call iso_recolorise_condensation(zqn(i),zcondn(i), &
+     &           zxtn(1,i),zxtcondn(1,i),zqs(i),1.0,zxtres, &
+     &           seuil_tag_tmin_ls)
+            endif !if (option_traceurs.eq.17) then
+            do ixt=1+niso,ntraciso
+               zxt(ixt,i)=rneb(i,k)*zxtres(ixt) &
+     &           +(1.0-rneb(i,k))*zxtcs(ixt,i)
+            enddo         
+#ifdef ISOVERIF
+         call iso_verif_traceur(zxt(1,i),'ilp 1025')
+         do izone=izone_oce,izone_ddft
+          if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
+           ! avec condensation, on ne peut que perdre des traceurs à
+           ! part le tag résuel et le condensat
+           if (iso_verif_positif_choix_nostop( &
+     &           zxt_ancien(itZonIso(izone,iso_eau),i) &
+     &          -zxt(itZonIso(izone,iso_eau),i),1e-8,'ilp 1270') &
+     &          .eq.1) then
+            write(*,*) 'i,izone,rneb=',i,izone,rneb(i,k)
+            write(*,*) 'zxt(:,i)=',zxt(:,i)
+            write(*,*) 'zxt_ancien(:,i)=',zxt_ancien(:,i)
+            write(*,*) 'zxtcs(:,i)=',zxtcs(:,i)
+            write(*,*) 'zxtres(:)=',zxtres(:)
+            write(*,*) 'zxtcond(:,i)=',zxtcond(:,i)
+            write(*,*) 'zxtn(:,i)=',zxtn(:,i)
+            write(*,*) 'zxtcondn(:,i)=',zxtcondn(:,i)
+            stop
+           endif !if (iso_verif_positif_nostop(
+          endif !if ((izone.eq.izone_oce).or.(izone.eq.izone_ddft)) then
+         enddo !do izone=izone_oce,izone_revap
+#endif                     
+          endif !if (cond.gt.0.0) then
+        enddo !do i=1,klon
+
+#ifdef ISOVERIF
+        ! vérifier qu'on recolorise bien le condensat
+        call iso_verif_positif(float(option_cond-1),'ilp 1310')
+#endif        
+        endif ! if (option_tmin.ge.1) then
+
+        if (option_cond.eq.1) then
+        ! colorier le condensat en un tag spécifique
+        do i=1,klon        
+          do ixt=1+niso,ntraciso
+            if (index_zone(ixt).eq.izone_cond) then
+              zxtcond(ixt,i)=zxtcond(index_iso(ixt),i)
+            else
+              zxtcond(ixt,i)=0.0
+            endif
+          enddo !do ixt=1,ntraciso
+        enddo !do i=1,klon
+#ifdef ISOVERIF
+        do i=1,klon    
+          call iso_verif_egalite_choix(zxtcond(iso_eau,i),zcond(i), &
+     &                   'il pleut 1139',errmax,errmaxrel)
+          call iso_verif_traceur(zxtcond(1,i),'cv3_routines 2302')
+        enddo !do i=1,klon        
+#endif         
+      endif ! if (option_cond.ge.1) then
+#endif        
+     
+      
+#ifdef ISOVERIF
+!        write(*,*) 'ilp tmp 2066'
+      do i=1,klon
+         do ixt=1,niso
+          call iso_verif_positif(zxt(ixt,i),'il pleut 778')
+         enddo !do ixt=1,niso
+#ifdef ISOTRAC  
+         ! on sépare pour mieux débugger 8 fev 2010
+         do ixt=1,ntraciso
+          call iso_verif_positif(zxt(ixt,i),'il pleut 1135')
+          call iso_verif_traceur(zxt(1,i),'il pleut 1485')
+         enddo !do ixt=1,niso
+#endif
+         if (iso_eau.gt.0) then
+            call iso_verif_egalite_choix(zxt(iso_eau,i),zq(i), &
+     &                   'il pleut 644',errmax,errmaxrel)
+         endif !if (iso_eau.gt.0) then
+         if (iso_HDO.gt.0) then  
+            if (zq(i).gt.ridicule) then
+             call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
+     &                   'il pleut 648')                 
+              if (iso_O18.gt.0) then                
+                if (iso_verif_O18_aberrant_nostop(zxt(iso_HDO,i)/zq(i), &
+     &           zxt(iso_O18,i)/zq(i),'il pleut 2059').eq.1) then
+                   write(*,*) 'i,k,zq=',i,k,zq(i)
+                   stop
+                endif !if (iso_verif_O18_aberrant_nostop
+              endif !if (iso_O18.gt.0) then
+            endif !if (zq(i).gt.ridicule) then
+         endif  ! if (iso_HDO.gt.0) then 
+        enddo  ! i=1,klon
+#ifdef ISOTRAC
+        do i=1,klon
+             call iso_verif_traceur(zxt(1,i),'il pleut 1485')
+        enddo     
+#endif         
+#endif
+#ifdef ISOVERIF   
+       do i=1,klon   
+           do ixt=1,ntraciso
+             call iso_verif_noNAN(zxtcond(ixt,i),'il pleut 1560')
+             call iso_verif_noNAN(zxt(ixt,i),'il pleut 1561')
+           enddo   ! ixt=1,niso
+        enddo  !do i=1,klon         
+#endif   
+#endif    
+     ! ----------------------------------------------------------------
+     ! P3> Formation des precipitations
+     ! ----------------------------------------------------------------
+     !
+     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
+     !
+
+     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           zoliq(i) = zcond(i)
+           zrho(i) = pplay(i,k) / zt(i) / RD
+           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+        ENDIF
+     ENDDO
+!AJ<
+     IF (.NOT. ice_thermo) THEN
+       IF (iflag_t_glace.EQ.0) THEN
+         DO i = 1, klon
+            IF (rneb(i,k).GT.0.0) THEN
+               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
+               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+               zfice(i) = zfice(i)**exposant_glace_old
+!              zfice(i) = zfice(i)**nexpo
+         !!      zfice(i)=0.
+            ENDIF
+         ENDDO
+       ELSE ! of IF (iflag_t_glace.EQ.0)
+         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
+!         DO i = 1, klon
+!            IF (rneb(i,k).GT.0.0) THEN
+! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
+!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
+!                                     t_glace_max, exposant_glace)
+!            ENDIF
+!         ENDDO
+       ENDIF
+     ENDIF
+
+     ! Calcul de radliq (eau condensee pour le rayonnement)
+     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
+     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
+     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
+     ! transmise au rayonnement;
+     ! ----------------------------------------------------------------
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           zneb(i) = MAX(rneb(i,k), seuil_neb)
+     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i))  
+!           print*,zt(i),'fractionglace'
+!>AJ
+           radliq(i,k) = zoliq(i)/REAL(ninter+1)
+        ENDIF
+     ENDDO
+#ifdef ISO
+      ! initilsation des quantités totales d'eau liquide et solide qui
+      ! tombent
+      do i=1,klon
+         if (rneb(i,k).gt.0.0) then
+            do ixt=1,ntraciso
+               zxtoliq(ixt,i)=zxtcond(ixt,i)
+            enddo
+         endif ! if (rneb(i,k).gt.0.0)
+       enddo
+          
+#ifdef ISOVERIF
+!       write(*,*) 'ilp tmp 1626: k,zoliq,zxtoliq(1,363)=',k,zoliq(363),zxtoliq(1,363)
+       do i=1,klon   
+         do ixt=1,ntraciso
+           if (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
+     &           'il pleut 1628').eq.1) then
+             write(*,*) 'i,k,ixt=',i,k,ixt
+             write(*,*) 'rneb(i,k)=',rneb(i,k)
+             write(*,*) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
+             stop
+           endif
+         enddo   ! ixt=1,niso
+        enddo  !do i=1,klon  
+#endif                
+#ifdef ISOVERIF
+       !i=519
+       !write(*,*) 'k,i,zoliq(i)=',k,i,zoliq(i)
+       do i=1,klon   
+           ! cam verif
+           if (iso_eau.gt.0) then
+            call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
+     &            'il pleut 678',errmax,errmaxrel)
+           endif !if (iso_eau.gt.0) then
+           if (zoliq(i).gt.ridicule) then
+            if (iso_HDO.gt.0) then  
+                ! Camille 9 mars 2023: on est moins stricte avec le condensat
+                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
+     &           ridicule_rain,deltalim_snow, 'il pleut 895a')
+              if (iso_O18.gt.0) then
+                if (iso_verif_O18_aberrant_nostop(zxtoliq(iso_HDO,i)/zoliq(i), &
+     &           zxtoliq(iso_O18,i)/zoliq(i),'il pleut 895b').eq.1) then
+                   write(*,*) 'i,k,zoliq,zfice=',i,k,zoliq(i),zfice(i)
+                   !stop
+                endif
+              endif ! if (iso_HDO.gt.0) then
+            endif  !if (iso_HDO.gt.0) then
+           endif !if (zoliq(i).gt.ridicule) then
+          ! end cam verif
+         enddo ! do i=1,klon
+#endif
+#endif
+
+     !
+     DO n = 1, ninter
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0) THEN
+              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+              ! Initialization of zpluie and zice:
+              zpluie=0
+              zice=0
+              IF (zneb(i).EQ.seuil_neb) THEN
+                 ztot = 0.0
+              ELSE
+                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
+                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
+                 if (ptconv(i,k)) then
+                    zcl   =cld_lc_con
+                    zct   =1./cld_tau_con
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvc(zrhol(i)) * zfice(i)
+                 else
+                    zcl   =cld_lc_lsc
+                    zct   =1./cld_tau_lsc
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvs(zrhol(i)) * zfice(i)
+                 endif
+
+                 ! si l'heterogeneite verticale est active, on utilise 
+                 ! la fraction volumique "vraie" plutot que la fraction
+                 ! surfacique modifiee, qui est plus grande et reduit
+                 ! sinon l'eau in-cloud de facon artificielle
+                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
+                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
+                 else
+                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
+                 endif
+!AJ<
+                 IF (.NOT. ice_thermo) THEN
+                   ztot    = zchau + zfroi
+                 ELSE
+                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
+                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
+                   ztot    = zpluie    + zice
+                 ENDIF
+!>AJ
+                 ztot    = MAX(ztot   ,0.0)
+              ENDIF
+              ztot    = MIN(ztot,zoliq(i))
+!AJ<
+     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
+     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
+!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
+!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
+!      si iflag_bergeron <> 2
+!      A SUPPRIMER A TERME
+              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
+              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
+              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
+!>AJ
+              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+           ENDIF
+        ENDDO  ! i = 1,klon
+     ENDDO     ! n = 1,ninter
+
+     ! ----------------------------------------------------------------
+     !
+#ifdef ISO
+      do i=1,klon
+       IF (rneb(i,k).GT.0.0) THEN
+        if (zcond(i).gt.ridicule**2) then
+            ! le 21 dec: on change .gt.0 en .gt.ridicule**2 pour éviter valeur
+            ! ridiculement petites
+           do ixt=1,ntraciso
+             zxtoliq(ixt,i)=zoliq(i)*(zxtcond(ixt,i)/zcond(i))
+           enddo !do ixt=1,niso           
+        else !if (zcond(i).gt.0.0) then            
+           if (zoliq(i).lt.ridicule) then
+               do ixt=1,ntraciso
+                zxtoliq(ixt,i)=0.0
+               enddo
+           else !if (zoliq(i).lt.ridicule) then
+               ! modif 28 oct 2008
+              do ixt=1,niso
+                zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
+              enddo !do ixt=1,niso
+#ifdef ISOTRAC
+              do ixt=niso+1,ntraciso
+                if (index_zone(ixt).eq.izone_poubelle) then
+                    zxtoliq(ixt,i)=zoliq(i)*Rdefault(ixt)
+                else
+                    zxtoliq(ixt,i)=0.0
+                endif
+              enddo  
+#endif    
+#ifdef ISOVERIF                   
+              write(*,*) 'il pleut 1412: zcond(i)=',zcond(i)
+              write(*,*) 'zoliq(i)=',zoliq(i)
+              stop
+#endif                   
+            endif   !if (zoliq(i).lt.ridicule) then         
+        endif !if (zcond(i).gt.0.0) then
+       endif !IF (rneb(i,k).GT.0.0) THEN
+      enddo ! do i=1,klon 
+           
+         ! cam verif
+#ifdef ISOVERIF
+      do i=1,klon  
+        IF (rneb(i,k).GT.0.0) THEN 
+           do ixt=1,niso
+             call iso_verif_noNAN(zxtoliq(ixt,i),'il pleut 718')
+           ! end cam verif
+           enddo ! do ixt=1,niso
+          ! cam verif
+        endif !IF (rneb(i,k).GT.0.0) THEN 
+      enddo !do i=1,klon  
+#endif
+#ifdef ISOVERIF
+      do i=1,klon  
+        IF (rneb(i,k).GT.0.0) THEN
+          if (iso_eau.gt.0) then
+            call iso_verif_egalite_choix(zxtoliq(iso_eau,i),zoliq(i), &
+     &           'il pleut 749',errmax,errmaxrel)
+          endif !if (iso_eau.gt.0) then
+          if (iso_HDO.gt.0) then
+                call iso_verif_aberrant_choix(zxtoliq(iso_HDO,i),zoliq(i), &
+     &           ridicule_rain, deltalim_snow, 'il pleut 963')
+          endif !if (iso_HDO.gt.0) then
+        ! end cam verif
+#ifdef ISOTRAC
+!             write(*,*) 'ilp 1201 tmp: i=',i
+!             write(*,*) 'zoliq(i)=',zoliq(i)
+!             write(*,*) 'zxtoliq(i)=',zxtoliq(iso_eau:ntraciso:3,i)
+!             write(*,*) 'zcond(i)=',zcond(i)
+!             write(*,*) 'zxtcond(i)=',zxtcond(iso_eau:ntraciso:3,i)
+             call iso_verif_traceur(zxtoliq(1,i),'il pleut 1634')
+#endif        
+         endif !IF (rneb(i,k).GT.0.0) THEN
+       enddo ! do i=1,klon
+#endif        
+#endif      
+     IF (.NOT. ice_thermo) THEN
+       DO i = 1, klon
+         IF (rneb(i,k).GT.0.0) THEN
+           d_ql(i,k) = zoliq(i)
+           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
+                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+
+#ifdef ISO   
+         do ixt=1,ntraciso
+            d_xtl(ixt,i,k)=zxtoliq(ixt,i)
+         enddo ! do ixt=1,niso   
+            ! Rq: 28 oct 2008: si zoliq a la même compo que zcond, alors
+            ! on évite pb numériques si on prend directement zrfl avec
+            ! la compo de zcond?
+!          do ixt=1,niso  
+!            zxtrfl(ixt,i)=zxtrfl(ixt,i) &
+!     &           +max(zxtcond(ixt,i)-zxtoliq(ixt,i),0.0) &
+!     &           * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+!          enddo  
+            if (zcond(i).gt.ridicule**2) then
+                ! modif le 19 dec 2012 pour éviter zcond pathologiquement petits
+                do ixt=1,ntraciso
+                  zxtrfl(ixt,i)=zxtrfl(ixt,i) &
+     &                  +(zxtcond(ixt,i)/zcond(i)) &
+     &                  * MAX(zcond(i)-zoliq(i),0.0) &
+     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+                enddo !do ixt=1,niso
+            else !if (zcond(i).gt.0.0) then
+                if (zcond(i)-zoliq(i).gt.0.0) then
+                    do ixt=1,niso
+                        zxtrfl(ixt,i)=zxtrfl(ixt,i) &
+     &                  +Rdefault(ixt) &
+     &                  * MAX(zcond(i)-zoliq(i),0.0) &
+     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+                    enddo !do ixt=1,niso  
+#ifdef ISOVERIF                   
+                    write(*,*) 'il pleut 1023: zcond(i)=',zcond(i)
+                    write(*,*) 'zcond(i),zoliq(i)=',zcond(i),zoliq(i)
+                    stop
+#endif        
+#ifdef ISOTRAC
+                   do ixt=niso+1,ntraciso
+                      if (index_zone(ixt).eq.izone_poubelle) then
+                        zxtrfl(ixt,i)=zxtrfl(ixt,i) &
+     &                  +Rdefault(ixt) &
+     &                  * MAX(zcond(i)-zoliq(i),0.0) &
+     &                  * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+                      endif
+                    enddo !do ixt=1,niso
+#endif                 
+                else !if (zcond(i)-zoliq(i).gt.0.0) then
+                    ! zrfl non modifié
+                endif !if (zcond(i)-zoliq(i).gt.0.0) then
+            endif !if (zcond(i).gt.0.0) then
+
+         ! cam verif         
+!#ifdef ISOVERIF
+#ifdef ISOVERIF
+            do ixt=1,ntraciso
+             if ((iso_verif_noNAN_nostop(zxtrfl(ixt,i), &
+     &                   'il pleut 772').eq.1).or. &
+     &                   (iso_verif_noNAN_nostop(zxtoliq(ixt,i), &
+     &                   'il pleut 772b').eq.1)) then
+               write(*,*) 'zcond(i)=',zcond(i)
+               write(*,*) 'zxtcond(ixt,i)=',zxtcond(ixt,i)
+               write(*,*) 'zoliq(i)=',zoliq(i)
+               write(*,*) 'Rdefault(ixt)=',Rdefault(ixt)
+               write(*,*) 'paprs(i,k),paprs(i,k+1)=',paprs(i,k),paprs(i,k+1)
+               stop
+             endif
+            enddo
+#endif
+#ifdef ISOVERIF
+            if (iso_eau.gt.0) then
+              call iso_verif_egalite_choix(d_xtl(iso_eau,i,k), &
+     &           d_ql(i,k),'il pleut 771',errmax,errmaxrel)
+              if (iso_verif_egalite_choix_nostop( &
+     &           zxtrfl(iso_eau,i),zrfl(i), &
+     &           'il pleut 771b',errmax,errmaxrel).eq.1)then
+                 write(*,*) 'trace_cas(i)=',trace_cas(i)
+                 write(*,*) 'zxtcond(iso_eau,i)',zxtcond(iso_eau,i)
+                 write(*,*) 'zcond(i)',zcond(i)
+                 write(*,*) 'zxtoliq(iso_eau,i)',zxtoliq(iso_eau,i)
+                 write(*,*) 'zoliq(i)',zoliq(i)
+                 stop
+              endif
+            endif !if (iso_eau.gt.0) then
+            if (iso_HDO.gt.0) then
+                IF ((t(i,1)) .LT. RTT) THEN
+                call iso_verif_aberrant_choix(zxtrfl(iso_hdo,i), &
+     &           zrfl(i),ridicule_rain,deltalim_snow,'il pleut 1458')
+                endif
+            endif !if (iso_HDO.gt.0) then
+#ifdef ISOTRAC
+             call iso_verif_traceur(zxtrfl(1,i),'il pleut 1729')
+#endif            
+#endif
+         ! bidouille
+         if ((bidouille_anti_divergence).and.(iso_eau.gt.0)) then 
+               zxtrfl(iso_eau,i)= zrfl(i)   
+!             if (zxtrfl(iso_eau,i).gt.0.0) then
+!             do ixt=1,niso
+!               zxtrfl(ixt,i)=zxtrfl(ixt,i)*zrfl(i)/zxtrfl(iso_eau,i)
+!             enddo
+!           endif !if (zxtrfl(iso_eau,i).gt.0.0) then
+         endif ! if (bidouille_anti_divergence) then
+         ! end bidouille
+         ! end cam verif
+#endif
+         ENDIF !IF (rneb(i,k).GT.0.0) THEN
+       ENDDO
+     ELSE
+
+#ifdef ISO
+        CALL abort_physic('ilp 2347', 'isotopes pas encore prevus ici', 1)
+#endif
+!
+!CR&JYG<
+! On prend en compte l'effet Bergeron dans les flux de precipitation :
+! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
+! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
+! et les precipitations est grossierement pris en compte en linearisant les equations 
+! et en approximant le processus de precipitation liquide par un processus a seuil.
+! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
+! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
+! Le condensat precipitant solide est augmente.
+! La vapeur d'eau est augmentee.
+!
+       IF ((iflag_bergeron .EQ. 2)) THEN
+         DO i = 1, klon
+           IF (rneb(i,k) .GT. 0.0) THEN
+             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
+             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
+           if (fl_cor_ebil .GT. 0) then
+             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
+!            Calcul de DT si toute les precips liquides congelent
+             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
+!            On ne veut pas que T devienne superieur a la temp. de congelation.
+!            donc que Delta > RTT-zt(i
+             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
+             zt(i)      = zt(i)      + DeltaT
+!            Eau vaporisee du fait de l'augmentation de T
+             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
+!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
+             zq(i) = zq(i) +  Deltaq
+!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
+             zcond(i)   = max( zcond(i)- Deltaq, 0. )
+!            precip liquide qui congele ou qui s'evapore
+             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
+             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
+!            bilan eau glacee
+             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
+           else ! if (fl_cor_ebil .GT. 0)
+!            ancien calcul
+             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
+             coef1 = RLMLT*zdqs(i)/RLVTT
+             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
+             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
+             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
+             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
+             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
+             zt(i)      = zt(i)      + DeltaT
+           end if ! if (fl_cor_ebil .GT. 0)
+           ENDIF  ! rneb(i,k) .GT. 0.0
+         ENDDO
+         DO i = 1, klon
+           IF (rneb(i,k).GT.0.0) THEN
+             d_ql(i,k) = (1-zfice(i))*zoliq(i)
+             d_qi(i,k) = zfice(i)*zoliq(i)
+             zrfl(i) = zrfl(i)+ zqprecl(i) &
+                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+             zifl(i) = zifl(i)+ zqpreci(i) &
+                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+           ENDIF                      
+         ENDDO
+!!
+       ELSE  ! iflag_bergeron
+!>CR&JYG
+!!
+       DO i = 1, klon
+         IF (rneb(i,k).GT.0.0) THEN
+!CR on prend en compte la phase glace
+!JLD inutile car on ne passe jamais ici si .not.ice_thermo
+!           if (.not.ice_thermo) then
+!           d_ql(i,k) = zoliq(i)
+!           d_qi(i,k) = 0.
+!           else
+           d_ql(i,k) = (1-zfice(i))*zoliq(i)
+           d_qi(i,k) = zfice(i)*zoliq(i)
+!           endif
+!AJ<
+           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
+               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
+                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+     !      zrfl(i) = zrfl(i)+  zpluie                         &
+     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+     !      zifl(i) = zifl(i)+  zice                    &
+     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                    
+
+!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
+           IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
+              zsolid = zrfl(i)
+              zifl(i) = zifl(i)+zrfl(i)
+              zrfl(i) = 0.
+           if (fl_cor_ebil .GT. 0) then
+              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+           else
+              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
+           end if
+           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
+!RC   
+
+         ENDIF ! rneb(i,k).GT.0.0
+       ENDDO
+
+       ENDIF  ! iflag_bergeron .EQ. 2
+     ENDIF  ! .NOT. ice_thermo
+
+!CR: la fonte est faite au debut
+!      IF (ice_thermo) THEN
+!       DO i = 1, klon
+!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
+!           zmelt = MIN(MAX(zmelt,0.),1.)
+!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
+!           zifl(i)=zifl(i)*(1.-zmelt)
+!           print*,zt(i),'octavio1'
+!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
+!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
+!       ENDDO
+!     ENDIF
+
+       
+     IF (.NOT. ice_thermo) THEN
+       DO i = 1, klon
+         IF (zt(i).LT.RTT) THEN
+           psfl(i,k)=zrfl(i)
+#ifdef ISO
+         do ixt=1,ntraciso
+           pxtsnowfl(ixt,i,k)=zxtrfl(ixt,i)
+         enddo
+#endif
+         ELSE
+           prfl(i,k)=zrfl(i)
+#ifdef ISO
+         do ixt=1,ntraciso
+           pxtrainfl(ixt,i,k)=zxtrfl(ixt,i)
+         enddo         
+#endif
+         ENDIF
+       ENDDO
+     ELSE ! if ice_thermo
+     ! JAM*************************************************
+     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
+     ! *****************************************************
+
+#ifdef ISO
+     CALL abort_physic('ilp 2465', 'isos pas prevus ici', 1)   
+#endif
+       DO i = 1, klon
+     !   IF (zt(i).LT.RTT) THEN
+           psfl(i,k)=zifl(i) 
+     !   ELSE
+           prfl(i,k)=zrfl(i)
+     !   ENDIF
+!>AJ
+       ENDDO
+     ENDIF
+
+#ifdef ISO
+#ifdef ISOVERIF
+      DO i = 1, klon
+        if (iso_eau.gt.0) then
+          call iso_verif_egalite_choix(pxtrainfl(iso_eau,i,k),prfl(i,k), & 
+     &           'il pleut 804',errmax,errmaxrel)
+          call iso_verif_egalite_choix(pxtsnowfl(iso_eau,i,k),psfl(i,k), &  
+     &           'il pleut 805',errmax,errmaxrel) 
+          call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &  
+     &           'il pleut 1612',errmax,errmaxrel) 
+          if (iso_verif_positif_choix_nostop(pxtrainfl(iso_eau,i,k), &
+     &           ridicule_rain,'ilp 1547').eq.1) then
+            write(*,*) 'prfl(i,k)=',prfl(i,k)
+            stop
+          endif  
+        endif !if (iso_eau.gt.0) then
+        if (iso_HDO.gt.0) then
+          if (prfl(i,k).gt.ridicule_rain) then
+            call iso_verif_aberrant(pxtrainfl(iso_HDO,i,k)/prfl(i,k), &
+     &            'il pleut 1020')
+          endif !if (prfl(i,k).gt.ridicule_rain) then           
+          call iso_verif_aberrant_choix(pxtsnowfl(iso_HDO,i,k),psfl(i,k), &
+     &           ridicule_rain, deltalim_snow, 'il pleut 1024')
+          if (zq(i).gt.ridicule) then
+              call iso_verif_aberrant_encadre(zxt(iso_HDO,i)/zq(i), &
+     &           'il pleut 1204')
+          endif !if (zq(i).gt.ridicule) then          
+         endif !if (iso_HDO.gt.0) then
+#ifdef ISOTRAC
+          call iso_verif_traceur(zxt(1,i),'il pleut 1792')
+          call iso_verif_traceur(pxtsnowfl(1,i,k),'il pleut 1793')
+          call iso_verif_traceur(pxtrainfl(1,i,k),'il pleut 1794')
+#endif         
+       ENDDO !DO i = 1, klon  
+#endif
+#endif 
+     ! ----------------------------------------------------------------
+     ! Fin de formation des precipitations
+     ! ----------------------------------------------------------------
+     !
+     ! Calculer les tendances de q et de t:
+     !
+     DO i = 1, klon
+        d_q(i,k) = zq(i) - q(i,k)
+#ifdef ISO
+         do ixt=1,ntraciso
+             d_xt(ixt,i,k)=zxt(ixt,i)-xt(ixt,i,k)   
+         enddo
+         ! cam verif
+#ifdef ISOVERIF
+           if (iso_eau.gt.0) then
+             if (iso_verif_egalite_choix_nostop(d_xt(iso_eau,i,k), &
+     &            d_q(i,k),'il pleut 838',errmax,errmaxrel).eq.1) then
+                write(*,*) 'i,k=',i,k
+                write(*,*) 'zxt,zq=',zxt(iso_eau,i),zq(i)
+                write(*,*) 'xt,q=',xt(iso_eau,i,k),q(i,k)
+                stop
+             endif
+           endif !if (iso_eau.gt.0) then
+           if (1.eq.0) then
+           if (iso_HDO.gt.0) then
+            if (abs(d_q(i,k)).gt.max(errmaxrel*q(i,k),errmax)) then
+              if (iso_verif_aberrant_nostop( &
+     &           d_xt(iso_HDO,i,k)/d_q(i,k),'il pleut 1229').eq.1) then
+!                write(*,*) 'd_q(i,k)=',d_q(i,k)
+!                write(*,*) 'q(i,k)=',q(i,k)
+!                write(*,*) 'zq(i)=',zq(i)
+!                write(*,*) 'deltaD_prec=',deltaD(xt(iso_HDO,i,k)/q(i,k))
+!                write(*,*) 'deltaD=',deltaD(zxt(iso_HDO,i)/zq(i))
+!                stop       
+              endif !   if (iso_verif_aberrant_nostop(
+            endif !if (zq(i).gt.ridicule) then          
+           endif !if (iso_HDO.gt.0) then
+          endif !if (1.eq.0) then
+#endif
+         ! end cam verif
+#endif
+        d_t(i,k) = zt(i) - t(i,k)
+     ENDDO
+     !
+     !AA--------------- Calcul du lessivage stratiforme  -------------
+
+     DO i = 1,klon
+        !
+        if(zcond(i).gt.zoliq(i)+1.e-10) then
+         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
+        else
+         beta(i,k) = 0.
+        endif
+        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
+             * (paprs(i,k)-paprs(i,k+1))/RG
+        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+           !AA lessivage nucleation LMD5 dans la couche elle-meme
+          IF (iflag_t_glace.EQ.0) THEN
+           if (t(i,k) .GE. t_glace_min_old) THEN
+              zalpha_tr = a_tr_sca(3)
+           else
+              zalpha_tr = a_tr_sca(4)
+           endif
+          ELSE ! of IF (iflag_t_glace.EQ.0)
+           if (t(i,k) .GE. t_glace_min) THEN
+              zalpha_tr = a_tr_sca(3)
+           else
+              zalpha_tr = a_tr_sca(4)
+           endif
+          ENDIF
+           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
+           !
+           ! nucleation avec un facteur -1 au lieu de -0.5
+           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
+           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+        ENDIF
+        !
+     ENDDO      ! boucle sur i
+     !
+     !AA Lessivage par impaction dans les couches en-dessous
+     DO kk = k-1, 1, -1
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+             IF (iflag_t_glace.EQ.0) THEN
+              if (t(i,kk) .GE. t_glace_min_old) THEN
+                 zalpha_tr = a_tr_sca(1)
+              else
+                 zalpha_tr = a_tr_sca(2)
+              endif
+             ELSE ! of IF (iflag_t_glace.EQ.0)
+              if (t(i,kk) .GE. t_glace_min) THEN
+                 zalpha_tr = a_tr_sca(1)
+              else
+                 zalpha_tr = a_tr_sca(2)
+              endif
+             ENDIF
+              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
+              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
+           ENDIF
+        ENDDO
+     ENDDO !DO kk = k-1, 1, -1
+
+#ifdef ISO      
+      ! verif en fin de boucle
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        do i=1,klon       
+        call iso_verif_egalite_choix(zxtrfl(iso_eau,i) &
+     &      ,zrfl(i), 'il pleut 1089',errmax,errmaxrel)
+        call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
+     &           'il pleut 1205',errmax,errmaxrel) 
+#ifdef ISOTRAC
+        call iso_verif_traceur(zxtrfl(1,i),'il pleut 1894')
+        call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1895')
+#endif        
+        enddo !do i=1,klon
+      endif !  if (iso_eau.gt.0) then        
+#endif
+#endif      
+      ! end verif en fin de boucle
+
+     !
+     !AA===============================================================
+     !                     FIN DE LA BOUCLE VERTICALE  
+  end DO
+  !
+  !AA==================================================================
+  !
+  ! Pluie ou neige au sol selon la temperature de la 1ere couche
+  !
+!CR: si la thermo de la glace est active, on calcule zifl directement
+  IF (.NOT.ice_thermo) THEN
+  DO i = 1, klon
+     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
+!AJ<
+!        snow(i) = zrfl(i)
+        snow(i) = zrfl(i)+zifl(i)
+!>AJ
+        zlh_solid(i) = RLSTT-RLVTT
+#ifdef ISO
+        do ixt=1,ntraciso
+          xtsnow(ixt,i) = zxtrfl(ixt,i)
+        enddo !do ixt=1,ntraciso
+#endif
+     ELSE
+        rain(i) = zrfl(i)
+        zlh_solid(i) = 0.
+#ifdef ISO
+        do ixt=1,ntraciso
+          xtrain(ixt,i) = zxtrfl(ixt,i)
+        enddo !do ixt=1,ntraciso
+#endif
+     ENDIF
+#ifdef ISO
+      ! cam verif 
+#ifdef ISOVERIF
+       do ixt=1,ntraciso
+        call iso_verif_noNAN(xtrain(ixt,i),'il pleut 927')
+        call iso_verif_noNAN(xtsnow(ixt,i),'il pleut 927')
+       enddo ! do ixt=1,niso
+#endif         
+#ifdef ISOVERIF
+       if (iso_eau.gt.0) then
+        call iso_verif_egalite_choix(xtrain(iso_eau,i),rain(i), &
+     &           'il pleut 926a',errmax,errmaxrel)  
+        if (snow(i)*dtime.gt.ridicule) then
+           call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
+     &            'il pleut 938a',errmax,errmaxrel)
+        endif !if (snow(i)*dtime.gt.ridicule) then
+       endif !if (iso_eau.gt.0) then
+       if (iso_HDO.gt.0) then
+           call iso_verif_aberrant_choix(xtrain(iso_hdo,i),rain(i), &
+     &           ridicule_rain,deltalim,'il pleut 926b')  
+           call iso_verif_aberrant_choix(xtsnow(iso_hdo,i),snow(i), &
+     &            ridicule_rain,deltalim_snow,'il pleut 938b')
+       endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+       call iso_verif_traceur(xtrain(1,i),'il pleut 1952')
+       call iso_verif_traceur(xtsnow(1,i),'il pleut 1953')
+#endif       
+#endif
+       ! end cam verif   
+#endif
+  ENDDO
+
+  ELSE
+#ifdef ISO
+        CALL abort_physic('ilp 2708', 'isos pas prevus ici', 1)
+#endif
+     DO i = 1, klon
+        snow(i) = zifl(i)
+        rain(i) = zrfl(i)
+     ENDDO
+   
+   ENDIF
+  !
+  ! For energy conservation : when snow is present, the solification
+  ! latent heat is considered.
+!CR: si thermo de la glace, neige deja prise en compte
+  IF (.not.ice_thermo) THEN
+  DO k = 1, klev
+     DO i = 1, klon
+        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
+        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
+        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
+     END DO
+  END DO
+  ENDIF
+  !
+
+  if (ncoreczq>0) then
+     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
+  endif
+
+
+#ifdef ISO
+      ! verif en sortie de il pleut 
+#ifdef ISOVERIF
+       do k=1,klev
+        do i=1,klon
+         if (iso_eau.gt.0) then
+           call iso_verif_egalite_choix(d_xtl(iso_eau,i,k),d_ql(i,k), &
+     &            'il pleut 1289',errmax,errmaxrel)
+         endif !if (iso_eau.gt.0) then
+#ifdef ISOTRAC
+         call iso_verif_traceur(d_xtl(1,i,k),'il pleut 1982')
+#endif
+        enddo !do i=1,klon
+       enddo !do k=1,klev
+#endif
+      ! end verif
+#endif
+
+END SUBROUTINE fisrtilp
+END MODULE lmdz_lscp_old
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_tools.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_tools.F90	(revision 4664)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_tools.F90	(revision 4664)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_lscp_tools.F90
Index: LMDZ6/trunk/libf/phylmdiso/lscp_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lscp_ini_mod.F90	(revision 4653)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/lscp_ini_mod.F90
Index: LMDZ6/trunk/libf/phylmdiso/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lscp_mod.F90	(revision 4653)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/lscp_mod.F90
Index: LMDZ6/trunk/libf/phylmdiso/lscp_tools_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lscp_tools_mod.F90	(revision 4653)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/lscp_tools_mod.F90
Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4653)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4664)
@@ -77,5 +77,6 @@
     USE write_field_phy
     USE wxios, ONLY: g_ctx, wxios_set_context, using_xios
-    USE lscp_mod, ONLY : lscp
+    USE lmdz_lscp, ONLY : lscp
+    USE lmdz_lscp_old, ONLY : fisrtilp
     USE lmdz_wake_ini, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
@@ -84,5 +85,5 @@
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
     USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
-    USE lscp_ini_mod, ONLY : lscp_ini
+    USE lmdz_lscp_ini, ONLY : lscp_ini
     USE lmdz_ratqs_main, ONLY : ratqs_main
     USE lmdz_ratqs_ini, ONLY : ratqs_ini
@@ -885,5 +886,4 @@
     !KE43
     EXTERNAL conema3  ! convect4.3
-    EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
     !AA 
     ! JBM (3/14) fisrtilp_tr not loaded
