Index: /LMDZ5/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/field_def_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/field_def_lmdz.xml	(revision 2945)
@@ -515,4 +515,5 @@
         <field id="rnebcon"    long_name="Convective Cloud Fraction"    unit="-" />
         <field id="rnebls"    long_name="LS Cloud fraction"    unit="-" />
+        <field id="rneblsvol" long_name="LS Cloud fraction by volume"    unit="-" />
         <field id="rhum"    long_name="Relative humidity"    unit="-" />
         <field id="ozone"    long_name="Ozone mole fraction"    unit="-" />
Index: /LMDZ5/trunk/DefLists/file_def_histLES_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histLES_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histLES_lmdz.xml	(revision 2945)
@@ -474,4 +474,5 @@
                 <field field_ref="rnebcon" level="10" />
                 <field field_ref="rnebls" level="10" />
+                <field field_ref="rneblsvol" level="10" />
                 <field field_ref="rhum" level="10" />
                 <field field_ref="ozone" level="10" />
Index: /LMDZ5/trunk/DefLists/file_def_histday_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histday_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histday_lmdz.xml	(revision 2945)
@@ -481,4 +481,5 @@
                 <field field_ref="rnebcon" level="5" />
                 <field field_ref="rnebls" level="5" />
+                <field field_ref="rneblsvol" level="5" />
                 <field field_ref="rhum" level="5" />
                 <field field_ref="ozone" level="10" />
Index: /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histhf_lmdz.xml	(revision 2945)
@@ -504,4 +504,5 @@
                 <field field_ref="rnebcon" level="10" />
                 <field field_ref="rnebls" level="10" />
+                <field field_ref="rneblsvol" level="10" />
                 <field field_ref="rhum" level="10" />
                 <field field_ref="ozone" level="10" />
Index: /LMDZ5/trunk/DefLists/file_def_histins_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histins_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histins_lmdz.xml	(revision 2945)
@@ -474,4 +474,5 @@
                 <field field_ref="rnebcon" level="10" />
                 <field field_ref="rnebls" level="10" />
+                <field field_ref="rneblsvol" level="10" />
                 <field field_ref="rhum" level="10" />
                 <field field_ref="ozone" level="10" />
Index: /LMDZ5/trunk/DefLists/file_def_histmth_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histmth_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histmth_lmdz.xml	(revision 2945)
@@ -514,4 +514,5 @@
                 <field field_ref="rnebcon" level="2" />
                 <field field_ref="rnebls" level="2" />
+                <field field_ref="rneblsvol" level="2" />
                 <field field_ref="rhum" level="2" />
                 <field field_ref="ozone" level="2" />
Index: /LMDZ5/trunk/DefLists/file_def_histstn_lmdz.xml
===================================================================
--- /LMDZ5/trunk/DefLists/file_def_histstn_lmdz.xml	(revision 2944)
+++ /LMDZ5/trunk/DefLists/file_def_histstn_lmdz.xml	(revision 2945)
@@ -474,4 +474,5 @@
                 <field field_ref="rnebcon" level="10" />
                 <field field_ref="rnebls" level="10" />
+                <field field_ref="rneblsvol" level="10" />
                 <field field_ref="rhum" level="10" />
                 <field field_ref="ozone" level="10" />
Index: /LMDZ5/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/cloudth_mod.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/cloudth_mod.F90	(revision 2945)
@@ -587,5 +587,5 @@
        SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -624,14 +624,16 @@
       REAL zqsatenv(ngrid,klev)
       
-      
-      REAL sigma1(ngrid,klev)
+      REAL sigma1(ngrid,klev)                                                         
       REAL sigma2(ngrid,klev)
       REAL qlth(ngrid,klev)
       REAL qlenv(ngrid,klev)
       REAL qltot(ngrid,klev) 
-      REAL cth(ngrid,klev)  
+      REAL cth(ngrid,klev)
       REAL cenv(ngrid,klev)   
       REAL ctot(ngrid,klev)
-      REAL rneb(ngrid,klev)
+      REAL cth_vol(ngrid,klev)
+      REAL cenv_vol(ngrid,klev)
+      REAL ctot_vol(ngrid,klev)
+      REAL rneb(ngrid,klev)      
       REAL t(ngrid,klev)
       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
@@ -652,5 +654,5 @@
       CALL cloudth_vert_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
      &           ratqs,zqs,t)
       RETURN
@@ -672,4 +674,7 @@
       cenv(:,:)=0.
       ctot(:,:)=0.
+      cth_vol(:,:)=0.
+      cenv_vol(:,:)=0.
+      ctot_vol(:,:)=0.
       qsatmmussig1=0.
       qsatmmussig2=0.
@@ -747,4 +752,5 @@
       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
 
       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
@@ -780,4 +786,5 @@
       xenv=senv/(sqrt2*sigma1s)
       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
 
@@ -801,5 +808,5 @@
      SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
      &           ztv,po,zqta,fraca, & 
-     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+     &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
      &           ratqs,zqs,t)
 
@@ -838,5 +845,7 @@
       REAL zqsatenv(ngrid,klev)
       
-      
+   
+
+
       REAL sigma1(ngrid,klev)                                                         
       REAL sigma2(ngrid,klev)
@@ -844,7 +853,10 @@
       REAL qlenv(ngrid,klev)
       REAL qltot(ngrid,klev) 
-      REAL cth(ngrid,klev)  
+      REAL cth(ngrid,klev)
       REAL cenv(ngrid,klev)   
       REAL ctot(ngrid,klev)
+      REAL cth_vol(ngrid,klev)
+      REAL cenv_vol(ngrid,klev)
+      REAL ctot_vol(ngrid,klev)
       REAL rneb(ngrid,klev)
       REAL t(ngrid,klev)                                                                  
@@ -854,5 +866,5 @@
       REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
       REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
-      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
+      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
       REAL Tbef,zdelta,qsatbef,zcor
       REAL qlbef  
@@ -883,4 +895,7 @@
       cenv(:,:)=0.
       ctot(:,:)=0.
+      cth_vol(:,:)=0.
+      cenv_vol(:,:)=0.
+      ctot_vol(:,:)=0.
       qsatmmussig1=0.
       qsatmmussig2=0.
@@ -1032,39 +1047,53 @@
       exp_xth2 = exp(-1.*xth2**2)
       
+      !CF_surfacique
       cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
       cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
-      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
-
-      
+      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 
+
+
+      !CF_volumique & eau condense
       !environnement
       IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
+      IntJ_CF=0.5*(1.-1.*erf(xenv2))
       if (deltasenv .lt. 1.e-10) then 
       qlenv(ind1,ind2)=IntJ
+      cenv_vol(ind1,ind2)=IntJ_CF
       else
       IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
       IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
       IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
+      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
+      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
       qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
       endif
-
 
       !thermique
       IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
-      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
+      IntJ_CF=0.5*(1.-1.*erf(xth2))
+      if (deltasth .lt. 1.e-10) then
       qlth(ind1,ind2)=IntJ
+      cth_vol(ind1,ind2)=IntJ_CF
       else
       IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
       IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
       IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
+      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
+      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
       qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
+      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
       endif
 
+
       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
-
+      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
 
       ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
+
 
       if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
       ctot(ind1,ind2)=0.
+      ctot_vol(ind1,ind2)=0.
       qcloud(ind1)=zqsatenv(ind1,ind2) 
 
@@ -1100,4 +1129,5 @@
       xenv=senv/(sqrt2*sigma1s)
       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
+      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
       
Index: /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 2945)
@@ -169,4 +169,5 @@
     INTEGER,SAVE :: iflag_t_glace_omp
     INTEGER,SAVE :: iflag_cloudth_vert_omp
+    INTEGER,SAVE :: iflag_rain_incloud_vol_omp
     REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
     REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
@@ -1248,4 +1249,13 @@
 
     !
+    !Config Key  = iflag_rain_incloud_vol
+    !Config Desc =  
+    !Config Def  = 0
+    !Config Help = 
+    !
+    iflag_rain_incloud_vol_omp = 0
+    CALL getin('iflag_rain_incloud_vol',iflag_rain_incloud_vol_omp)
+
+    !
     !Config Key  = iflag_ice_thermo
     !Config Desc =  
@@ -2163,4 +2173,5 @@
     iflag_t_glace = iflag_t_glace_omp
     iflag_cloudth_vert=iflag_cloudth_vert_omp
+    iflag_rain_incloud_vol=iflag_rain_incloud_vol_omp
     iflag_ice_thermo = iflag_ice_thermo_omp
     rei_min = rei_min_omp
@@ -2511,4 +2522,5 @@
     write(lunout,*)' iflag_t_glace = ',iflag_t_glace
     write(lunout,*)' iflag_cloudth_vert = ',iflag_cloudth_vert
+    write(lunout,*)' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol
     write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo
     write(lunout,*)' rei_min = ',rei_min
Index: /LMDZ5/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2945)
@@ -18,4 +18,5 @@
   USE ioipsl_getin_p_mod, ONLY : getin_p
   USE phys_local_var_mod, ONLY: ql_seri,qs_seri
+  USE phys_local_var_mod, ONLY: rneblsvol
   ! flag to include modifications to ensure energy conservation (if flag >0)
   USE add_phys_tend_mod, only : fl_cor_ebil 
@@ -132,4 +133,5 @@
   INTEGER ncoreczq  
   REAL ctot(klon,klev)
+  REAL ctot_vol(klon,klev)
   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
   REAL zdqsdT_raw(klon)
@@ -217,9 +219,11 @@
 !  ice_thermo = iflag_ice_thermo .GE. 1
 
-itap=itap+1
-znebprecip(:)=0.
-
-   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
+  itap=itap+1
+  znebprecip(:)=0.
+
+  ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
   zdelq=0.0
+  ctot_vol(1:klon,1:klev)=0.0
+  rneblsvol(1:klon,1:klev)=0.0
 
   if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
@@ -727,9 +731,10 @@
                call cloudth_v3(klon,klev,k,ztv, &
                    zq,zqta,fraca, &
-                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+                   qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
                    ratqs,zqs,t)
               endif
               do i=1,klon
                  rneb(i,k)=ctot(i,k)
+                 rneblsvol(i,k)=ctot_vol(i,k)
                  zqn(i)=qcloud(i)
               enddo
@@ -1148,6 +1153,16 @@
                          *fallvs(zrhol(i)) * zfice(i)
                  endif
-                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
-                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
+
+                 ! si l'heterogeneite verticale est active, on utilise 
+                 ! la fraction volumique "vraie" plutot que la fraction
+                 ! surfacique modifiee, qui est plus grande et reduit
+                 ! sinon l'eau in-cloud de facon artificielle
+                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
+                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
+                 else
+                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
+                 endif
 !AJ<
                  IF (.NOT. ice_thermo) THEN
Index: /LMDZ5/trunk/libf/phylmd/nuage.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/nuage.h	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/nuage.h	(revision 2945)
@@ -10,4 +10,5 @@
 
       INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
+      INTEGER iflag_rain_incloud_vol
 
       common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
@@ -15,4 +16,5 @@
      &                  tau_cld_cv,coefw_cld_cv,                        &
      &                  tmax_fonte_cv,                                  &
-     &                  iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv
+     &                  iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv,  &
+     &                  iflag_rain_incloud_vol
 !$OMP THREADPRIVATE(/nuagecom/)
Index: /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2945)
@@ -417,6 +417,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: beta_prec
 !$OMP THREADPRIVATE(beta_prec)
-      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rneb,rnebjn
-!$OMP THREADPRIVATE(rneb,rnebjn)
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rneb,rnebjn,rneblsvol
+!$OMP THREADPRIVATE(rneb,rnebjn,rneblsvol)
 
 ! variables de sorties MM
@@ -732,5 +732,5 @@
       ALLOCATE(wdtrainA(klon,klev),wdtrainM(klon,klev))
       ALLOCATE(beta_prec(klon,klev))
-      ALLOCATE(rneb(klon,klev),rnebjn(klon,klev))
+      ALLOCATE(rneb(klon,klev),rnebjn(klon,klev),rneblsvol(klon,klev))
 
 
Index: /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2945)
@@ -1328,4 +1328,6 @@
   TYPE(ctrl_out), SAVE :: o_rnebls = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rnebls', 'LS Cloud fraction', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_rneblsvol = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'rneblsvol', 'LS Cloud fraction by volume', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
Index: /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2944)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2945)
@@ -124,5 +124,5 @@
          o_vitu, o_vitv, o_vitw, o_pres, o_paprs, &
          o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
-         o_rnebls, o_rhum, o_ozone, o_ozone_light, &
+         o_rnebls, o_rneblsvol, o_rhum, o_ozone, o_ozone_light, &
          o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
          o_dqsphy, o_dqsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
@@ -269,5 +269,5 @@
          ql_seri, qs_seri, tr_seri, &
          zphi, u_seri, v_seri, omega, cldfra, &
-         rneb, rnebjn, zx_rh, d_t_dyn,  & 
+         rneb, rnebjn, rneblsvol, zx_rh, d_t_dyn,  & 
          d_q_dyn,  d_ql_dyn, d_qs_dyn, &
          d_q_dyn2d,  d_ql_dyn2d, d_qs_dyn2d, &
@@ -1344,4 +1344,5 @@
        CALL histwrite_phy(o_rnebcon, rnebcon)
        CALL histwrite_phy(o_rnebls, rneb)
+       CALL histwrite_phy(o_rneblsvol, rneblsvol)
        IF (vars_defined)  THEN
           DO k=1, klev
