Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F	(revision 373)
@@ -1,8 +1,9 @@
 c $Header$
-      SUBROUTINE conema (dtime,paprs,pplay,t,q,u,v,tra,ntra,
+      SUBROUTINE conema3 (dtime,paprs,pplay,t,q,u,v,tra,ntra,
      .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
      .             rain, snow, kbas, ktop,
      .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
-     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
+     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
+     .             qcond_incld)
 
       IMPLICIT none
@@ -52,4 +53,5 @@
 #include "dimensions.h"
 #include "dimphy.h"
+#include "conema3.h"
       INTEGER i, l,m,itra
       INTEGER ntra,ntrac !number of tracers; if no tracer transport
@@ -58,5 +60,8 @@
       REAL dtime
 c
-
+      REAL d_t2(klon,klev), d_q2(klon,klev) ! sbl
+      REAL d_u2(klon,klev), d_v2(klon,klev) ! sbl
+      REAL em_d_t2(klev), em_d_q2(klev)     ! sbl   
+      REAL em_d_u2(klev), em_d_v2(klev)     ! sbl   
 c 
       REAL paprs(klon,klev+1), pplay(klon,klev)
@@ -73,4 +78,7 @@
       REAL dplcldt(klon), dplcldr(klon)
       INTEGER kbas(klon), ktop(klon)
+
+      REAL wd(klon)
+      REAL qcond_incld(klon,klev)
 c
       REAL em_t(klev)
@@ -91,4 +99,8 @@
       INTEGER em_bas, em_top
       SAVE em_bas, em_top
+
+      REAL em_wd
+      REAL em_qcond(klev)
+      REAL em_qcondc(klev)
 c
       REAL zx_t, zx_qs, zdelta, zcor
@@ -115,4 +127,21 @@
       real em_pbase, em_bbase
       integer iw,j,k,ix,iy
+
+c -- sb: pour schema nuages:
+
+       integer iflagcon
+       integer em_ifc(klev)
+     
+       real em_pradj
+       real em_cldf(klev), em_cldq(klev)
+       real em_ftadj(klev), em_fradj(klev)
+
+       integer ifc(klon,klev)
+       real pradj(klon)
+       real cldf(klon,klev), cldq(klon,klev)
+       real ftadj(klon,klev), fqadj(klon,klev)
+
+c sb --
+
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c
@@ -121,4 +150,5 @@
 #include "FCTTRE.h"
  
+      qcond_incld(:,:) = 0.
 c
 c$$$      print*,'debut conema'
@@ -190,5 +220,8 @@
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c 
-      CALL convect(dtime, em_t, em_q, em_qs,em_u ,em_v ,
+
+c     print*,'avant convect i=',i
+      CALL convect3(dtime,epmax,ok_adj_ema,
+     .              em_t, em_q, em_qs,em_u ,em_v ,
      .              em_tra, em_p, em_ph,
      .              klev, klev+1, klev-1,ntra, dtime, iflag,
@@ -197,10 +230,35 @@
      .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
      .              em_work1, em_work2,emmip,emMke,emMa,Ment,
-c SB 11sept98     .              Qent,TPS,TLS,SIJ)
-c 19oct98     .              Qent,TPS,TLS,SIJ,em_CAPE,em_TVP)
      .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
-     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr)
+     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
+     .  em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
+c     print*,'apres convect '
 c
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c -- sb: Appel schema statistique de nuages couple a la convection
+c (Bony et Emanuel 2001):
+
+c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
+
+        iflagcon = 3
+c       CALL cv_thermo(iflagcon)
+
+c -- appel schema de nuages:
+
+c       CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
+c    i          ,em_p,em_ph,dtime,em_qcondc
+c    o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
+
+        do k = 1, klev 
+         cldf(i,k)  = em_cldf(k)  ! cloud fraction (0-1)
+         cldq(i,k)  = em_cldq(k)  ! in-cloud water content (kg/kg)
+         ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
+         fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
+         ifc(i,k)   = em_ifc(k)   ! flag convergence clouds_gno (1 ou 2)
+        enddo
+        pradj(i) = em_pradj       ! precip from LS supersat adj (mm/day)
+
+c sb --
 c
 c SB:
@@ -251,4 +309,5 @@
       snow(i) = 0.0
       cape(i) = em_CAPE
+      wd(i) = em_wd
       rflag(i) = float(iflag)
 c SB      kbas(i) = em_bas
@@ -257,4 +316,9 @@
       dplcldr(i) = em_dplcldr
       DO l = 1, klev
+         d_t2(i,l) = dtime * em_d_t2(l) 
+         d_q2(i,l) = dtime * em_d_q2(l)
+         d_u2(i,l) = dtime * em_d_u2(l)
+         d_v2(i,l) = dtime * em_d_v2(l)
+
          d_t(i,l) = dtime * em_d_t(l) 
          d_q(i,l) = dtime * em_d_q(l)
@@ -273,6 +337,31 @@
          dtvpdt1(i,l) = em_dtvpdt1(l)
          dtvpdq1(i,l) = em_dtvpdq1(l)
+
+         if (iflag_clw.eq.0) then
+            qcond_incld(i,l) = em_qcondc(l)
+         else if (iflag_clw.eq.1) then
+            qcond_incld(i,l) = em_qcond(l)
+         endif
       ENDDO
   999 CONTINUE
+
+c   On calcule une eau liquide diagnostique en fonction de la 
+c  precip.
+      if ( iflag_clw.eq.2 ) then
+      do l=1,klev
+         do i=1,klon
+            if (ktop(i)-kbas(i).gt.0.and.
+     s         l.ge.kbas(i).and.l.le.ktop(i)) then
+               qcond_incld(i,l)=rain(i)*8.e4
+c    s         *(pplay(i,l      )-paprs(i,ktop(i)+1))
+     s         /(pplay(i,kbas(i))-pplay(i,ktop(i)))
+c    s         **2
+            else
+               qcond_incld(i,l)=0.
+            endif
+         enddo
+         print*,'l=',l,',   qcond_incld=',qcond_incld(1,l)
+      enddo
+      endif
  
 
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 373)
@@ -3,9 +3,13 @@
 !
 
-  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan)
+  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
+ &                     fact_cldcon, facttemps,ok_newmicro,iflag_cldcon, &
+ &                     ratqsbas,ratqshaut)
 
    use IOIPSL
    implicit none
 
+#include "conema3.h"
+#include "fisrtilp.inc"
 !
 ! Configuration de la "physique" de LMDZ a l'aide de la fonction
@@ -26,6 +30,8 @@
 ! Sortie:
   character (len = 6)  :: ocean
-  logical              :: ok_veget
+  logical              :: ok_veget, ok_newmicro
   logical              :: ok_journe, ok_mensuel, ok_instan
+  real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
+  integer              :: iflag_cldcon
 ! Local
   integer              :: numout = 6
@@ -78,4 +84,153 @@
   ok_instan = .false.
   call getin('OK_instan', ok_instan)
+!!
+!! KE
+!
+!Config Key  = epmax
+!Config Desc = Efficacite precip
+!Config Def  = 0.993
+!Config Help = 
+!
+  epmax = .993
+  call getin('epmax', epmax)
+!
+!Config Key  = ok_adj_ema
+!Config Desc =  
+!Config Def  = false
+!Config Help = 
+!
+  ok_adj_ema = .false.
+  call getin('ok_adj_ema',ok_adj_ema)
+!
+!Config Key  = iflag_clw
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  iflag_clw = 0
+  call getin('iflag_clw',iflag_clw)
+!
+!Config Key  = cld_lc_lsc 
+!Config Desc =  
+!Config Def  = 2.6e-4
+!Config Help = 
+!
+  cld_lc_lsc = 2.6e-4
+  call getin('cld_lc_lsc',cld_lc_lsc)
+!
+!Config Key  = cld_lc_con
+!Config Desc =  
+!Config Def  = 2.6e-4
+!Config Help = 
+!
+  cld_lc_con = 2.6e-4
+  call getin('cld_lc_con',cld_lc_con)
+!
+!Config Key  = cld_tau_lsc
+!Config Desc =  
+!Config Def  = 3600.
+!Config Help = 
+!
+  cld_tau_lsc = 3600.
+  call getin('cld_tau_lsc',cld_tau_lsc)
+!
+!Config Key  = cld_tau_con
+!Config Desc =  
+!Config Def  = 3600.
+!Config Help = 
+!
+  cld_tau_con = 3600.
+  call getin('cld_tau_con',cld_tau_con)
+!
+!Config Key  = ffallv_lsc
+!Config Desc =  
+!Config Def  = 1.
+!Config Help = 
+!
+  ffallv_lsc = 1.
+  call getin('ffallv_lsc',ffallv_lsc)
+!
+!Config Key  = ffallv_con
+!Config Desc =  
+!Config Def  = 1.
+!Config Help = 
+!
+  ffallv_con = 1.
+  call getin('ffallv_con',ffallv_con)
+!
+!Config Key  = coef_eva
+!Config Desc =  
+!Config Def  = 2.e-5
+!Config Help = 
+!
+  coef_eva = 2.e-5
+  call getin('coef_eva',coef_eva)
+!
+!Config Key  = reevap_ice
+!Config Desc =  
+!Config Def  = .false.
+!Config Help = 
+!
+  reevap_ice = .false.
+  call getin('reevap_ice',reevap_ice)
+!
+!Config Key  = iflag_cldcon 
+!Config Desc =  
+!Config Def  = 1
+!Config Help = 
+!
+  iflag_cldcon = 1
+  call getin('iflag_cldcon',iflag_cldcon)
+
+!
+!Config Key  = iflag_pdf 
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  iflag_pdf = 0
+  call getin('iflag_pdf',iflag_pdf)
+!
+!Config Key  = fact_cldcon
+!Config Desc =  
+!Config Def  = 0.375
+!Config Help = 
+!
+  fact_cldcon = 0.375
+  call getin('fact_cldcon',fact_cldcon)
+
+!
+!Config Key  = facttemps
+!Config Desc =  
+!Config Def  = 1.e-4
+!Config Help = 
+!
+  facttemps = 1.e-4
+  call getin('facttemps',facttemps)
+
+!
+!Config Key  = ok_newmicro
+!Config Desc =  
+!Config Def  = .true.
+!Config Help = 
+!
+  ok_newmicro = .true.
+  call getin('ok_newmicro',ok_newmicro)
+!
+!Config Key  = ratqsbas
+!Config Desc =  
+!Config Def  = 0.01
+!Config Help = 
+!
+  ratqsbas = 0.01
+  call getin('ratqsbas',ratqsbas)
+!
+!Config Key  = ratqshaut
+!Config Desc =  
+!Config Def  = 0.3
+!Config Help = 
+!
+  ratqshaut = 0.3
+  call getin('ratqshaut',ratqshaut)
 
 !
@@ -90,5 +245,22 @@
   write(numout,*)' Sortie mensuelle = ', ok_mensuel
   write(numout,*)' Sortie instantanee = ', ok_instan
-
+  write(numout,*)' epmax = ', epmax
+  write(numout,*)' ok_adj_ema = ', ok_adj_ema
+  write(numout,*)' iflag_clw = ', iflag_clw
+  write(numout,*)' cld_lc_lsc = ', cld_lc_lsc
+  write(numout,*)' cld_lc_con = ', cld_lc_con
+  write(numout,*)' cld_tau_lsc = ', cld_tau_lsc
+  write(numout,*)' cld_tau_con = ', cld_tau_con
+  write(numout,*)' ffallv_lsc = ', ffallv_lsc
+  write(numout,*)' ffallv_con = ', ffallv_con
+  write(numout,*)' coef_eva = ', coef_eva
+  write(numout,*)' reevap_ice = ', reevap_ice
+  write(numout,*)' iflag_pdf = ', iflag_pdf
+  write(numout,*)' iflag_cldcon = ', iflag_cldcon
+  write(numout,*)' fact_cldcon = ', fact_cldcon
+  write(numout,*)' facttemps = ', facttemps
+  write(numout,*)' ok_newmicro = ',ok_newmicro 
+  write(numout,*)' ratqsbas = ',ratqsbas 
+  write(numout,*)' ratqshaut = ',ratqshaut 
 
   return
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F	(revision 373)
@@ -1,10 +1,13 @@
 c $Header$
-      SUBROUTINE CONVECT
-     *    (DTIME,  T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
+      SUBROUTINE CONVECT3	
+     *    (DTIME,EPMAX,ok_adj,
+     *     T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
      *     ND,       NDP1,     NL, NTRA,  DELT,  IFLAG,
      *     FT, FR, FU,  FV,  FTRA,  PRECIP,
      *     icb,inb,   upwd,dnwd,dnwd0,SIG, W0,Mike,Mke,
      *     Ma,MENTS,QENTS,TPS,TLS,SIGIJ,CAPE,TVP,PBASE,BUOYBASE,
-     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
+cccc     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
+     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR,   ! sbl
+     *     FT2,FR2,FU2,FV2,WD,QCOND,QCONDC)   ! sbl
 C
 C    ***  THE PARAMETER NA SHOULD IN GENERAL EQUAL ND   ***
@@ -19,4 +22,6 @@
       integer NTRAC
       PARAMETER (NTRAC=nqmx-2)
+      REAL DELTAC              ! cld
+      PARAMETER (DELTAC=0.01)  ! cld
 
       INTEGER NENT(NA)
@@ -34,4 +39,7 @@
       REAL T(NA),RR(NA)
 C
+      REAL FT2(ND),FR2(ND),FU2(ND),FV2(ND) ! added sbl
+      REAL U1(ND),V1(ND) ! added sbl
+C
       REAL BUOY(NA)     !  Lifted parcel buoyancy
       REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
@@ -39,5 +47,9 @@
       REAL CLW_NEW(NA),QI(NA)
 C
-      LOGICAL ICE_CONV
+      REAL WD, BETAD ! for gust factor (sb)
+      REAL QCONDC(ND)  ! interface cld param (sb)
+      REAL QCOND(ND),NQCOND(NA),WA(NA),MAA(NA),SIGA(NA),AXC(NA) ! cld
+C
+      LOGICAL ICE_CONV,ok_adj
       PARAMETER (ICE_CONV=.TRUE.)
  
@@ -110,9 +122,22 @@
          FU(I)=0.0
          FV(I)=0.0
+
+         FT2(I)=0.0
+         FR2(I)=0.0
+         FU2(I)=0.0
+         FV2(I)=0.0
+
          DO 4 J=1,NTRA
           FTRA(I,J)=0.0
     4    CONTINUE
+
+         QCONDC(I)=0.0  ! cld
+         QCOND(I)=0.0   ! cld
+         NQCOND(I)=0.0  ! cld
+
          T(I)=T1(I)
          RR(I)=R1(I)
+         U1(I)=U(I) ! added sbl
+         V1(I)=V(I) ! added sbl
     5 CONTINUE
       DO 7 I=1,NL
@@ -137,4 +162,5 @@
  
       PRECIP=0.0
+      WD=0.0 ! sb
       IFLAG=1
 C
@@ -160,5 +186,5 @@
       SPFAC=0.15
 c sb:
-      EPMAX=0.993 ! precip efficiency less than unity
+c     EPMAX=0.993 ! precip efficiency less than unity
 c      EPMAX=1. ! precip efficiency less than unity
 C
@@ -196,4 +222,5 @@
 C
       NOPT=0
+c!      NOPT=1 ! sbl
 C
 C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
@@ -202,57 +229,67 @@
 C     ***                BOUNDARY LAYER SCHEME !!!               ***
 C
-c test sb      DO 30 I=NL-1,1,-1
-c test sb         JN=0
-c test sb         DO 10 J=I+1,NL
-c test sb   10    IF(TH(J).LT.TH(I))JN=J
-c test sb         IF(JN.EQ.0)GOTO 30
-c test sb         AHM=0.0
-c test sb         RM=0.0
-c test sb         UM=0.0
-c test sb         VM=0.0
-c test sb         DO K=1,NTRA
-c test sb          TRATM(K)=0.0
-c test sb         END DO
-c test sb         DO 15 J=I,JN
-c test sb          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
-c test sb          RM=RM+RR(J)*(PH(J)-PH(J+1))
-c test sb          UM=UM+U(J)*(PH(J)-PH(J+1))
-c test sb          VM=VM+V(J)*(PH(J)-PH(J+1))
-c test sb          DO K=1,NTRA
-c test sb           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
-c test sb          END DO
-c test sb   15    CONTINUE
-c test sb         DPHINV=1./(PH(I)-PH(JN+1))
-c test sb         RM=RM*DPHINV
-c test sb         UM=UM*DPHINV
-c test sb         VM=VM*DPHINV
-c test sb         DO K=1,NTRA
-c test sb          TRATM(K)=TRATM(K)*DPHINV
-c test sb         END DO
-c test sb         A2=0.0
-c test sb         DO 20 J=I,JN
-c test sb            RR(J)=RM
-c test sb          U(J)=UM
-c test sb          V(J)=VM
-c test sb          DO K=1,NTRA
-c test sb           TRA(J,K)=TRATM(K)
-c test sb          END DO
-c test sb            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
-c test sb            X=(0.001*P(J))**RDCP
-c test sb            T(J)=X
-c test sb            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
-c test sb   20    CONTINUE
-c test sb         DO 25 J=I,JN
-c test sb            TH(J)=AHM/A2
-c test sb            T(J)=T(J)*TH(J)
-c test sb   25    CONTINUE
-c test sb   30 CONTINUE
-C
-C   ***   RESET INPUT ARRAYS IF NOPT .NE. 0   ***
-C
-      IF (NOPT.NE.0)THEN
+      IF (ok_adj) THEN ! added sbl
+
+      DO 30 I=NL-1,1,-1
+         JN=0
+         DO 10 J=I+1,NL
+   10    IF(TH(J).LT.TH(I))JN=J
+         IF(JN.EQ.0)GOTO 30
+         AHM=0.0
+         RM=0.0
+         UM=0.0
+         VM=0.0
+         DO K=1,NTRA
+          TRATM(K)=0.0
+         END DO
+         DO 15 J=I,JN
+          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
+          RM=RM+RR(J)*(PH(J)-PH(J+1))
+          UM=UM+U(J)*(PH(J)-PH(J+1))
+          VM=VM+V(J)*(PH(J)-PH(J+1))
+          DO K=1,NTRA
+           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
+          END DO
+   15    CONTINUE
+         DPHINV=1./(PH(I)-PH(JN+1))
+         RM=RM*DPHINV
+         UM=UM*DPHINV
+         VM=VM*DPHINV
+         DO K=1,NTRA
+          TRATM(K)=TRATM(K)*DPHINV
+         END DO
+         A2=0.0
+         DO 20 J=I,JN
+            RR(J)=RM
+          U(J)=UM
+          V(J)=VM
+          DO K=1,NTRA
+           TRA(J,K)=TRATM(K)
+          END DO
+            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
+            X=(0.001*P(J))**RDCP
+            T(J)=X
+            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
+   20    CONTINUE
+         DO 25 J=I,JN
+            TH(J)=AHM/A2
+            T(J)=T(J)*TH(J)
+   25    CONTINUE
+   30 CONTINUE
+
+      ENDIF ! added sbl
+C
+C   ***   RESET INPUT ARRAYS IF ok_adj 0   ***
+C
+      IF (ok_adj)THEN
          DO 35 I=1,ND
-            T1(I)=T(I)
-            R1(I)=RR(I)
+
+           FT2(I)=(T(I)-T1(I))/DELT  ! sbl
+           FR2(I)=(RR(I)-R1(I))/DELT  ! sbl
+           FU2(I)=(U(I)-U1(I))/DELT  ! sbl
+           FV2(I)=(V(I)-V1(I))/DELT  ! sbl
+
+c!            T1(I)=T(I)      ! commente sbl
+c!            R1(I)=RR(I)     ! commente sbl
    35    CONTINUE
       END IF
@@ -401,4 +438,12 @@
      $        +TV(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
 C
+c test sb:
+c@      write(*,*) '++++++++++++++++++++++++++++++++++++++++'
+c@      write(*,*) 'plcl,dpbas,tvpbas,tvbas,tvp(icb),tvp(icb1)
+c@     :             ,tv(icb),tv(icb1)'
+c@      write(*,*) plcl,dpbase,tvpbase,tvbase,tvp(icb)
+c@     L          ,tvp(icb+1),tv(icb),tv(icb+1)
+c@      write(*,*) '++++++++++++++++++++++++++++++++++++++++'
+c fin test sb
       BUOYBASE = TVPBASE - TVBASE
 C
@@ -418,4 +463,16 @@
 c sb3d      print *, 'P ',(p(i),i=1,nl)
 c sb3d      print *,'ICB ',icb
+c test sb:
+c@      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
+c@      write(*,*) 'icb,icbs,inb,buoybase:'
+c@      write(*,*) icb,icb+1,inb,buoybase
+c@      write(*,*) 'k,tvp,tv,tp,buoy,ep: '
+c@      do k=1,nl
+c@      write(*,*) k,tvp(k),tv(k),tp(k),buoy(k),ep(k)
+c@      enddo
+c@      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
+c fin test sb
+
+
 C
 C   ***   MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE  ***
@@ -552,4 +609,11 @@
 Cjyg2
 c sb3d         print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape
+c test sb:
+c@       write(*,*) '############################################'
+c@         write(*,*) 'cape,rrd,buoy,deltap,p,pbase,ph:'
+c@     :     ,cape,rd,buoy(i-1),deltap,p(i-1),pbase,ph(i)
+c@       write(*,*) '############################################'
+
+c fin test sb
          CAPE=AMAX1(0.0,CAPE)
 C
@@ -573,4 +637,13 @@
          AMU=0.5*(SIG(I)+SIGOLD)*W
          M(I)=AMU*0.007*P(I)*(PH(I)-PH(I+1))/TV(I)
+
+c --------- test sb:
+c       write(*,*) '############################################'
+c       write(*,*) 'k,amu,buoy(k-1),deltap,w,beta,fac,cape,w0(k)'
+c       write(*,*) i,amu,buoy(i-1),deltap
+c     :           ,w,beta,fac,cape,w0(i)
+c       write(*,*) '############################################'
+c ---------
+
          W0(I)=W
    98 CONTINUE
@@ -750,4 +823,11 @@
  
 ***********************************************************
+c--- test sb:
+c@       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
+c@       write(*,*) 'inb,m(inb),ment(inb,inb),sigij(inb,inb):'
+c@       write(*,*) inb,m(inb),ment(inb,inb),sigij(inb,inb)
+c@       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
+c---
+
  
  
@@ -936,4 +1016,13 @@
 C
         PRECIP=WT(1)*SIGD*WATER(1)*8640.0
+
+c sb  ***  Calculate downdraft velocity scale and surface temperature and  ***
+c sb  ***                    water vapor fluctuations                      ***
+c sb		(inspire de convect 4.3)
+
+c       BETAD=10.0         
+       BETAD=5.0         
+       WD=BETAD*ABS(MP(ICB))*0.01*RD*T(ICB)/(SIGD*P(ICB))
+
   405   CONTINUE
 C
@@ -1021,4 +1110,7 @@
          FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
          FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
+C (saturated updrafts resulting from mixing)      ! cld   
+         QCOND(I)=QCOND(I)+(ELIJ(K,I)-AWAT)       ! cld
+         NQCOND(I)=NQCOND(I)+1.                   ! cld
          DO J=1,NTRA
           FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
@@ -1045,4 +1137,17 @@
      1    MP(I)*(TRAP(I,J)-TRAP(I-1,J)))
         END DO
+C (saturated downdrafts resulting from mixing)    ! cld
+        DO K=I+1,INB                              ! cld
+         QCOND(I)=QCOND(I)+ELIJ(K,I)              ! cld
+         NQCOND(I)=NQCOND(I)+1.                   ! cld
+        ENDDO                                     ! cld
+C (particular case: no detraining level is found) ! cld
+        IF (NENT(I).EQ.0) THEN                    ! cld
+         QCOND(I)=QCOND(I)+(1-EP(I))*CLW(I)       ! cld
+         NQCOND(I)=NQCOND(I)+1.                   ! cld
+        ENDIF                                     ! cld
+        IF (NQCOND(I).NE.0.) THEN                 ! cld
+         QCOND(I)=QCOND(I)/NQCOND(I)              ! cld
+        ENDIF                                     ! cld
   500   CONTINUE
  
@@ -1054,4 +1159,14 @@
 C   ***          INTEGRATED ENTHALPY AND WATER TENDENCIES         ***
 C
+c test sb:
+c@      write(*,*) '--------------------------------------------'
+c@      write(*,*) 'inb,ft,hp,h,t,rr,qent,ment,water,waterp,wt,mp,b'
+c@      write(*,*) inb,ft(inb),hp(inb),h(inb)
+c@     :   ,t(inb),rr(inb),qent(inb,inb)
+c@     :   ,ment(inb,inb),water(inb)
+c@     :   ,water(inb+1),wt(inb),mp(inb),b(inb)
+c@      write(*,*) '--------------------------------------------'
+c fin test sb:
+
       AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)*
      1    (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)*
@@ -1187,4 +1302,37 @@
        Mke(i)=upwd(i)+dnwd(i)
        enddo
+
+C
+C   *** Diagnose the in-cloud mixing ratio   ***              ! cld
+C   ***           of condensed water         ***              ! cld
+C                                                             ! cld
+       DO I=1,ND                                              ! cld
+        MAA(I)=0.0                                            ! cld
+        WA(I)=0.0                                             ! cld
+        SIGA(I)=0.0                                           ! cld
+       ENDDO                                                  ! cld
+       DO I=NK,INB                                            ! cld
+       DO K=I+1,INB+1                                         ! cld
+        MAA(I)=MAA(I)+M(K)                                    ! cld
+       ENDDO                                                  ! cld
+       ENDDO                                                  ! cld
+       DO I=ICB,INB-1                                         ! cld
+        AXC(I)=0.                                             ! cld
+        DO J=ICB,I                                            ! cld
+         AXC(I)=AXC(I)+RD*(TVP(J)-TV(J))*(PH(J)-PH(J+1))/P(J) ! cld
+        ENDDO                                                 ! cld
+        IF (AXC(I).GT.0.0) THEN                               ! cld
+         WA(I)=SQRT(2.*AXC(I))                                ! cld
+        ENDIF                                                 ! cld
+       ENDDO                                                  ! cld
+       DO I=1,NL                                              ! cld
+        IF (WA(I).GT.0.0)                                     ! cld
+     1    SIGA(I)=MAA(I)/WA(I)*RD*TVP(I)/P(I)/100./DELTAC     ! cld
+        SIGA(I) = MIN(SIGA(I),1.0)                            ! cld
+        QCONDC(I)=SIGA(I)*CLW(I)*(1.-EP(I))                   ! cld
+     1          + (1.-SIGA(I))*QCOND(I)                       ! cld
+       ENDDO                                                  ! cld
+
+
 c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c$$$         call writeg1d(1,klev,ma,'ma  ','ma  ')
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F	(revision 373)
@@ -1,5 +1,10 @@
-      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,
+c $Header$
+c
+      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
      s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
-     s                   prfl, psfl)
+     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
+     s                   frac_impa, frac_nucl,
+     s                   prfl, psfl, rhcl)
+
 c
       IMPLICIT none
@@ -14,4 +19,6 @@
 #include "dimphy.h"
 #include "YOMCST.h"
+#include "tracstoke.h"
+#include "fisrtilp.h"
 c
 c Arguments:
@@ -27,8 +34,23 @@
       REAL rneb(klon,klev) ! fraction nuageuse
       REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
+      REAL rhcl(klon,klev) ! humidite relative en ciel clair
       REAL rain(klon) ! pluies (mm/s)
       REAL snow(klon) ! neige (mm/s)
       REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
       REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
+cAA
+c Coeffients de fraction lessivee : pour OFF-LINE
+c
+      REAL pfrac_nucl(klon,klev)
+      REAL pfrac_1nucl(klon,klev)
+      REAL pfrac_impa(klon,klev)
+c
+c Fraction d'aerosols lessivee par impaction et par nucleation
+c POur ON-LINE
+c
+      REAL frac_impa(klon,klev)
+      REAL frac_nucl(klon,klev)
+      real zct(klon),zcl(klon)
+cAA
 c
 c Options du programme:
@@ -36,24 +58,15 @@
       REAL seuil_neb ! un nuage existe vraiment au-dela
       PARAMETER (seuil_neb=0.001)
-      REAL ct ! inverse du temps pour qu'un nuage precipite
-      PARAMETER (ct=1./1800.)
-      REAL cl ! seuil de precipitation
-      PARAMETER (cl=2.6e-4)
-ccc      PARAMETER (cl=2.3e-4)
-ccc      PARAMETER (cl=2.0e-4)
+
       INTEGER ninter ! sous-intervals pour la precipitation
       PARAMETER (ninter=5)
       LOGICAL evap_prec ! evaporation de la pluie
       PARAMETER (evap_prec=.TRUE.)
-      REAL coef_eva
-      PARAMETER (coef_eva=2.0E-05)
-      LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
-      REAL ratqs ! determine la largeur de distribution de vapeur
-      PARAMETER (calcrat=.TRUE.)
-      REAL zx_min, rat_max
-      PARAMETER (zx_min=1.0, rat_max=0.01)
-      REAL zx_max, rat_min
-      PARAMETER (zx_max=0.1, rat_min=0.3)
-      REAL zx
+      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
+      logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
+
+      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
+      real erf
 c
       LOGICAL cpartiel ! condensation partielle
@@ -64,5 +77,5 @@
 c Variables locales:
 c
-      INTEGER i, k, n
+      INTEGER i, k, n, kk
       REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
       REAL zrfl(klon), zrfln(klon), zqev, zqevt
@@ -76,18 +89,34 @@
       SAVE appel1er
 c
+c---------------------------------------------------------------
+c
+cAA Variables traceurs:
+cAA  Provisoire !!! Parametres alpha du lessivage
+cAA  A priori on a 4 scavenging # possibles
+c
+      REAL a_tr_sca(4)
+      save a_tr_sca
+c
+c Variables intermediaires
+c
+      REAL zalpha_tr
+      REAL zfrac_lessi
+      REAL zprec_cond(klon)
+cAA
+c---------------------------------------------------------------
+c
 c Fonctions en ligne:
 c
-      REAL fallv ! vitesse de chute pour crystaux de glace
+      REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
       REAL zzz
 #include "YOETHF.h"
 #include "FCTTRE.h"
-      fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
-ccc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
-ccc      fallv (zzz) = 3.29 * ((zzz)**0.16)
+      fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
+      fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
 c
       DATA appel1er /.TRUE./
-c
+
       IF (appel1er) THEN
-         PRINT*, 'fisrtilp, calcrat:', calcrat
+c
          PRINT*, 'fisrtilp, ninter:', ninter
          PRINT*, 'fisrtilp, evap_prec:', evap_prec
@@ -96,10 +125,32 @@
           PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
           PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
-          CALL abort
+c         CALL abort
          ENDIF
          appel1er = .FALSE.
-      ENDIF
-c
+c
+cAA initialiation provisoire
+       a_tr_sca(1) = -0.5
+       a_tr_sca(2) = -0.5
+       a_tr_sca(3) = -0.5
+       a_tr_sca(4) = -0.5
+c
+cAA Initialisation a 1 des coefs des fractions lessivees 
+c
+      DO k = 1, klev
+       DO i = 1, klon
+          pfrac_nucl(i,k)=1.
+          pfrac_1nucl(i,k)=1.
+          pfrac_impa(i,k)=1.
+       ENDDO 
+      ENDDO 
+
+      ENDIF          !  test sur appel1er
+c
+cMAf Initialisation a 0 de zoliq
+       DO i = 1, klon
+          zoliq(i)=0.
+       ENDDO 
 c Determiner les nuages froids par leur temperature
+c  nexpo regle la raideur de la transition eau liquide / eau glace.
 c
       ztglace = RTT - 15.0
@@ -115,4 +166,5 @@
       ENDDO
       ENDDO
+
       DO k = 1, klev
       DO i = 1, klon
@@ -122,4 +174,6 @@
          rneb(i,k) = 0.0
          radliq(i,k) = 0.0
+         frac_nucl(i,k) = 1. 
+         frac_impa(i,k) = 1. 
       ENDDO
       ENDDO
@@ -136,7 +190,17 @@
       ENDDO
 c
+c
+cAA Pour plus de securite 
+
+      zalpha_tr   = 0.
+      zfrac_lessi = 0.
+
+cAA----------------------------------------------------------
+c
 c Boucle verticale (du haut vers le bas)
 c
       DO 9999 k = klev, 1, -1
+c
+cAA----------------------------------------------------------
 c
       DO i = 1, klon
@@ -171,4 +235,11 @@
          zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
      .                            /RG/dtime
+
+c pour la glace, on réévapore toute la précip dans la couche du dessous
+c la glace venant de la couche du dessus est simplement dans la couche
+c 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
@@ -210,22 +281,73 @@
 c
       IF (cpartiel) THEN
-         DO i = 1, klon
-c
-            zx = pplay(i,k)/paprs(i,1)
-            zx = (zx_max-zx)/(zx_max-zx_min)
-            zx = MIN(MAX(zx,0.0),1.0)
-            zx = zx * zx * zx
-            ratqs = zx * (rat_max-rat_min) + rat_min
-            IF (.NOT.calcrat) ratqs=0.05
-ccc            IF (.NOT.calcrat) ratqs=0.2
-c
-            zdelq = ratqs * zq(i)
+
+c        print*,'Dans partiel k=',k
+c
+c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
+c   nuageuse a partir des PDF de Sandrine Bony.
+c   rneb  : fraction nuageuse
+c   zqn   : eau totale dans le nuage
+c   zcond : eau condensee moyenne dans la maille.
+c           on prend en compte le réchauffement qui diminue la partie condensee
+c
+c   Version avec les raqts
+
+         if (iflag_pdf.eq.0) then
+
+           do i=1,klon
+            zdelq = min(ratqs(i,k),0.99) * zq(i)
             rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
             zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
+           enddo
+
+         else
+c
+c   Version avec les nouvelles PDFs.
+           do i=1,klon
+              if(zq(i).lt.1.e-15) then
+                print*,'ZQ(',i,',',k,')=',zq(i)
+                zq(i)=1.e-15
+              endif
+           enddo
+           do i=1,klon
+            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
+            
+           enddo
+
+        endif ! iflag_pdf
+
+         do i=1,klon
             IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
             IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
             rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
-         ENDDO
+c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
+c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
+c  la convection.
+c  ATTENTION !!! Il va falloir verifier tout ca.
+            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+c           print*,'ZDQS ',zdqs(i)
+c--Olivier
+            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
+            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
+c--fin
+           ENDDO
       ELSE
          DO i = 1, klon
@@ -262,9 +384,24 @@
       DO i = 1, klon
       IF (rneb(i,k).GT.0.0) THEN
-         zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
-     .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
          zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-         zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
-     .              *fallv(zrhol(i)) * zfice(i)
+
+         if (ptconv(i,k)) then
+            zcl(i)=cld_lc_con
+            zct(i)=1./cld_tau_con
+         else
+            zcl(i)=cld_lc_lsc
+            zct(i)=1./cld_tau_lsc
+         endif
+c  quantité d'eau à élminier.
+         zchau(i) = zct(i)*dtime/FLOAT(ninter) * zoliq(i)
+     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))**2)) *(1.-zfice(i))
+c  meme chose pour la glace.
+         if (ptconv(i,k)) then
+            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
+     .              *fallvc(zrhol(i)) * zfice(i)
+         else
+            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
+     .              *fallvs(zrhol(i)) * zfice(i)
+         endif
          ztot(i) = zchau(i) + zfroi(i)
          IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
@@ -296,5 +433,49 @@
       ENDDO
 c
+cAA--------------- Calcul du lessivage stratiforme  -------------
+
+      DO i = 1,klon
+c
+         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
+     .                * (paprs(i,k)-paprs(i,k+1))/RG
+         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+cAA lessivage nucleation LMD5 dans la couche elle-meme
+            if (t(i,k) .GE. ztglace) THEN
+               zalpha_tr = a_tr_sca(3)
+            else
+               zalpha_tr = a_tr_sca(4)
+            endif
+            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
+c
+c nucleation avec un facteur -1 au lieu de -0.5
+            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
+            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+         ENDIF
+c
+      ENDDO      ! boucle sur i
+c
+cAA Lessivage par impaction dans les couches en-dessous
+      DO kk = k-1, 1, -1
+        DO i = 1, klon
+          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+            if (t(i,kk) .GE. ztglace) THEN
+              zalpha_tr = a_tr_sca(1)
+            else
+              zalpha_tr = a_tr_sca(2)
+            endif
+            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
+            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
+          ENDIF
+        ENDDO
+      ENDDO
+c
+cAA----------------------------------------------------------
+c                     FIN DE BOUCLE SUR K   
  9999 CONTINUE
+c
+cAA-----------------------------------------------------------
 c
 c Pluie ou neige au sol selon la temperature de la 1ere couche
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 373)
@@ -7,5 +7,5 @@
      .           fder,radsol,frugs,agesno,clesphy0,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,tabcntr0,
-     .           t_ancien,q_ancien,ancien_ok)
+     .           t_ancien,q_ancien,ancien_ok, rnebcon, ratqs,clwcon)
       IMPLICIT none
 c======================================================================
@@ -54,4 +54,5 @@
 
       REAL t_ancien(klon,klev), q_ancien(klon,klev)
+      real rnebcon(klon,klev),clwcon(klon,klev),ratqs(klon,klev)
       LOGICAL ancien_ok
 
@@ -1170,4 +1171,88 @@
       ENDIF
 c
+      ierr = NF_INQ_VARID (nid, "CLWCON", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, "phyetat0: Le champ CLWCON est absent"
+         PRINT*, "Depart legerement fausse. Mais je continue"
+         clwcon = 0.
+      ELSE
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, clwcon)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, clwcon)
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <CLWCON>"
+            CALL abort
+         ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      xmin = MINval(clwcon)
+      xmax = MAXval(clwcon)
+      PRINT*,'Eau liquide convective (ecart-type) clwcon:', xmin, xmax
+c
+      ierr = NF_INQ_VARID (nid, "RNEBCON", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, "phyetat0: Le champ RNEBCON est absent"
+         PRINT*, "Depart legerement fausse. Mais je continue"
+         rnebcon = 0.
+      ELSE
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rnebcon)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, rnebcon)
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <RNEBCON>"
+            CALL abort
+         ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      xmin = MINval(rnebcon)
+      xmax = MAXval(rnebcon)
+      PRINT*,'Nebulosite convective (ecart-type) rnebcon:', xmin, xmax
+
+c
+      ierr = NF_INQ_VARID (nid, "QANCIEN", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, "phyetat0: Le champ <QANCIEN> est absent"
+         PRINT*, "Depart legerement fausse. Mais je continue"
+         ancien_ok = .FALSE.
+      ELSE
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q_ancien)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, q_ancien)
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <QANCIEN>"
+            CALL abort
+         ENDIF
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "RATQS", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, "phyetat0: Le champ <RATQS> est absent"
+         PRINT*, "Depart legerement fausse. Mais je continue"
+         ratqs = 0.
+      ELSE
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ratqs)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, ratqs)
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <RATQS>"
+            CALL abort
+         ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      xmin = MINval(ratqs)
+      xmax = MAXval(ratqs)
+      PRINT*,'(ecart-type) ratqs:', xmin, xmax
+c
 c Fermer le fichier:
 c
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 373)
@@ -7,5 +7,5 @@
      .           radsol,frugs,agesno,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
-     .           t_ancien, q_ancien)
+     .           t_ancien, q_ancien, rnebcon, ratqs, clwcon)
       IMPLICIT none
 c======================================================================
@@ -53,4 +53,5 @@
       REAL pctsrf(klon, nbsrf)
       REAL t_ancien(klon,klev), q_ancien(klon,klev)
+      real clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
 c
       INTEGER nid, nvarid, idim1, idim2, idim3
@@ -96,6 +97,4 @@
       IF(     ok_orolf ) tab_cntrl(11 ) = 1.
 
-      tab_cntrl(13) = dayref
-      tab_cntrl(14) = anneeref
       tab_cntrl(13) = day_end
       tab_cntrl(14) = annee_ref
@@ -670,4 +669,49 @@
 #endif
 c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "CLWCON", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "CLWCON", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Eau liquide convective")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,clwcon)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,clwcon)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RNEBCON", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RNEBCON", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Nebulosite convective")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rnebcon)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rnebcon)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RATQS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RATQS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Ratqs")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ratqs)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ratqs)
+#endif
+c
 c
       ierr = NF_CLOSE(nid)
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 367)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 373)
@@ -388,5 +388,5 @@
       EXTERNAL conlmd    ! convection (schema LMD)
 cKE43
-      EXTERNAL conema  ! convect4.3
+      EXTERNAL conema3  ! convect4.3
       EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
 cAA 
@@ -412,4 +412,9 @@
 c Variables locales
 c
+      real clwcon(klon,klev),rnebcon(klon,klev)
+      real clwcon0(klon,klev),rnebcon0(klon,klev)
+      save rnebcon, clwcon
+
+      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
       REAL dialiq(klon,klev)  ! eau liquide nuageuse
       REAL diafra(klon,klev)  ! fraction nuageuse
@@ -464,6 +469,7 @@
 c
       REAL za, zb
-      REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
-      INTEGER i, k, iq, nsrf, ll
+      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
+      real zqsat(klon,klev)
+      INTEGER i, k, iq, ig, j, nsrf, ll
       REAL t_coup
       PARAMETER (t_coup=234.0)
@@ -534,7 +540,20 @@
       REAL d_t_lif(klon,klev)
 
-      REAL ratqs(klon,klev)
-      integer flag_ratqs
+      REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
+      real ratqsbas,ratqshaut
+      save ratqsbas,ratqshaut, ratqs
       real zpt_conv(klon,klev)
+
+c Parametres lies au nouveau schema de nuages (SB, PDF)
+      real fact_cldcon
+      real facttemps
+      logical ok_newmicro
+      save ok_newmicro
+      save fact_cldcon,facttemps
+
+      integer iflag_cldcon
+      save iflag_cldcon
+
+      logical ptconv(klon,klev)
 
 c
@@ -630,5 +649,6 @@
 c
          call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
-     .                  ok_instan)
+     .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
+     .                  iflag_cldcon,ratqsbas,ratqshaut)
 
          DO k = 2, nvm          ! pas de vegetation
@@ -655,5 +675,5 @@
      .       dlw,radsol,frugs,agesno,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
-     .       t_ancien, q_ancien, ancien_ok )
+     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
 
 c
@@ -1307,4 +1327,8 @@
      .                "ave(X)", zsto,zout)
 c
+         CALL histdef(nid_mth, "ducon", "Convection du", "m/s2",
+     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+c
          CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
      .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
@@ -1774,5 +1798,6 @@
       DO i = 1, klon
          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-         zlsdcp=RLSTT/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))
          zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
          zb = MAX(0.0,ql_seri(i,k))
@@ -1981,5 +2006,7 @@
           else
 
-          CALL conema (dtime,paprs,pplay,t_seri,q_seri,
+c          print*,'Avant conema OUI'
+          CALL conema3 (dtime,
+     .        paprs,pplay,t_seri,q_seri,
      .        u_seri,v_seri,tr_seri,nbtr,
      .        ema_work1,ema_work2,
@@ -1988,6 +2015,36 @@
      .        upwd,dnwd,dnwd0,bas,top,
      .        Ma,cape,tvp,rflag,
-     .       pbase
-     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
+     .        pbase
+     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
+     .        ,clwcon0)
+c          print*,'Apres conema3 '
+
+c Calculer l'humidite relative pour diagnostique
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         zx_t = t_seri(i,k)
+         IF (thermcep) THEN
+            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
+            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
+            zx_qs  = MIN(0.5,zx_qs)
+            zcor   = 1./(1.-retv*zx_qs)
+            zx_qs  = zx_qs*zcor
+         ELSE
+           IF (zx_t.LT.t_coup) THEN
+              zx_qs = qsats(zx_t)/pplay(i,k)
+           ELSE
+              zx_qs = qsatl(zx_t)/pplay(i,k)
+           ENDIF
+         ENDIF
+         zqsat(i,k)=zx_qs
+      ENDDO
+      ENDDO
+
+c   calcul des propriétés des nuages convectifs
+             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
+             call clouds_gno
+     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
+
           endif
           DO i = 1, klon
@@ -2054,10 +2111,10 @@
       IF (nqmax.GT.2) THEN !--melange convectif de traceurs
 c
-          IF (iflag_con .LT. 2 .AND.  iflag_con .GT. 4 ) THEN 
+          IF (iflag_con .NE. 2 .AND. debut) THEN 
               PRINT*, 'Pour l instant, seul conflx fonctionne ',
      $            'avec traceurs', iflag_con
               PRINT*,' Mettre iflag_con', 
-     $            ' = 2, 3 ou 4 dans run.def et repasser'
-              CALL abort
+     $            ' = 2 dans run.def et repasser'
+c              CALL abort
               ENDIF 
 c
@@ -2074,24 +2131,63 @@
       ENDDO
 
-c   RATQS
-      if (iflag_con.eq.2) then
-          flag_ratqs=0
+
+c-------------------------------------------------------------------------
+c  Caclul des ratqs
+c-------------------------------------------------------------------------
+
+c      print*,'calcul des ratqs'
+c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
+c   ----------------
+c   on ecrase le tableau ratqsc calcule par clouds_gno
+      if (iflag_cldcon.eq.1) then
+         do k=1,klev
+         do i=1,klon
+            if(ptconv(i,k)) then
+              ratqsc(i,k)=ratqsbas
+     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
+            else
+               ratqsc(i,k)=0.
+            endif
+         enddo
+         enddo
+      endif
+
+c   ratqs stables
+c   -------------
+      do k=1,klev
+         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
+     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)        
+      enddo
+
+
+c  ratqs final
+c  -----------
+      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
+c   les ratqs sont une conbinaison de ratqss et ratqsc
+c   ratqs final
+c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
+c   relaxation des ratqs
+         facttemps=exp(-pdtphys/1.e4)
+         ratqs(:,:)=max(ratqs(:,:)*facttemps,ratqss(:,:))
+         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
+c         print*,'calcul des ratqs fini'
       else
-          flag_ratqs=1
+c   on ne prend que le ratqs stable pour fisrtilp
+         ratqs(:,:)=ratqss(:,:)
       endif
-      call calcratqs (flag_ratqs,
-     I            paprs,pplay,q_seri,d_t_con,d_t_ajs
-     O           ,ratqs,zpt_conv)
+
+
 c
 c Appeler le processus de condensation a grande echelle
 c et le processus de precipitation
-c
-      CALL fisrtilp_tr(dtime,paprs,pplay,
-     .           t_seri, q_seri,ratqs,
+c-------------------------------------------------------------------------
+      CALL fisrtilp(dtime,paprs,pplay,
+     .           t_seri, q_seri,ptconv,ratqs,
      .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
      .           rain_lsc, snow_lsc,
      .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
      .           frac_impa, frac_nucl,
-     .           prfl, psfl)
+     .           prfl, psfl, rhcl)
+
       WHERE (rain_lsc < 0) rain_lsc = 0.
       WHERE (snow_lsc < 0) snow_lsc = 0.
@@ -2118,7 +2214,13 @@
       ENDIF
 c
-c Nuages diagnostiques:
-c
-      IF (iflag_con.EQ.2) THEN ! seulement pour Tiedtke
+c-------------------------------------------------------------------
+c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
+c-------------------------------------------------------------------
+
+c 1. NUAGES CONVECTIFS
+c
+      IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
+
+c Nuages diagnostiques pour Tiedtke
       CALL diagcld1(paprs,pplay,
      .             rain_con,snow_con,ibas_con,itop_con,
@@ -2132,7 +2234,29 @@
       ENDDO
       ENDDO
+
+      ELSE IF (iflag_cldcon.eq.3) THEN
+c  On prend pour les nuages convectifs le max du calcul de la
+c  convection et du calcul du pas de temps précédent diminué d'un facteur
+c  facttemps
+      facttemps=pdtphys/1.e4
+      do k=1,klev
+         do i=1,klon
+            rnebcon(i,k)=rnebcon(i,k)*facttemps
+            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
+     s      then
+                rnebcon(i,k)=rnebcon0(i,k)
+                clwcon(i,k)=clwcon0(i,k)
+            endif
+         enddo
+      enddo
+
+c   On prend la somme des fractions nuageuses et des contenus en eau
+      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
+      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
+
+
       ENDIF
 c
-c Nuages stratus artificiels:
+c 2. NUAGES STARTIFORMES
 c
       IF (ok_stratus) THEN
@@ -2174,4 +2298,5 @@
          ENDIF
          zx_rh(i,k) = q_seri(i,k)/zx_qs
+         zqsat(i,k)=zx_qs
       ENDDO
       ENDDO
@@ -2180,7 +2305,13 @@
 c parametres pour diagnostiques:
 c
+      if (ok_newmicro) then
+      CALL newmicro (paprs, pplay,ok_newmicro,
+     .            t_seri, cldliq, cldfra, cldtau, cldemi,
+     .            cldh, cldl, cldm, cldt, cldq)
+      else
       CALL nuage (paprs, pplay,
      .            t_seri, cldliq, cldfra, cldtau, cldemi,
      .            cldh, cldl, cldm, cldt, cldq)
+      endif
 c
 c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
@@ -3361,5 +3492,5 @@
      .      radsol,frugs,agesno,
      .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
-     .      t_ancien, q_ancien)
+     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
       ENDIF
 
