Index: LMDZ5/trunk/libf/phylmd/clesphys.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/clesphys.h	(revision 1849)
@@ -74,5 +74,5 @@
        REAL freq_COSP
        LOGICAL :: ok_cosp,ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP
-       INTEGER :: ip_ebil_phy, iflag_rrtm
+       INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo
        LOGICAL :: ok_strato
        LOGICAL :: ok_hines
@@ -102,5 +102,5 @@
      &     , ok_lic_melt, cvl_corr, aer_type                            &
      &     , qsol0, iflag_rrtm, ok_strato,ok_hines,ecrit_LES            &
-     &     , co2_ppm0
+     &     , co2_ppm0, iflag_ice_thermo
      
        save /clesphys/
Index: LMDZ5/trunk/libf/phylmd/concvl.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/concvl.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/concvl.F	(revision 1849)
@@ -1,3 +1,3 @@
-      SUBROUTINE concvl (iflag_con,iflag_clos,
+      SUBROUTINE concvl (iflag_clos,
      .             dtime,paprs,pplay,
      .             t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
@@ -70,7 +70,9 @@
 c======================================================================
 c
+
+#include "clesphys.h"
 #include "dimensions.h"
 c
-       INTEGER iflag_con,iflag_clos
+       INTEGER iflag_clos
 c
        REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
@@ -376,5 +378,6 @@
       nloc=klon
       CALL cva_driver(klon,klev,klev+1,ntra,nloc,
-     $              iflag_con,iflag_mix,iflag_clos,dtime,
+     $              iflag_con,iflag_mix,iflag_ice_thermo,
+     $              iflag_clos,dtime,
      :              t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
      $              em_p,em_ph,
Index: LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 1849)
@@ -141,4 +141,5 @@
   LOGICAL,SAVE :: reevap_ice_omp
   INTEGER,SAVE :: iflag_pdf_omp
+  INTEGER,SAVE :: iflag_ice_thermo_omp
   REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
   REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
@@ -966,4 +967,13 @@
   t_glace_max_omp = 273.13
   call getin('t_glace_max',t_glace_max_omp)
+
+!
+!Config Key  = iflag_ice_thermo
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  iflag_ice_thermo_omp = 0
+  call getin('iflag_ice_thermo',iflag_ice_thermo_omp)
 
 !Config Key  = rei_min
@@ -1694,4 +1704,5 @@
     t_glace_min = t_glace_min_omp
     t_glace_max = t_glace_max_omp
+    iflag_ice_thermo = iflag_ice_thermo_omp
     rei_min = rei_min_omp
     rei_max = rei_max_omp
@@ -1911,4 +1922,5 @@
   write(lunout,*)' t_glace_min = ',t_glace_min
   write(lunout,*)' t_glace_max = ',t_glace_max
+  write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
   write(lunout,*)' rei_min = ',rei_min
   write(lunout,*)' rei_max = ',rei_max
Index: LMDZ5/trunk/libf/phylmd/cv3_routines.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/cv3_routines.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cv3_routines.F	(revision 1849)
@@ -143,5 +143,5 @@
 
       SUBROUTINE cv3_prelim(len,nd,ndp1,t,q,p,ph
-     :                    ,lv,cpn,tv,gz,h,hm,th)
+     :                    ,lv,lf,cpn,tv,gz,h,hm,th)
       implicit none
 
@@ -157,5 +157,5 @@
 
 c outputs:
-      real lv(len,nd), cpn(len,nd), tv(len,nd)
+      real lv(len,nd), lf(len,nd), cpn(len,nd), tv(len,nd)
       real gz(len,nd), h(len,nd), hm(len,nd)
       real th(len,nd)
@@ -178,4 +178,5 @@
 cdebug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
           lv(i,k)= lv0-clmcpv*(t(i,k)-273.15)
+          lf(i,k)= lf0-clmci*(t(i,k)-273.15)
           cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
           cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
@@ -918,8 +919,40 @@
       end
 
+      SUBROUTINE Icefrac(t,clw,qi,nl,len)
+      implicit none
+
+
+cJAM--------------------------------------------------------------------
+C    Calcul de la quantité d'eau sous forme de glace
+C--------------------------------------------------------------------
+      Real qi(len,nl)
+      Real t(len,nl), clw(len,nl)
+      Real fracg
+      Integer nl, len, k, i
+
+      do k=3, nl
+           do i = 1, len
+           if (t(i,k).gt.263.15) then
+           qi(i,k)=0.
+           else
+           if (t(i,k).lt.243.15) then
+           qi(i,k)=clw(i,k)
+           else
+           fracg=(263.15-t(i,k))/20
+           qi(i,k)=clw(i,k)*fracg
+           endif
+           endif
+C           print*,t(i,k),qi(i,k),'temp,testglace'
+           enddo
+       enddo
+
+      return
+
+      end
+
       SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk
      :                       ,tnk,qnk,gznk,hnk,t,q,qs,gz
-     :                       ,p,h,tv,lv,pbase,buoybase,plcl
-     o                       ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
+     :                       ,p,h,tv,lv,lf,pbase,buoybase,plcl
+     o                       ,inb,tp,tvp,clw,hp,ep,sigp,buoy,frac)
       implicit none
 
@@ -945,7 +978,8 @@
 #include "cv3param.h"
 #include "conema3.h"
+#include "cvflag.h"
 
 c inputs:
-      integer ncum, nd, nloc
+      integer ncum, nd, nloc, j
       integer icb(nloc), icbs(nloc), nk(nloc)
       real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
@@ -953,5 +987,5 @@
       real tnk(nloc), qnk(nloc), gznk(nloc)
       real hnk(nloc)
-      real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
+      real lv(nloc,nd), lf(nloc,nd), tv(nloc,nd), h(nloc,nd)
       real pbase(nloc), buoybase(nloc), plcl(nloc)
 
@@ -964,9 +998,13 @@
 c local variables:
       integer i, k
-      real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
-      real by, defrac, pden
+       real tg,qg,ahg,alv,alf,s,tc,es,esi,denom,rg,tca,elacrit
+      real als 
+      real qsat_new,snew, qi(nloc,nd)
+      real by, defrac, pden, tbis
       real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
+      real frac(nloc,nd)
       logical lcape(nloc)
       integer iposit(nloc)
+      Real fracg
 
 !=====================================================================
@@ -978,4 +1016,5 @@
        ep(i,k)=0.0
        sigp(i,k)=spfac
+       qi(i,k)=0.
  160  continue
  170  continue
@@ -1060,6 +1099,11 @@
 
 c convect3: no approximation:
-           tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
-
+           if (cvflag_ice) then
+              tp(i,k)=Max(0.,(ah0(i)-gz(i,k)-alv*qg)
+     &                /(cpd+(cl-cpd)*qnk(i)))
+           else
+              tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
+           endif
+c          
                clw(i,k)=qnk(i)-qg
                clw(i,k)=max(0.0,clw(i,k))
@@ -1068,5 +1112,83 @@
 c convect3: (qg utilise au lieu du vrai mixing ratio rg):
                tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing
+               if (cvflag_ice) then
+               if(clw(i,k).lt.1.e-11) then
+               tp(i,k)=tv(i,k)
+               tvp(i,k)=tv(i,k)
+               endif
+               endif
             endif
+
+        if (cvflag_ice) then
+cCR:attention boucle en klon dans Icefrac          
+c          Call Icefrac(t,clw,qi,nl,nloc)
+           if (t(i,k).gt.263.15) then
+           qi(i,k)=0.
+           else
+           if (t(i,k).lt.243.15) then
+           qi(i,k)=clw(i,k)
+           else
+           fracg=(263.15-t(i,k))/20
+           qi(i,k)=clw(i,k)*fracg
+           endif
+           endif     
+cCR: fin test     
+          if(t(i,k).lt.263.15) then
+cCR: on commente les calculs d'Arnaud car division par zero
+cnouveau calcul propose par JYG
+c           alv=lv0-clmcpv*(t(i,k)-273.15)
+c           alf=lf0-clmci*(t(i,k)-273.15)
+c           tg=tp(i,k)
+c           tc=tp(i,k)-273.15
+c           denom=243.5+tc
+c             do j=1,3
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c         il faudra que esi vienne en argument de la convection
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c             tbis=t(i,k)+(tp(i,k)-tg)
+c             esi=exp(23.33086-(6111.72784/tbis)
+c     :               +0.15215*log(tbis))
+c             qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
+c             snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/
+c     :                               (rrv*tbis*tbis)
+c             snew=1./snew
+c             print*,esi,qsat_new,snew,'esi,qsat,snew'
+c             tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
+c             print*,k,tp(i,k),qnk(i),'avec glace'
+c             print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
+c             enddo
+
+           alv=lv0-clmcpv*(t(i,k)-273.15)
+           alf=lf0+clmci*(t(i,k)-273.15)
+           als=alf+alv
+           tg=tp(i,k)
+           tp(i,k) = t(i,k)
+           do j=1,3
+             esi=exp(23.33086-(6111.72784/tp(i,k))
+     :               +0.15215*log(tp(i,k)))
+             qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
+             snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*als*qsat_new/
+     :                               (rrv*tp(i,k)*tp(i,k))
+             snew=1./snew
+cc             print*,esi,qsat_new,snew,'esi,qsat,snew'
+             tp(i,k)=tp(i,k)+
+     :         ( (cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) +
+     :         alv*(qg-qsat_new) + alf*qi(i,k) )*snew
+c             print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k),
+c     :             'k,tp,q,qt,qi avec glace'
+           enddo
+
+cCR:reprise du code AJ
+           clw(i,k)=qnk(i)-qsat_new
+           clw(i,k)=max(0.0,clw(i,k))
+           tvp(i,k)=Max(0.,tp(i,k)*(1.+qsat_new/eps-qnk(i)))
+c           print*,tvp(i,k),'tvp'
+          endif
+          if(clw(i,k).lt.1.e-11) then
+           tp(i,k)=tv(i,k)
+           tvp(i,k)=tv(i,k)
+          endif
+        endif  ! (cvflag_ice)
+
   290     continue
   300   continue
@@ -1303,5 +1425,15 @@
         do 590 i=1,ncum
         if((k.ge.icb(i)).and.(k.le.inb(i)))then
+
+          if (cvflag_ice) then
+          frac(i,k)=1.-(t(i,k)-243.15)/(263.15-243.15)
+          frac(i,k)=min(max(frac(i,k),0.0),1.0)
+          hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))
+     :            *ep(i,k)*clw(i,k)
+
+          else
           hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
+          endif
+
         endif
  590    continue
@@ -1555,5 +1687,5 @@
 
       SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
-     :                    ,ph,t,rr,rs,u,v,tra,h,lv,qnk
+     :                    ,ph,t,rr,rs,u,v,tra,h,lv,lf,frac,qnk
      :                    ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
      :   ,ment,qent,uent,vent,nent,sij,elij,ments,qents,traent)
@@ -1567,4 +1699,5 @@
 #include "cvthermo.h"
 #include "cv3param.h"
+#include "cvflag.h"
 
 c inputs:
@@ -1578,4 +1711,5 @@
       real tra(nloc,nd,ntra) ! input of convect3
       real lv(nloc,na), h(nloc,na), hp(nloc,na)
+      real lf(nloc,na), frac(nloc,na)
       real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
       real m(nloc,na)        ! input of convect3
@@ -1659,4 +1793,14 @@
           rti=qnk(il)-ep(il,i)*clw(il,i)
           bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
+
+
+          if (cvflag_ice) then
+c          print*,cvflag_ice,'cvflag_ice dans do 700'
+          if (t(il,j).le.263.15) then
+          bf2=1.+(lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)
+     :          *lf(il,j))*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
+          endif
+          endif
+
           anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
           denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
@@ -1671,6 +1815,14 @@
           if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
      :                 .and.j.gt.i)then
+
+          if (cvflag_ice) then
+          anum=anum-(lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)
+     :          -cwat*bf2)
+          denom=denom+(lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
+          else
            anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
            denom=denom+lv(il,j)*(rr(il,i)-rti)
+          endif
+
            if(abs(denom).lt.0.01)denom=0.01
            sij(il,i,j)=anum/denom
@@ -1780,8 +1932,19 @@
         lwork(il)=(nent(il,i).ne.0)
         qp=qnk(il)-ep(il,i)*clw(il,i)
+
+        if(cvflag_ice) then
+
+        anum=h(il,i)-hp(il,i)-(lv(il,i)+frac(il,i)*lf(il,i))*
+     :       (qp-rs(il,i))+(cpv-cpd)*t(il,i)*(qp-rr(il,i))
+        denom=h(il,i)-hp(il,i)+(lv(il,i)+frac(il,i)*lf(il,i))*
+     :        (rr(il,i)-qp)+(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
+        else
+
         anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
      :           +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
         denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
      :           +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
+        endif
+
         if(abs(denom).lt.0.01)denom=0.01
         scrit(il)=anum/denom
@@ -1948,7 +2111,8 @@
       SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb,iflag
      :              ,t,rr,rs,gz,u,v,tra,p,ph
-     :              ,th,tv,lv,cpn,ep,sigp,clw
+     :              ,th,tv,lv,lf,cpn,ep,sigp,clw
      :              ,m,ment,elij,delt,plcl,coef_clos
-     o              ,mp,rp,up,vp,trap,wt,water,evap,b,sigd
+     o              ,mp,rp,up,vp,trap,wt,water,evap,fondue,ice
+     :              ,faci,b,sigd
      o              ,wdtrainA,wdtrainM)                                ! RomP
       implicit none
@@ -1969,4 +2133,5 @@
       real ep(nloc,na), sigp(nloc,na), clw(nloc,na)
       real th(nloc,na),tv(nloc,na),lv(nloc,na),cpn(nloc,na)
+      real lf(nloc,na)
       real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
       real coef_clos(nloc)
@@ -1978,4 +2143,5 @@
       real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
       real water(nloc,na), evap(nloc,na), wt(nloc,na)
+      real ice(nloc,na), fondue(nloc,na),faci(nloc,na)
       real trap(nloc,na,ntra)
       real b(nloc,na), sigd(nloc)
@@ -1988,12 +2154,14 @@
 c local variables
       integer i,j,k,il,num1,ndp1
-      real tinv, delti
+      real tinv, delti, coef
       real awat, afac, afac1, afac2, bfac
-      real pr1, pr2, sigt, b6, c6, revap, delth
+      real pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
       real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
-      real ampmax
+      real ampmax, thaw
       real tevap(nloc)
-      real lvcp(nloc,na)
+      real lvcp(nloc,na),lfcp(nloc,na)
       real h(nloc,na),hm(nloc,na)
+      real frac(nloc,na)
+      real fraci(nloc,na),prec(nloc,na)
       real wdtrain(nloc)
       logical lwork(nloc),mplus(nloc)
@@ -2015,7 +2183,14 @@
           wt(il,i)=0.001
           water(il,i)=0.0
+          frac(il,i)=0.0
+          faci(il,i)=0.0
+          fraci(il,i)=0.0
+          ice(il,i)=0.0
+          prec(il,i)=0.0
+          fondue(il,i)=0.0
           evap(il,i)=0.0
           b(il,i)=0.0
           lvcp(il,i)=lv(il,i)/cpn(il,i)
+          lfcp(il,i)=lf(il,i)/cpn(il,i)
          enddo
         enddo
@@ -2115,9 +2290,28 @@
       wt(il,i)=45.0
 
+      if (cvflag_ice) then
+      frac(il,inb(il)) = 1. -(t(il,inb(il))-243.15)/(263.15-243.15)
+      frac(il,inb(il)) = min(max(frac(il,inb(il)),0.),1.)
+      fraci(il,inb(il)) = frac(il,inb(il))
+      else
+      continue
+      endif
+
       if(i.lt.inb(il))then
+
+      if (cvflag_ice) then
+      thaw = (t(il,i)-273.15)/(275.15-273.15)
+      thaw = min(max(thaw,0.0),1.0)
+      frac(il,i)=frac(il,i)*(1.-thaw)
+      else
+      continue
+      endif
+
        rp(il,i)=rp(il,i+1)
      :       +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
        rp(il,i)=0.5*(rp(il,i)+rr(il,i))
       endif
+      fraci(il,i)=1.-(t(il,i)-243.15)/(263.15-243.15)
+      fraci(il,i)=min(max(fraci(il,i),0.0),1.0)
       rp(il,i)=max(rp(il,i),0.0)
       rp(il,i)=amin1(rp(il,i),rs(il,i))
@@ -2126,4 +2320,7 @@
       if(i.eq.1)then
        afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))
+       if (cvflag_ice) then
+       afac1=p(il,i)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))
+       endif
       else
        rp(il,i-1)=rp(il,i)
@@ -2161,7 +2358,9 @@
 c       b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
 c       c6 = water(il,i+1) + wdtrain(il)*bfac
+c       c6 = prec(il,i+1) + wdtrain(il)*bfac
 c        revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
 c        evap(il,i)=sigt*afac*revap
 c        water(il,i)=revap*revap
+c        prec(il,i)=revap*revap
 cc        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ',
 cc     $            i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
@@ -2169,4 +2368,87 @@
 c
 c--------retour à la formulation originale d'Emanuel.
+       if (cvflag_ice) then
+
+c      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
+c      c6=prec(il,i+1)+bfac*wdtrain(il)
+c     :    -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
+c      if(c6.gt.0.0)then
+c       revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
+c
+cJAM  Attention: evap=sigt*E
+c     Modification: evap devient l'évaporation en milieu de couche
+c     car nécessaire dans cv3_yield
+c     Du coup, il faut modifier pas mal d'équations...
+c     et l'expression de afac qui devient afac1
+c     revap=sqrt((prec(i+1)+prec(i))/2)
+
+      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
+      c6=prec(il,i+1)+0.5*bfac*wdtrain(il)
+c      print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
+c      print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
+c      print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
+      if (c6 .gt. b6*b6+1.e-20) then
+       revap=2.*c6/(b6+sqrt(b6*b6+4.*c6))
+      else
+       revap=(-b6+sqrt(b6*b6+4.*c6))/2.
+      endif
+      prec(il,i)=Max(0.,2.*revap*revap-prec(il,i+1))
+c      print*,prec(il,i),'neige'
+
+cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
+cc             evap(il,i)=sigt*afac*revap
+c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
+c     Ici,l'evaporation evap est simplement calculee par l'equation de
+c     conservation.
+c       prec(il,i)=revap*revap
+c      else
+cjyg----   Correction : si c6 <= 0, water(il,i)=0.
+c       prec(il,i)=0.
+c      endif
+
+cjyg---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
+c                                    moins [tt ce qui sort de la couche i]
+c       print *, 'evap avec ice'
+       evap(il,i)=
+     :       (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i)))
+     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
+
+      d6=bfac*wdtrain(il)-100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))
+     :   *evap(il,i)
+      e6=bfac*wdtrain(il)
+      f6=-100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*
+     :   evap(il,i)
+
+      thaw=(t(il,i)-273.15)/(275.15-273.15)
+      thaw=min(max(thaw,0.0),1.0)
+      water(il,i)=water(il,i+1)+(1-fraci(il,i))*d6
+      water(il,i)=max(water(il,i),0.)
+      ice(il,i)=ice(il,i+1)+fraci(il,i)*d6
+      ice(il,i)=max(ice(il,i),0.)
+      fondue(il,i)=ice(il,i)*thaw
+      water(il,i)=water(il,i)+fondue(il,i)
+      ice(il,i)=ice(il,i)-fondue(il,i)
+
+      if(water(il,i)+ice(il,i).lt.1.e-30)then
+      faci(il,i)=0.
+      else
+      faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
+      endif
+
+c      water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
+c      water(il,i)=max(water(il,i),0.)
+c      ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
+c      ice(il,i)=max(ice(il,i),0.)
+c      fondue(il,i)=ice(il,i)*thaw
+c      water(il,i)=water(il,i)+fondue(il,i)
+c      ice(il,i)=ice(il,i)-fondue(il,i)
+c
+c      if((water(il,i)+ice(il,i)).lt.1.e-30)then
+c       faci(il,i)=0.
+c      else
+c       faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
+c      endif
+
+      else
       b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
       c6=water(il,i+1)+bfac*wdtrain(il)
@@ -2174,32 +2456,14 @@
       if(c6.gt.0.0)then
        revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
-cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
-cc             evap(il,i)=sigt*afac*revap
-c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
-c     Ici,l'evaporation evap est simplement calculee par l'equation de
-c     conservation.
        water(il,i)=revap*revap
       else
-cjyg----   Correction : si c6 <= 0, water(il,i)=0.
-       water(il,i) = 0.
+       water(il,i)= 0.
       endif
-cJYG/IM : ci-dessous formulation originale de KE
-c      evap(il,i)=-evap(il,i+1)
-c    :            +(wdtrain(il)+sigd(il)*wt(il,i)*water(il,i+1))
-c    :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*50.)
-c
-cJYG/IM : ci-dessous modification formulation originale de KE
-c        pour eliminer oscillations verticales de pluie se produisant
-c        lorsqu'il y a evaporation totale de la pluie
-c
-c       evap(il,i)= +(wdtrain(il)+sigd(il)*wt(il,i)*water(il,i+1)) !itlmd(jyg)
-c     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
-c      end if  !itlmd(jyg)
-cjyg---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
-c                                    moins [tt ce qui sort de la couche i]
+c       print *, 'evap sans ice'
        evap(il,i)=
      :       (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))
      :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
-c
+
+       endif
        endif !(i.le.inb(il) .and. lwork(il))
 995   Continue
@@ -2215,4 +2479,22 @@
       tevap(il)=max(0.0,evap(il,i))
       delth=max(0.001,(th(il,i)-th(il,i-1)))
+      if (cvflag_ice) then
+      if (cvflag_grav) then
+       mp(il,i)=100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)
+     :          *(p(il,i-1)-p(il,i))/delth
+     :          +lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)
+     :          *(p(il,i-1)-p(il,i))/delth
+     :          +lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)
+     :          *(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
+      else
+       mp(il,i)=10.*(lvcp(il,i)*sigd(il)*tevap(il)
+     :          *(p(il,i-1)-p(il,i))/delth
+     :          +lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)
+     :          *(p(il,i-1)-p(il,i))/delth
+     :          +lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)
+     :          *(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
+
+      endif
+      else
       if (cvflag_grav) then
        mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il)
@@ -2223,4 +2505,6 @@
       endif
 c
+      endif
+c
        endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
 996   Continue
@@ -2243,6 +2527,17 @@
      :               /(lvcp(il,i)*sigd(il)*th(il,i))
        af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
+
+       if (cvflag_ice) then
+       bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
+     :            +50.*(p(il,i-1)-p(il,i))*xf*
+     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
+     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/
+     :            (ph(il,i)-ph(il,i+1)))
+       else
+
        bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
      :            +50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
+       endif
+
        fac2=1.0
        if(bf.lt.0.0)fac2=-1.0
@@ -2262,9 +2557,28 @@
        mp(il,i)=max(0.0,mp(il,i))
 
+       if (cvflag_ice) then
        if (cvflag_grav) then
 Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
 C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
 C Et il faut bien revoir les facteurs 100.
-        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
+         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*
+     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
+     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)
+     :            /(ph(il,i)-ph(il,i+1)))
+     2            /(mp(il,i)+sigd(il)*0.1)
+     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
+     : *sigd(il)*th(il,i))
+       else
+        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*
+     :            (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i))+
+     :            (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)
+     :            /(ph(il,i)-ph(il,i+1)))
+     2            /(mp(il,i)+sigd(il)*0.1)
+     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
+     : *sigd(il)*th(il,i))
+       endif
+       else
+       if (cvflag_grav) then
+         b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
      2   /(mp(il,i)+sigd(il)*0.1)
      3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
@@ -2276,4 +2590,5 @@
      : *sigd(il)*th(il,i))
        endif
+      endif
        b(il,i-1)=max(b(il,i-1),0.0)
 c
@@ -2391,7 +2706,7 @@
      :                    ,icb,inb,delt
      :                    ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
-     :                    ,gz,p,ph,h,hp,lv,cpn,th,th_wake
+     :                    ,gz,p,ph,h,hp,lv,lf,cpn,th,th_wake
      :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
-     :                    ,wt,water,evap,b,sigd
+     :                    ,wt,water,ice,evap,fondue,faci,b,sigd
      :                    ,ment,qent,hent,iflag_mix,uent,vent
      :                    ,nent,elij,traent,sig
@@ -2422,7 +2737,9 @@
       real th(nloc,na), p(nloc,nd), tp(nloc,na)
       real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
+      real lf(nloc,na)
       real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
       real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
       real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc)
+      real fondue(nloc,na),faci(nloc,na), ice(nloc,na)
       real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
       real hent(nloc,na,na)
@@ -2455,5 +2772,5 @@
       real cpinv, rdcp, dpinv
       real awat(nloc)
-      real lvcp(nloc,na), mke(nloc,na)
+      real lvcp(nloc,na), lfcp(nloc,na), mke(nloc,na)
       real am(nloc), work(nloc), ad(nloc), amp1(nloc)
 c!!      real up1(nloc), dn1(nloc)
@@ -2513,4 +2830,5 @@
        do il=1,ncum
          lvcp(il,i)=lv(il,i)/cpn(il,i)
+         lfcp(il,i)=lf(il,i)/cpn(il,i)
        enddo
       enddo
@@ -2522,4 +2840,12 @@
       do il=1,ncum
        if(ep(il,inb(il)).ge.0.0001 .and. iflag(il) .le. 1)then
+        if (cvflag_ice) then
+        if (cvflag_grav) then
+           precip(il)=wt(il,1)*sigd(il)*(water(il,1)+ice(il,1))*86400.
+     :               *1000./(rowl*grav)
+        else
+         precip(il)=wt(il,1)*sigd(il)*(water(il,1)+ice(il,1))*8640.
+        endif
+        else
         if (cvflag_grav) then
            precip(il)=wt(il,1)*sigd(il)*water(il,1)*86400.*1000.
@@ -2527,4 +2853,5 @@
         else
          precip(il)=wt(il,1)*sigd(il)*water(il,1)*8640.
+        endif
         endif
        endif
@@ -2539,8 +2866,18 @@
        if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il)
      :    .and. iflag(il) .le. 1)then
+        if (cvflag_ice) then
+        if (cvflag_grav) then
+           VPrecip(il,i) = wt(il,i)*sigd(il)*(water(il,i)+ice(il,i))
+     :                     /grav
+        else
+           VPrecip(il,i) = wt(il,i)*sigd(il)*(water(il,i)+ice(il,i))
+     :                     /10.
+        endif
+        else
         if (cvflag_grav) then
            VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/grav
         else
            VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/10.
+        endif
         endif
        endif
@@ -2585,5 +2922,12 @@
 cjyg  Correction pour conserver l'eau
 ccc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2)) !precip
-       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)                  !precip
+       if (cvflag_ice) then
+       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)
+     :           -lfcp(il,1)*sigd(il)*evap(il,1)*faci(il,1)
+     :           -lfcp(il,1)*sigd(il)*(fondue(il,1)*wt(il,1))/
+     :           (100.*(ph(il,1)-ph(il,2)))                               !precip
+       else
+       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)
+       endif
 
       if (cvflag_grav) then
@@ -2595,6 +2939,13 @@
       endif
 
+      if (cvflag_ice) then
+        ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2)
+     :     *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)+
+     :     0.01*sigd(il)*wt(il,1)*(ci-cpd)*ice(il,2)*
+     :     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)
+      else
       ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2)
      :     *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)
+      endif
 
       ftd(il,1) = ft(il,1)                        ! fin precip
@@ -2791,5 +3142,13 @@
        ! precip
 ccc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
+       if (cvflag_ice) then
        ft(il,i)= -sigd(il)*lvcp(il,i)*evap(il,i)
+     :           -sigd(il)*lfcp(il,i)*evap(il,i)*faci(il,i)
+     :           -sigd(il)*lfcp(il,i)*fondue(il,i)*wt(il,i)/
+     :           (100.*(p(il,i-1)-p(il,i)))
+       else
+       ft(il,i)= -sigd(il)*lvcp(il,i)*evap(il,i)
+       endif
+
         rat=cpn(il,i-1)*cpinv
 c
@@ -2798,6 +3157,14 @@
      :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
      :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
+       if (cvflag_ice) then
        ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
      :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+     :           +0.01*sigd(il)*wt(il,i)*(ci-cpd)*ice(il,i+1)*
+     :           (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+      else
+       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
+     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+      endif
+
        ftd(il,i)=ft(il,i)
         ! fin precip
@@ -2818,6 +3185,15 @@
      :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
      :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
+
+       if (cvflag_ice) then
        ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
      :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+     :           +0.01*sigd(il)*wt(il,i)*(ci-cpd)*ice(il,i+1)*
+     :           (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+       else
+       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
+     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+       endif
+
        ftd(il,i)=ft(il,i)
         ! fin precip
@@ -3744,3 +4120,2 @@
         return
         end
-
Index: LMDZ5/trunk/libf/phylmd/cv3a_compress.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/cv3a_compress.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cv3a_compress.F	(revision 1849)
@@ -6,6 +6,6 @@
      :    ,u1,v1,gz1,th1,th1_wake
      :    ,tra1
-     :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
-     :    ,h1_wake,lv1_wake,cpn1_wake     ,tv1_wake
+     :    ,h1     ,lv1, lf1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
+     :    ,h1_wake,lv1_wake,lf1_wake,cpn1_wake     ,tv1_wake
      :    ,sig1,w01,ptop21
      :    ,Ale1,Alp1
@@ -16,6 +16,6 @@
      o    ,u,v,gz,th,th_wake
      o    ,tra
-     o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
-     o    ,h_wake,lv_wake,cpn_wake    ,tv_wake
+     o    ,h     ,lv, lf     ,cpn    ,p,ph,tv    ,tp,tvp,clw
+     o    ,h_wake,lv_wake,lf_wake,cpn_wake    ,tv_wake
      o    ,sig,w0,ptop2
      o    ,Ale,Alp  )
@@ -45,9 +45,9 @@
       real gz1(len,nd),th1(len,nd),th1_wake(len,nd)
       real tra1(len,nd,ntra)
-      real h1(len,nd),lv1(len,nd),cpn1(len,nd)
+      real h1(len,nd),lv1(len,nd),lf1(len,nd),cpn1(len,nd)
       real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
       real tvp1(len,nd),clw1(len,nd)
       real h1_wake(len,nd),lv1_wake(len,nd),cpn1_wake(len,nd)
-      real tv1_wake(len,nd)
+      real tv1_wake(len,nd),lf1_wake(len,nd)
       real sig1(len,nd), w01(len,nd), ptop21(len)
       real Ale1(len),Alp1(len)
@@ -65,9 +65,9 @@
       real gz(len,nd),th(len,nd),th_wake(len,nd)
       real tra(len,nd,ntra)
-      real h(len,nd),lv(len,nd),cpn(len,nd)
+      real h(len,nd),lv(len,nd),lf(len,nd),cpn(len,nd)
       real p(len,nd),ph(len,nd+1),tv(len,nd),tp(len,nd)
       real tvp(len,nd),clw(len,nd)
       real h_wake(len,nd),lv_wake(len,nd),cpn_wake(len,nd)
-      real tv_wake(len,nd)
+      real tv_wake(len,nd),lf_wake(len,nd)
       real sig(len,nd), w0(len,nd), ptop2(len)
       real Ale(len),Alp(len)
@@ -99,4 +99,5 @@
         h(nn,k)=h1(i,k)
         lv(nn,k)=lv1(i,k)
+        lf(nn,k)=lf1(i,k)
         cpn(nn,k)=cpn1(i,k)
         p(nn,k)=p1(i,k)
@@ -108,4 +109,5 @@
         h_wake(nn,k)=h1_wake(i,k)
         lv_wake(nn,k)=lv1_wake(i,k)
+        lf_wake(nn,k)=lf1_wake(i,k)
         cpn_wake(nn,k)=cpn1_wake(i,k)
         tv_wake(nn,k)=tv1_wake(i,k)
Index: LMDZ5/trunk/libf/phylmd/cv_driver.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/cv_driver.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cv_driver.F	(revision 1849)
@@ -345,5 +345,5 @@
 c   (common cvflag)
 
-       CALL cv_flag
+       CALL cv_flag(0)
 
 c -- set thermodynamical constants:
@@ -701,6 +701,10 @@
 
 !==================================================================
-      SUBROUTINE cv_flag
+      SUBROUTINE cv_flag(iflag_ice_thermo)
       implicit none
+
+c Argument : iflag_ice_thermo : ice thermodynamics is taken into account if
+c                               iflag_ice_thermo >=1
+      INTEGER iflag_ice_thermo
 
 #include "cvflag.h"
@@ -709,4 +713,5 @@
 c differente de 10.0 dans convect3: 
       cvflag_grav = .TRUE.
+      cvflag_ice = iflag_ice_thermo .GE. 1
 
       return
@@ -744,7 +749,9 @@
        cpv = RCPV
        cl  = RCW
+       ci  = RCS
        rrv = RV
        rrd = RD
        lv0 = RLVTT
+       lf0 = RLSTT-RLVTT
        g   = RG     ! not used in convect3
 c ori      t0  = RTT
@@ -758,4 +765,5 @@
       clmcpv=cl-cpv
       clmcpd=cl-cpd
+      clmci=cl-ci
       cpdmcp=cpd-cpv
       cpvmcpd=cpv-cpd
Index: LMDZ5/trunk/libf/phylmd/cva_driver.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/cva_driver.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cva_driver.F	(revision 1849)
@@ -3,5 +3,5 @@
 !
       SUBROUTINE cva_driver(len,nd,ndp1,ntra,nloc,
-     &                   iflag_con,iflag_mix,
+     &                   iflag_con,iflag_mix,iflag_ice_thermo,
      &                   iflag_clos,delt,
      &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake,
@@ -52,4 +52,5 @@
 C      iflag_con     Integer        Input        version of convect (3/4)
 C      iflag_mix     Integer        Input        version of mixing  (0/1/2)
+C   iflag_ice_thermo Integer        Input        accounting for ice thermodynamics (0/1)
 C      iflag_clos    Integer        Input        version of closure (0/1)
 C      delt          Real           Input        time step
@@ -155,4 +156,5 @@
       integer iflag_con
       integer iflag_mix
+      integer iflag_ice_thermo
       integer iflag_clos
       real delt
@@ -381,4 +383,5 @@
       real buoybase1(klon)
 
+      real lf1(klon,klev) ,lf1_wake(klon,klev)
       real lv1(klon,klev) ,lv1_wake(klon,klev)
       real cpn1(klon,klev),cpn1_wake(klon,klev)
@@ -419,6 +422,6 @@
       real gz(nloc,klev),h(nloc,klev)     ,hm(nloc,klev)
       real               h_wake(nloc,klev),hm_wake(nloc,klev)
-      real lv(nloc,klev)     ,cpn(nloc,klev)
-      real lv_wake(nloc,klev),cpn_wake(nloc,klev)
+      real lv(nloc,klev), lf(nloc,klev), cpn(nloc,klev)
+      real lv_wake(nloc,klev), lf_wake(nloc,klev), cpn_wake(nloc,klev)
       real p(nloc,klev),ph(nloc,klev+1),tv(nloc,klev)    ,tp(nloc,klev)
       real                              tv_wake(nloc,klev)
@@ -430,5 +433,5 @@
       real sig(nloc,klev), w0(nloc,klev), ptop2(nloc)
       real hp(nloc,klev), ep(nloc,klev), sigp(nloc,klev)
-      real frac(nloc), buoy(nloc,klev)
+      real buoy(nloc,klev)
       real cape(nloc)
       real cin(nloc)
@@ -444,10 +447,12 @@
       real sigd(nloc)
 !      real mp(nloc,klev), qp(nloc,klev), up(nloc,klev), vp(nloc,klev)
-!      real wt(nloc,klev), water(nloc,klev), evap(nloc,klev)
+!      real wt(nloc,klev), water(nloc,klev), evap(nloc,klev), ice(nloc,klev)
 !      real b(nloc,klev), sigd(nloc)
 !      save mp,qp,up,vp,wt,water,evap,b
       real, save, allocatable :: mp(:,:),qp(:,:),up(:,:),vp(:,:)
-      real, save, allocatable :: wt(:,:),water(:,:),evap(:,:), b(:,:)
-c$OMP THREADPRIVATE(mp,qp,up,vp,wt,water,evap,b)
+      real, save, allocatable :: wt(:,:),water(:,:),evap(:,:)
+      real, save, allocatable :: ice(:,:),fondue(:,:), b(:,:)
+      real, save, allocatable :: frac(:,:), faci(:,:)
+c$OMP THREADPRIVATE(mp,qp,up,vp,wt,water,evap,ice,fondue,b,frac,faci)
       real  ft(nloc,klev), fq(nloc,klev)
       real ftd(nloc,klev), fqd(nloc,klev)
@@ -496,5 +501,7 @@
          allocate(mp(nloc,klev), qp(nloc,klev), up(nloc,klev))
          allocate(vp(nloc,klev), wt(nloc,klev), water(nloc,klev))
+         allocate(ice(nloc,klev), fondue(nloc,klev))
          allocate(evap(nloc,klev), b(nloc,klev))
+         allocate(frac(nloc,klev), faci(nloc,klev))
          first=.false.
        endif
@@ -502,5 +509,5 @@
 c   (common cvflag)
 
-       CALL cv_flag
+       CALL cv_flag(iflag_ice_thermo)
 
 c -- set thermodynamical constants:
@@ -610,9 +617,9 @@
 !       print*,'t1, q1 ',t1,q1
        CALL cv3_prelim(len,nd,ndp1,t1,q1,p1,ph1      ! nd->na
-     o               ,lv1,cpn1,tv1,gz1,h1,hm1,th1)
+     o               ,lv1,lf1,cpn1,tv1,gz1,h1,hm1,th1)
     
 c
        CALL cv3_prelim(len,nd,ndp1,t1_wake,q1_wake,p1,ph1 ! nd->na
-     o               ,lv1_wake,cpn1_wake,tv1_wake,gz1_wake
+     o               ,lv1_wake,lf1_wake,cpn1_wake,tv1_wake,gz1_wake
      o               ,h1_wake,bid,th1_wake)
     
@@ -755,6 +762,6 @@
      :    ,u1,v1,gz1,th1,th1_wake
      :    ,tra1
-     :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
-     :    ,h1_wake,lv1_wake,cpn1_wake     ,tv1_wake
+     :    ,h1     ,lv1, lf1,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
+     :    ,h1_wake,lv1_wake,lf1_wake,cpn1_wake     ,tv1_wake
      :    ,sig1,w01,ptop21
      :    ,Ale1,Alp1
@@ -765,6 +772,6 @@
      o    ,u,v,gz,th,th_wake
      o    ,tra
-     o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
-     o    ,h_wake,lv_wake,cpn_wake    ,tv_wake
+     o    ,h     ,lv, lf     ,cpn    ,p,ph,tv    ,tp,tvp,clw
+     o    ,h_wake,lv_wake,lf_wake,cpn_wake    ,tv_wake
      o    ,sig,w0,ptop2
      o    ,Ale,Alp  ) 
@@ -800,6 +807,7 @@
        CALL cv3_undilute2(nloc,ncum,nd,icb,icbs,nk        !na->nd
      :                        ,tnk,qnk,gznk,hnk,t,q,qs,gz
-     :                        ,p,h,tv,lv,pbase,buoybase,plcl
-     o                        ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
+     :                        ,p,h,tv,lv,lf,pbase,buoybase,plcl
+     o                        ,inb,tp,tvp,clw,hp,ep,sigp,buoy
+     :                        ,frac)
    
       endif
@@ -816,4 +824,9 @@
 !-------------------------------------------------------------------
       IF (iflag_con .eq. 3) THEN
+       if ((iflag_ice_thermo.eq.1).and.(iflag_mix.ne.0)) then
+         write(*,*) " iflag_ice_thermo==1 requires iflag_mix==0",
+     &   " but iflag_mix=",iflag_mix, ". Might as well stop here."
+         stop
+       endif
        IF (iflag_mix .ge. 1 ) THEN
          CALL zilch(supmax,nloc*klev)    
@@ -877,5 +890,5 @@
         IF (iflag_mix.eq.0) THEN
          CALL cv3_mixing(nloc,ncum,nd,nd,ntra,icb,nk,inb    ! na->nd
-     :                       ,ph,t,q,qs,u,v,tra,h,lv,qnk
+     :                       ,ph,t,q,qs,u,v,tra,h,lv,lf,frac,qnk
      :                       ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
      o   ,ment,qent,uent,vent,nent,sigij,elij,ments,qents,traent)
@@ -913,8 +926,9 @@
        CALL cv3_unsat(nloc,ncum,nd,nd,ntra,icb,inb,iflag    ! na->nd
      :               ,t_wake,q_wake,qs_wake,gz,u,v,tra,p,ph
-     :               ,th_wake,tv_wake,lv_wake,cpn_wake
+     :               ,th_wake,tv_wake,lv_wake,lf_wake,cpn_wake
      :               ,ep,sigp,clw
      :               ,m,ment,elij,delt,plcl,coef_clos
-     o          ,mp,qp,up,vp,trap,wt,water,evap,b,sigd
+     o          ,mp,qp,up,vp,trap,wt,water,evap,fondue,ice
+     :          ,faci,b,sigd
      o          ,wdtrainA,wdtrainM)   ! RomP
       endif
@@ -943,7 +957,7 @@
      :                     ,icb,inb,delt
      :                     ,t,q,t_wake,q_wake,s_wake,u,v,tra
-     :                     ,gz,p,ph,h,hp,lv,cpn,th,th_wake
+     :                     ,gz,p,ph,h,hp,lv,lf,cpn,th,th_wake
      :                     ,ep,clw,m,tp,mp,qp,up,vp,trap
-     :                     ,wt,water,evap,b,sigd
+     :                     ,wt,water,ice,evap,fondue,faci,b,sigd
      :                    ,ment,qent,hent,iflag_mix,uent,vent
      :                    ,nent,elij,traent,sig
@@ -1043,3 +1057,2 @@
       return
       end
-
Index: LMDZ5/trunk/libf/phylmd/cvflag.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/cvflag.h	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cvflag.h	(revision 1849)
@@ -3,5 +3,6 @@
 !
       logical cvflag_grav
+      logical cvflag_ice
 
-      COMMON /cvflag/ cvflag_grav 
+      COMMON /cvflag/ cvflag_grav, cvflag_ice  
 c$OMP THREADPRIVATE(/cvflag/)
Index: LMDZ5/trunk/libf/phylmd/cvthermo.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/cvthermo.h	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/cvthermo.h	(revision 1849)
@@ -4,13 +4,13 @@
 c Thermodynamical constants for convectL:
 
-      real cpd, cpv, cl, rrv, rrd, lv0, g, rowl, t0
-      real clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl
+      real cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl, t0
+      real clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl, clmci
       real eps, epsi, epsim1
       real ginv, hrd
       real grav
 
-      COMMON /cvthermo/ cpd, cpv, cl, rrv, rrd, lv0, g, rowl, t0
-     :                 ,clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl 
-     :                 ,eps, epsi, epsim1, ginv, hrd, grav
+      COMMON /cvthermo/ cpd, cpv, cl, ci, rrv, rrd, lv0, lf0, g, rowl
+     :                 ,t0, clmcpv, clmcpd, cpdmcp, cpvmcpd, cpvmcl 
+     :                 ,clmci, eps, epsi, epsim1, ginv, hrd, grav
 
 c$OMP THREADPRIVATE(/cvthermo/)
Index: LMDZ5/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 1849)
@@ -8,5 +8,6 @@
      frac_impa, frac_nucl, beta,                        &
      prfl, psfl, rhcl, zqta, fraca,                     &
-     ztv, zpspsk, ztla, zthl, iflag_cldcon)
+     ztv, zpspsk, ztla, zthl, iflag_cldcon,             &
+     iflag_ice_thermo)
 
   !
@@ -51,6 +52,8 @@
   REAL zpspsk(klon,klev),ztla(klon,klev)
   REAL zthl(klon,klev)
+  REAL ztfondue, qsl, qsi
 
   logical lognormale(klon)
+  logical ice_thermo
 
   !AA
@@ -77,4 +80,5 @@
   INTEGER ncoreczq  
   INTEGER iflag_cldcon
+  INTEGER iflag_ice_thermo
   PARAMETER (ninter=5)
   LOGICAL evap_prec ! evaporation de la pluie
@@ -97,10 +101,14 @@
   INTEGER i, k, n, kk
   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
-  REAL zrfl(klon), zrfln(klon), zqev, zqevt
-  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
+  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 ztglace, zt(klon)
   INTEGER nexpo ! exponentiel pour glace/eau
   REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
   REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
+  REAL zmelt, zpluie, zice, zcondold
+  PARAMETER (ztfondue=278.15)
   !
   LOGICAL appel1er
@@ -145,4 +153,5 @@
   DATA appel1er /.TRUE./
   !ym
+  ice_thermo = iflag_ice_thermo .GE. 1
   zdelq=0.0
 
@@ -188,5 +197,13 @@
   !
   ztglace = RTT - 15.0
-  nexpo = 6
+!AJ<
+  IF (ice_thermo) THEN
+    nexpo = 2
+  ELSE
+    nexpo = 6
+  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
   !
@@ -223,4 +240,5 @@
      !     DO i = 1, klon
      zrfl(i) = 0.0
+     zifl(i) = 0.0
      zneb(i) = seuil_neb
   ENDDO
@@ -272,7 +290,14 @@
 
 
+     ! Calculer l'evaporation de la precipitation
+     !
+
+
      IF (evap_prec) THEN
         DO i = 1, klon
-           IF (zrfl(i) .GT.0.) THEN
+!AJ<
+!!           IF (zrfl(i) .GT.0.) THEN
+           IF (zrfl(i)+zifl(i).GT.0.) THEN
+!>AJ
               IF (thermcep) THEN
                  zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
@@ -288,29 +313,129 @@
                  ENDIF
               ENDIF
-              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
-              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))
-              zqev = MIN (zqev, zqevt)
-              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
-                   /RG/dtime
-
-              ! pour la glace, on ré-évapore toute la précip dans la
-              ! couche du dessous
-              ! la glace venant de la couche du dessus est simplement
-              ! dans la couche du dessous.
-
-              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-
-              zq(i) = zq(i) - (zrfln(i)-zrfl(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))
-              zrfl(i) = zrfln(i)
-           ENDIF
-        ENDDO
-     ENDIF
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+        ENDDO
+!AJ<
+       IF (.NOT. ice_thermo) THEN
+        DO i = 1, klon
+!AJ<
+!!           IF (zrfl(i) .GT.0.) THEN
+           IF (zrfl(i)+zifl(i).GT.0.) THEN
+!>AJ
+                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+                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))
+                zqev = MIN (zqev, zqevt)
+                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
+                     /RG/dtime
+       
+                ! pour la glace, on ré-évapore toute la précip dans la
+                ! couche du dessous
+                ! la glace venant de la couche du dessus est simplement
+                ! dans la couche du dessous.
+       
+                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
+       
+                zq(i) = zq(i) - (zrfln(i)-zrfl(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))
+                zrfl(i) = zrfln(i)
+                zifl(i) = 0.
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+        ENDDO
+!
+       ELSE ! (.NOT. ice_thermo) 
+!
+        DO i = 1, klon
+!AJ<
+!!           IF (zrfl(i) .GT.0.) THEN
+           IF (zrfl(i)+zifl(i).GT.0.) THEN
+!>AJ
+     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! Modification de l'évaporation avec la glace
+     ! Différentiation entre précipitation liquide et solide 
+     ! On suppose que coef_evai=2*coef_eva
+     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
+     
+         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+     !    zqev0 = MAX (0.0, zqs(i)-zq(i) ) 
+
+     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! On différencie qsat pour l'eau et la glace
+     ! Si zdelta=1. --> glace
+     ! Si zdelta=0. --> eau liquide
+     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+         
+         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
+         qsl= MIN(0.5,qsl)
+         zcor= 1./(1.-RETV*qsl)
+         qsl= qsl*zcor
+         
+         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
+         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+              *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
+
+         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
+         qsi= MIN(0.5,qsi)
+         zcor= 1./(1.-RETV*qsi)
+         qsi= qsi*zcor
+
+         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
+         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
+              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
+
+              
+     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! Vérification sur l'évaporation
+     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     
+         IF (zqevt+zqevti.GT.zqev0) THEN
+	          zqev=zqev0*zqevt/(zqevt+zqevti)
+	          zqevi=zqev0*zqevti/(zqevt+zqevti)
+             
+         ELSE
+             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
+     
+         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)
+         
+     ! Pour la glace, on révapore toute la précip dans la couche du dessous
+     ! la glace venant de la couche du dessus est simplement dans la couche
+     ! du dessous.
+     
+     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
+!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
+         zq(i) = zq(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)) &
+                  + (zifln(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
+     
+         zrfl(i) = zrfln(i)
+         zifl(i) = zifln(i)
+
+           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+        ENDDO
+
+      ENDIF ! (.NOT. ice_thermo) 
+     
+     ENDIF ! (evap_prec)
      !
      ! Calculer Qs et L/Cp*dQs/dT:
@@ -478,6 +603,21 @@
         zq(i) = zq(i) - zcond(i)
         !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
-        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-     ENDDO
+     ENDDO
+!AJ<
+     IF (.NOT. ice_thermo) THEN
+       DO i = 1, klon
+         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+       ENDDO
+     ELSE
+       DO i = 1, klon
+         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
+         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+         zfice(i) = zfice(i)**nexpo
+         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))
+!         print*,zt(i),zrfl(i),zifl(i),'temp1'
+       ENDDO
+     ENDIF
+!>AJ
      !
      ! Partager l'eau condensee en precipitation et eau liquide nuageuse
@@ -488,8 +628,23 @@
            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
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
            zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
            zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
            zfice(i) = zfice(i)**nexpo
+     !!      zfice(i)=0.
+        ENDIF
+     ENDDO
+     ENDIF
+     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
@@ -502,4 +657,6 @@
 
               IF (zneb(i).EQ.seuil_neb) THEN
+                 zpluie= 0.0
+                 zice = 0.0  
                  ztot = 0.0
               ELSE
@@ -519,9 +676,23 @@
                  zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
                       *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
-                 ztot    = zchau    + zfroi
+!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)
+              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
@@ -529,16 +700,67 @@
      ENDDO
      !
-     DO i = 1, klon
-        IF (rneb(i,k).GT.0.0) THEN
+         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)
-        ENDIF
-        IF (zt(i).LT.RTT) THEN
+         ENDIF
+       ENDDO
+     ELSE
+       DO i = 1, klon
+         IF (rneb(i,k).GT.0.0) THEN
+           d_ql(i,k) = zoliq(i)
+!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)                                    
+
+         ENDIF                      
+       ENDDO
+     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)
+!           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)
-        ELSE
+         ELSE
            prfl(i,k)=zrfl(i)
-        ENDIF
-     ENDDO
+         ENDIF
+       ENDDO
+     ELSE
+     ! JAM*************************************************
+     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
+     ! *****************************************************
+
+       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
+     !
      !
      ! Calculer les tendances de q et de t:
@@ -604,5 +826,8 @@
   DO i = 1, klon
      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
-        snow(i) = zrfl(i)
+!AJ<
+!!        snow(i) = zrfl(i)
+        snow(i) = zrfl(i)+zifl(i)
+!>AJ
         zlh_solid(i) = RLSTT-RLVTT
      ELSE
Index: LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1848)
+++ LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1849)
@@ -1878,6 +1878,13 @@
       DO i = 1, klon
          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
+cjyg<
+c      Attention : Arnaud a propose des formules completement differentes
+c                  A verifier !!!
+         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
+         IF (iflag_ice_thermo .EQ. 0) THEN
+           zlsdcp=zlvdcp
+         ENDIF 
+c>jyg
+
          zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
          zb = MAX(0.0,ql_seri(i,k))
@@ -2263,5 +2270,7 @@
                 nbtr_tmp=nbtr
              END IF
-          CALL concvl (iflag_con,iflag_clos,
+cjyg   iflag_con est dans clesphys
+cc          CALL concvl (iflag_con,iflag_clos,
+          CALL concvl (iflag_clos,
      .        dtime,paprs,pplay,t_undi,q_undi,
      .        t_wake,q_wake,wake_s,
@@ -2765,5 +2774,6 @@
      .           frac_impa, frac_nucl, beta_prec_fisrt,
      .           prfl, psfl, rhcl, 
-     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
+     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon,
+     .           iflag_ice_thermo)
 
       WHERE (rain_lsc < 0) rain_lsc = 0.
