Index: LMDZ5/trunk/libf/phylmd/cloudth_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cloudth_mod.F90	(revision 2911)
+++ 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))
       
