Index: LMDZ5/trunk/libf/dyn3d_common/infotrac.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/infotrac.F90	(revision 2085)
+++ LMDZ5/trunk/libf/dyn3d_common/infotrac.F90	(revision 2086)
@@ -5,4 +5,6 @@
 ! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
   INTEGER, SAVE :: nqtot
+!CR: on ajoute le nombre de traceurs de l eau
+  INTEGER, SAVE :: nqo
 
 ! nbtr : number of tracers not including higher order of moment or water vapor or liquid
@@ -228,4 +230,11 @@
          endif ! of if (planet_type=="earth")
        END IF
+
+!CR: nombre de traceurs de l eau
+       if (tnom_0(3) == 'H2Oi') then
+          nqo=3
+       else
+          nqo=2
+       endif
        
        WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
Index: LMDZ5/trunk/libf/phylmd/add_pbl_tend.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/add_pbl_tend.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/add_pbl_tend.F90	(revision 2086)
@@ -1,3 +1,3 @@
-SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, paprs, text)
+SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text)
   ! ======================================================================
   ! Ajoute les tendances de couche limite, soit determinees par la
@@ -27,5 +27,5 @@
   ! ------------
   REAL zdu(klon, klev), zdv(klon, klev)
-  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev)
+  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
   CHARACTER *(*) text
   REAL paprs(klon,klev+1)
@@ -46,7 +46,7 @@
     PRINT *, ' add_pbl_tend, zzdt ', zzdt
     PRINT *, ' add_pbl_tend, zzdq ', zzdq
-    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, paprs, text)
+    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text)
   ELSE
-    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, paprs, text)
+    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text)
   END IF
 
Index: LMDZ5/trunk/libf/phylmd/add_phys_tend.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/add_phys_tend.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/add_phys_tend.F90	(revision 2086)
@@ -2,5 +2,5 @@
 ! $Id$
 !
-SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,paprs,text)
+SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text)
 !======================================================================
 ! Ajoute les tendances des variables physiques aux variables 
@@ -25,5 +25,5 @@
 !------------
 REAL zdu(klon,klev),zdv(klon,klev)
-REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev)
+REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
 REAL paprs(klon,klev+1)
 CHARACTER*(*) text
@@ -66,4 +66,5 @@
      v_seri(:,:)=v_seri(:,:)+zdv(:,:)
      ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
+     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
 
 !======================================================================
Index: LMDZ5/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2086)
@@ -4,5 +4,5 @@
 !
 SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
-     d_t, d_q, d_ql, rneb, radliq, rain, snow,          &
+     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
      pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
      frac_impa, frac_nucl, beta,                        &
@@ -41,4 +41,5 @@
   REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
   REAL d_ql(klon,klev) ! incrementation de l'eau liquide
+  REAL d_qi(klon,klev) ! incrementation de l'eau glace
   REAL rneb(klon,klev) ! fraction nuageuse
   REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
@@ -107,5 +108,6 @@
   REAL DDT0
   PARAMETER (DDT0=.01)
-  INTEGER n_i, iter
+  INTEGER n_i(klon), iter
+  REAL cste
  
   REAL zrfl(klon), zrfln(klon), zqev, zqevt 
@@ -123,4 +125,5 @@
   REAL zmelt, zpluie, zice, zcondold
   PARAMETER (ztfondue=278.15)
+  REAL dzfice(klon)
   !
   LOGICAL appel1er
@@ -165,6 +168,7 @@
   DATA appel1er /.TRUE./
   !ym
-
-  ice_thermo = iflag_ice_thermo .GE. 1
+!CR: pour iflag_ice_thermo=2, on active que la convection
+!  ice_thermo = iflag_ice_thermo .GE. 1
+   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
   zdelq=0.0
 
@@ -209,4 +213,5 @@
   !  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
@@ -220,4 +225,5 @@
       exposant_glace_old = 6
     ENDIF
+    
   ENDIF
   
@@ -243,4 +249,5 @@
         d_q(i,k) = 0.0
         d_ql(i,k) = 0.0
+        d_qi(i,k) = 0.0
         rneb(i,k) = 0.0
         radliq(i,k) = 0.0
@@ -451,4 +458,17 @@
          zifl(i) = zifln(i)
 
+!CR ATTENTION: deplacement de la fonte de la glace
+           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'
+!fin CR
+
+
+
            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
         ENDDO
@@ -487,9 +507,9 @@
      !
 !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
+!       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
 
      IF (cpartiel) THEN
@@ -550,8 +570,39 @@
            end if
 
+!CR: variation de qsat avec T en présence de glace ou non
+!initialisations
            do i=1,klon
+              DT(i) = 0.
+              n_i(i)=0
               Tbef(i)=zt(i)
-              qlbef(i)=0. 
-              if (lognormale(i)) then
+              qlbef(i)=0.
+           enddo
+
+
+!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
+           if (iflag_fisrtilp_qsat.ge.0) then
+             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(i)=abs(DT(i)).gt.DDT0
+                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
+                 Tbef(i)=Tbef(i)+DT(i) 
+                 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.eq.1) then
+                    zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
+                    endif
+                 endif
+                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+                 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))
@@ -574,60 +625,75 @@
                  endif
 
+                 endif !convergence
+               enddo ! boucle en i
+
+                 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))
                  num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
                  denom=1.+rneb(i,k)*zdqs(i)
                  DT(i)=num/denom
-              endif
-           enddo
+                 n_i(i)=n_i(i)+1
+                 endif
+                 enddo
+
+                 else
+
+!calcul de la fraction de glace
+!CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
+!                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
+!                                     t_glace_max, exposant_glace)
+!                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
+!                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+!                 zfice(i) = zfice(i)**nexpo
+                 if (iflag_t_glace.eq.1) then
+                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
+                      t_glace_min,t_glace_max,exposant_glace,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.eq.1) 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))
+                 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)
+                 DT(i)=num/denom
+                 n_i(i)=n_i(i)+1
+
+                 endif ! fin convergence
+                 enddo ! fin boucle i
+
+                 endif !ice_thermo
+
+!                 endif
+!               enddo
   
-           n_i=1
-!               do while ((abs(DT(i)).gt.DDT0).and.((zqn(i)-zqs(i)).gt.0.))
-! iflag_fisrtilp_qsat=0: qsat ne varie pas avec T
-! iflag_fisrtilp_qsat > 1 : nombre d iterations max pour converger sur le calcul de qsat(T)
-!               do while (n_i.le.iflag_fisrtilp_qsat)
-           if (iflag_fisrtilp_qsat.ge.1) then
-           do iter=1,iflag_fisrtilp_qsat
-              do i=1,klon
-                 convergence(i)=abs(DT(i)).gt.DDT0
-                 if (convergence(i).and.lognormale(i)) then
-                 Tbef(i)=Tbef(i)+DT(i)                 
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(i)))
-                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
-                 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
-
-                 qlbef(i)=max(0.,zqn(i)-zqs(i))    
-                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
-                 denom=1.+rneb(i,k)*zdqs(i)
-                 DT(i)=num/denom
-                 n_i=n_i+1
-                 endif
-              enddo
-              enddo
-
+          
+             enddo
            endif
 
@@ -635,7 +701,9 @@
         endif ! iflag_pdf
 
-        if (iflag_fisrtilp_qsat.eq.-1) then 
-!CR: ATTENTION option fausse mais a existe
-        DO i=1,klon
+
+!        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
@@ -653,6 +721,7 @@
            ENDIF
         ENDDO
-
-        ELSE
+        ENDIF
+
+!        ELSE
 
         DO i=1,klon
@@ -674,5 +743,5 @@
 
 
-        ENDIF
+!        ENDIF
 
         !         do i=1,klon
@@ -723,11 +792,18 @@
         endif
      ELSE
-       IF (iflag_t_glace.EQ.0) THEN
+         if (iflag_t_glace.eq.1) then
+            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
+                 t_glace_min,t_glace_max,exposant_glace,zfice(:))
+         endif
          if (iflag_fisrtilp_qsat.lt.1) then
            DO i = 1, klon
-              zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-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
+! JBM: icefrac_lsc is now a function contained in microphys_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 - (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
+              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))
@@ -735,38 +811,22 @@
          else
            DO i=1, klon
-              zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-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
-!CR: ATTENTION zt different de Tbef: à corriger
+              if (lognormale(i)) then
+                 zt(i)=Tbef(i)
+              else
+! JBM: icefrac_lsc is now a function contained in microphys_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 - (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
+              endif
               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)))
-           ENDDO
-         endif
-!         print*,zt(i),zrfl(i),zifl(i),'temp1'
-       ELSE ! of IF (iflag_t_glace.EQ.0)
-         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
-&               t_glace_min,t_glace_max,exposant_glace,zfice(:))
-         if (iflag_fisrtilp_qsat.lt.1) then
-           DO i = 1, klon
-! JBM: icefrac_lsc is now a function contained in microphys_mod
-!             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                    t_glace_max, exposant_glace)
-              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 microphys_mod
-!             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                    t_glace_max, exposant_glace)
-!CR: ATTENTION zt different de Tbef: à corriger
-              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)))
+              endif
            ENDDO
          endif
 !         print*,zt(i),zrfl(i),zifl(i),'temp1'
        ENDIF
-     ENDIF
 !>AJ
      !
@@ -794,12 +854,12 @@
        ELSE ! of IF (iflag_t_glace.EQ.0)
          CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
-&               t_glace_min,t_glace_max,exposant_glace,zfice(:))
-!        DO i = 1, klon
-!           IF (rneb(i,k).GT.0.0) THEN
+               t_glace_min,t_glace_max,exposant_glace,zfice(:))
+!         DO i = 1, klon
+!            IF (rneb(i,k).GT.0.0) THEN
 ! JBM: icefrac_lsc is now a function contained in microphys_mod
-!             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
-!                                    t_glace_max, exposant_glace)
-!           ENDIF
-!        ENDDO
+!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
+!                                     t_glace_max, exposant_glace)
+!            ENDIF
+!         ENDDO
        ENDIF
      ENDIF
@@ -874,5 +934,12 @@
        DO i = 1, klon
          IF (rneb(i,k).GT.0.0) THEN
+!CR on prend en compte la phase glace
+           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) &
@@ -889,16 +956,17 @@
      ENDIF
 
-     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)
+!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))
+!           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
+!       ENDDO
+!     ENDIF
 
        
@@ -1003,8 +1071,10 @@
   ! 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)
         snow(i) = zrfl(i)+zifl(i)
 !>AJ
@@ -1015,7 +1085,17 @@
      ENDIF
   ENDDO
+
+  ELSE
+     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
@@ -1026,4 +1106,5 @@
      END DO
   END DO
+  ENDIF
   !
 
Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2086)
@@ -33,5 +33,5 @@
       REAL, SAVE, ALLOCATABLE :: d_t_wake(:,:),d_q_wake(:,:)
       !$OMP THREADPRIVATE( d_t_wake,d_q_wake)
-      REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:)
+      REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:),d_qi_lsc(:,:)
       !$OMP THREADPRIVATE(d_t_lsc,d_q_lsc,d_ql_lsc)
       REAL, SAVE, ALLOCATABLE :: d_t_ajsb(:,:), d_q_ajsb(:,:)
@@ -307,5 +307,5 @@
       allocate(d_t_wake(klon,klev),d_q_wake(klon,klev))
       allocate(d_t_lsc(klon,klev),d_q_lsc(klon,klev))
-      allocate(d_ql_lsc(klon,klev))
+      allocate(d_ql_lsc(klon,klev),d_qi_lsc(klon,klev))
       allocate(d_t_ajsb(klon,klev),d_q_ajsb(klon,klev))
       allocate(d_t_ajs(klon,klev),d_q_ajs(klon,klev))
@@ -467,5 +467,5 @@
       deallocate(d_t_wake,d_q_wake)
       deallocate(d_t_lsc,d_q_lsc)
-      deallocate(d_ql_lsc)
+      deallocate(d_ql_lsc,d_qi_lsc)
       deallocate(d_t_ajsb,d_q_ajsb)
       deallocate(d_t_ajs,d_q_ajs)
Index: LMDZ5/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 2086)
@@ -399,82 +399,85 @@
 #endif
 
-      IF (nqtot>=3) THEN
-            DO iq=3,nqtot  
+!CR: ajout d'une variable eau
+!      IF (nqtot>=3) THEN
+       IF (nqtot>=nqo+1) THEN
+!            DO iq=3,nqtot  
+            DO iq=nqo+1,nqtot 
             iiq=niadv(iq)
-            o_trac(iq-2) = ctrl_out((/ 4, 5, 1, 1, 1, 10, 11, 11, 11 /), &
+            o_trac(iq-nqo) = ctrl_out((/ 4, 5, 1, 1, 1, 10, 11, 11, 11 /), &
                            tname(iiq),'Tracer '//ttext(iiq), "-",  &
                            (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_vdf(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_vdf(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                               'd'//trim(tname(iq))//'_vdf',  &
                               'Tendance tracer '//ttext(iiq), "-" , &
                               (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_the(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_the(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                               'd'//trim(tname(iq))//'_the', &
                               'Tendance tracer '//ttext(iiq), "-", &
                               (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_con(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_con(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                               'd'//trim(tname(iq))//'_con', &
                               'Tendance tracer '//ttext(iiq), "-", &
                               (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_lessi_impa(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_lessi_impa(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                      'd'//trim(tname(iq))//'_lessi_impa', &
                                      'Tendance tracer '//ttext(iiq), "-", &
                                      (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_lessi_nucl(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_lessi_nucl(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                      'd'//trim(tname(iq))//'_lessi_nucl', &
                                      'Tendance tracer '//ttext(iiq), "-", &
                                      (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_insc(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_insc(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                'd'//trim(tname(iq))//'_insc', &
                                'Tendance tracer '//ttext(iiq), "-", &
                                (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_bcscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_bcscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                  'd'//trim(tname(iq))//'_bcscav', &
                                  'Tendance tracer '//ttext(iiq), "-", &
                                  (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_evapls(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_evapls(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                  'd'//trim(tname(iq))//'_evapls', &
                                  'Tendance tracer '//ttext(iiq), "-", &
                                  (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_ls(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_ls(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                              'd'//trim(tname(iq))//'_ls', &
                              'Tendance tracer '//ttext(iiq), "-", &
                              (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_trsp(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_trsp(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                'd'//trim(tname(iq))//'_trsp', &
                                'Tendance tracer '//ttext(iiq), "-", &
                                (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_sscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_sscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                 'd'//trim(tname(iq))//'_sscav', &
                                 'Tendance tracer '//ttext(iiq), "-", &
                                 (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_sat(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_sat(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                'd'//trim(tname(iq))//'_sat', &
                                'Tendance tracer '//ttext(iiq), "-", &
                                (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_uscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_uscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                                 'd'//trim(tname(iq))//'_uscav', &
                                 'Tendance tracer '//ttext(iiq), "-", &
                                  (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_dtr_dry(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
+            o_dtr_dry(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
                               'cum'//'d'//trim(tname(iq))//'_dry', &
                               'tracer tendency dry deposition'//ttext(iiq), "-", &
                               (/ '', '', '', '', '', '', '', '', '' /))
 
-            o_trac_cum(iq-2) = ctrl_out((/ 3, 4, 10, 10, 10, 10, 11, 11, 11 /), &
+            o_trac_cum(iq-nqo) = ctrl_out((/ 3, 4, 10, 10, 10, 10, 11, 11, 11 /), &
                                'cum'//tname(iiq),&
                                'Cumulated tracer '//ttext(iiq), "-", &
Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2086)
@@ -222,5 +222,5 @@
     USE ocean_slab_mod, only: tslab, slab_bils
     USE indice_sol_mod, only: nbsrf
-    USE infotrac, only: nqtot
+    USE infotrac, only: nqtot, nqo
     USE comgeomphy, only: airephy
     USE surface_data, only: type_ocean, ok_veget, ok_snow
@@ -1288,20 +1288,22 @@
        ENDDO !nfiles
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-       IF (nqtot.GE.3) THEN
-          DO iq=3,nqtot
-             CALL histwrite_phy(o_trac(iq-2), qx(:,:,iq))
-             CALL histwrite_phy(o_dtr_vdf(iq-2),d_tr_cl(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_the(iq-2),d_tr_th(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_con(iq-2),d_tr_cv(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_lessi_impa(iq-2),d_tr_lessi_impa(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_lessi_nucl(iq-2),d_tr_lessi_nucl(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_insc(iq-2),d_tr_insc(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_bcscav(iq-2),d_tr_bcscav(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_evapls(iq-2),d_tr_evapls(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_ls(iq-2),d_tr_ls(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_trsp(iq-2),d_tr_trsp(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_sscav(iq-2),d_tr_sscav(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_sat(iq-2),d_tr_sat(:,:,iq-2))
-             CALL histwrite_phy(o_dtr_uscav(iq-2),d_tr_uscav(:,:,iq-2))
+!       IF (nqtot.GE.3) THEN
+!          DO iq=3,nqtot
+          IF (nqtot.GE.nqo+1) THEN
+            DO iq=nqo+1,nqtot
+             CALL histwrite_phy(o_trac(iq-nqo), qx(:,:,iq))
+             CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_lessi_impa(iq-nqo),d_tr_lessi_impa(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_lessi_nucl(iq-nqo),d_tr_lessi_nucl(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_insc(iq-nqo),d_tr_insc(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_bcscav(iq-nqo),d_tr_bcscav(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_evapls(iq-nqo),d_tr_evapls(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_ls(iq-nqo),d_tr_ls(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_trsp(iq-nqo),d_tr_trsp(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_sscav(iq-nqo),d_tr_sscav(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_sat(iq-nqo),d_tr_sat(:,:,iq-nqo))
+             CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo))
              zx_tmp_fi2d=0.
              IF(vars_defined) THEN
@@ -1310,5 +1312,5 @@
                 ENDDO
              ENDIF
-             CALL histwrite_phy(o_trac_cum(iq-2), zx_tmp_fi2d)
+             CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
           ENDDO
        ENDIF
Index: LMDZ5/trunk/libf/phylmd/physiq.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq.F90	(revision 2085)
+++ LMDZ5/trunk/libf/phylmd/physiq.F90	(revision 2086)
@@ -201,5 +201,7 @@
   INTEGER iliq          ! indice de traceurs pour eau liquide
   PARAMETER (iliq=2)
-
+!CR: on ajoute la phase glace
+  INTEGER isol          ! indice de traceurs pour eau glace
+  PARAMETER (isol=3)
   !
   !
@@ -596,5 +598,5 @@
 
   ! tendance nulles
-  REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0
+  REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
 
   !
@@ -1377,4 +1379,5 @@
   dq0(:,:)=0.
   dql0(:,:)=0.
+  dqi0(:,:)=0.
   !
   ! Mettre a zero des variables de sortie (pour securite)
@@ -1422,13 +1425,22 @@
         q_seri(i,k)  = qx(i,k,ivap)
         ql_seri(i,k) = qx(i,k,iliq)
-        qs_seri(i,k) = 0.
+!CR: ATTENTION, on rajoute la variable glace
+        if (nqo.eq.2) then
+           qs_seri(i,k) = 0.
+        else if (nqo.eq.3) then
+           qs_seri(i,k) = qx(i,k,isol)
+        endif
      ENDDO
   ENDDO
   tke0(:,:)=pbl_tke(:,:,is_ave)
-  IF (nqtot.GE.3) THEN
-     DO iq = 3, nqtot
+!CR:Nombre de traceurs de l'eau: nqo
+!  IF (nqtot.GE.3) THEN
+   IF (nqtot.GE.(nqo+1)) THEN
+!     DO iq = 3, nqtot        
+     DO iq = nqo+1, nqtot  
         DO  k = 1, klev
            DO  i = 1, klon
-              tr_seri(i,k,iq-2) = qx(i,k,iq)
+!              tr_seri(i,k,iq-2) = qx(i,k,iq)
+              tr_seri(i,k,iq-nqo) = qx(i,k,iq)
            ENDDO
         ENDDO
@@ -1617,4 +1629,7 @@
         ENDIF
         !>jyg
+     
+        if (iflag_ice_thermo.eq.0) then   
+!pas necessaire a priori
 
         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
@@ -1627,4 +1642,25 @@
         d_t_eva(i,k) = za
         d_q_eva(i,k) = zb
+
+        else
+
+!CR: on ré-évapore eau liquide et glace
+
+!        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
+!        zb = MAX(0.0,ql_seri(i,k))
+!        za = - MAX(0.0,ql_seri(i,k)) &
+!             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
+        zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
+        za = - MAX(0.0,ql_seri(i,k))*zlvdcp & 
+             - MAX(0.0,qs_seri(i,k))*zlsdcp
+        t_seri(i,k) = t_seri(i,k) + za
+        q_seri(i,k) = q_seri(i,k) + zb
+        ql_seri(i,k) = 0.0
+!on évapore la glace
+        qs_seri(i,k) = 0.0
+        d_t_eva(i,k) = za
+        d_q_eva(i,k) = zb
+        endif
+
      ENDDO
   ENDDO
@@ -1764,8 +1800,8 @@
      IF (klon_glo==1) THEN
         CALL add_pbl_tend &
-             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
+             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
      ELSE
         CALL add_phys_tend &
-             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
+             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
      ENDIF
      !--------------------------------------------------------------------
@@ -2135,5 +2171,5 @@
   !-----------------------------------------------------------------------------------------
   ! ajout des tendances de la diffusion turbulente
-  CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,paprs,'con')
+  CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,dqi0,paprs,'con')
   !-----------------------------------------------------------------------------------------
 
@@ -2252,5 +2288,5 @@
      d_t_wake(:,:)=dt_wake(:,:)*dtime
      d_q_wake(:,:)=dq_wake(:,:)*dtime
-     CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,paprs,'wake')
+     CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
      !-----------------------------------------------------------------------------------------
 
@@ -2512,5 +2548,5 @@
         !-----------------------------------------------------------------------------------------
         ! ajout des tendances de l'ajustement sec ou des thermiques
-        CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,paprs,'ajsb')
+        CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
@@ -2558,5 +2594,5 @@
   CALL fisrtilp(dtime,paprs,pplay, &
        t_seri, q_seri,ptconv,ratqs, &
-       d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
+       d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
        rain_lsc, snow_lsc, &
        pfrac_impa, pfrac_nucl, pfrac_1nucl, &
@@ -2570,9 +2606,10 @@
   !-----------------------------------------------------------------------------------------
   ! ajout des tendances de la diffusion turbulente
-  CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,paprs,'lsc')
+  CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
   !-----------------------------------------------------------------------------------------
   DO k = 1, klev
      DO i = 1, klon
         cldfra(i,k) = rneb(i,k)
+!CR: a quoi ca sert? Faut-il ajouter qs_seri?
         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
      ENDDO
@@ -2845,5 +2882,9 @@
         zx_t = t_seri(i,k)
         IF (thermcep) THEN
+           if (iflag_ice_thermo.eq.0) then 
            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
+           else
+           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
+           endif
            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
            zx_qs  = MIN(0.5,zx_qs)
@@ -3288,5 +3329,5 @@
      !-----------------------------------------------------------------------------------------
      ! ajout des tendances de la trainee de l'orographie
-     CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,paprs,'oro')
+     CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
      !-----------------------------------------------------------------------------------------
      !
@@ -3334,5 +3375,5 @@
      !-----------------------------------------------------------------------------------------
      ! ajout des tendances de la portance de l'orographie
-     CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,paprs,'lif')
+     CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
      !-----------------------------------------------------------------------------------------
      !
@@ -3348,5 +3389,5 @@
      !
      !  ajout des tendances
-     CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,paprs,'hin')
+     CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
 
   ENDIF
@@ -3356,5 +3397,5 @@
           rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
           du_gwd_rando, dv_gwd_rando)
-     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,paprs, &
+     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
           'flott_gwd_rando')
   end if
@@ -3623,12 +3664,20 @@
         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
+!CR: on ajoute le contenu en glace
+        if (nqo.eq.3) then
+        d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
+        endif
      ENDDO
   ENDDO
   !
-  IF (nqtot.GE.3) THEN
-     DO iq = 3, nqtot
+!CR: nb de traceurs eau: nqo
+!  IF (nqtot.GE.3) THEN
+   IF (nqtot.GE.(nqo+1)) THEN
+!     DO iq = 3, nqtot
+     DO iq = nqo+1, nqtot
         DO  k = 1, klev
            DO  i = 1, klon
-              d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
+!              d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
+               d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
            ENDDO
         ENDDO
@@ -3653,9 +3702,13 @@
 
 !!! RomP >>>
-  IF (nqtot.GE.3) THEN
-     DO iq = 3, nqtot
+!CR: nb de traceurs eau: nqo
+!  IF (nqtot.GE.3) THEN
+   IF (nqtot.GE.(nqo+1)) THEN
+!     DO iq = 3, nqtot
+     DO iq = nqo+1, nqtot
         DO k = 1, klev
            DO i = 1, klon
-              tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
+!              tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
+              tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
            ENDDO
         ENDDO
