Index: /LMDZ.3.3/trunk/libf/phylmd/conemav.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/conemav.F	(revision 254)
+++ /LMDZ.3.3/trunk/libf/phylmd/conemav.F	(revision 254)
@@ -0,0 +1,147 @@
+      SUBROUTINE conemav (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,Ma,cape,tvp,iflag,
+     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
+ 
+c
+      IMPLICIT none
+c======================================================================
+c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
+c Objet: schema de convection de Emanuel (1991) interface
+c======================================================================
+c Arguments:
+c dtime--input-R-pas d'integration (s)
+c s-------input-R-la valeur "s" pour chaque couche
+c sigs----input-R-la valeur "sigma" de chaque couche
+c sig-----input-R-la valeur de "sigma" pour chaque niveau
+c psolpa--input-R-la pression au sol (en Pa)
+C pskapa--input-R-exponentiel kappa de psolpa
+c h-------input-R-enthalpie potentielle (Cp*T/P**kappa)
+c q-------input-R-vapeur d'eau (en kg/kg)
+c
+c work*: input et output: deux variables de travail,
+c                            on peut les mettre a 0 au debut
+c ALE-----input-R-energie disponible pour soulevement
+c
+C d_h-----output-R-increment de l'enthalpie potentielle (h)
+c d_q-----output-R-increment de la vapeur d'eau
+c rain----output-R-la pluie (mm/s)
+c snow----output-R-la neige (mm/s)
+c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
+c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
+c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
+c Cape----output-R-CAPE (J/kg)
+c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
+c                  adiabatiquement a partir du niveau 1 (K)
+c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
+c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
+c======================================================================
+c
+#include "dimensions.h"
+#include "dimphy.h"
+c
+      integer NTRAC
+      PARAMETER (NTRAC=nqmx-2)
+c
+       REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
+       REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
+       REAL tra(klon,klev,ntrac)
+       INTEGER ntra
+       REAL work1(klon,klev),work2(klon,klev)
+c
+       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
+       REAL d_tra(klon,klev,ntrac)
+       REAL rain(klon),snow(klon)
+c
+       INTEGER kbas(klon),ktop(klon)
+       REAL em_ph(klon,klev+1),em_p(klon,klev)
+       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
+       REAL Ma(klon,klev),cape(klon),tvp(klon,klev)
+       INTEGER iflag(klon)
+       REAL rflag(klon)
+       REAL pbase(klon),bbase(klon)
+       REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
+       REAL dplcldt(klon),dplcldr(klon)
+c
+       REAL zx_t,zdelta,zx_qs,zcor
+c
+       INTEGER noff, minorig
+       INTEGER i,k,itra
+       REAL qs(klon,klev)
+       REAL cbmf(klon)
+       SAVE cbmf
+       INTEGER ifrst
+       SAVE ifrst
+       DATA ifrst /0/
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+c
+c
+      IF (ifrst .EQ. 0) THEN
+         ifrst = 1
+         DO i = 1, klon
+          cbmf(i) = 0.
+         ENDDO
+      ENDIF
+
+      DO k = 1, klev+1
+         DO i=1,klon
+         em_ph(i,k) = paprs(i,k) / 100.0
+      ENDDO
+      ENDDO
+c
+      DO k = 1, klev
+         DO i=1,klon
+         em_p(i,k) = pplay(i,k) / 100.0
+      ENDDO
+      ENDDO
+
+c
+      DO k = 1, klev
+        DO i = 1, klon
+         zx_t = t(i,k)
+         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
+         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
+         zcor=1./(1.-retv*zx_qs)
+         qs(i,k)=zx_qs*zcor
+        ENDDO
+      ENDDO
+c
+      noff = 2
+      minorig = 2
+      CALL convect1(klon,klev,klev+1,noff,minorig,t,q,qs,u,v,
+     $              em_p,em_ph,iflag,
+     $              d_t,d_q,d_u,d_v,rain,cbmf,dtime,Ma)
+c
+      DO i = 1,klon
+        rain(i) = rain(i)/86400.
+        rflag(i)=iflag(i)
+      ENDDO
+c      call dump2d(iim,jjm-1,rflag(2:klon-1),'FLAG CONVECTION   ')
+c     if (klon.eq.1) then
+c        print*,'IFLAG ',iflag
+c     else
+c        write(*,'(96i1)') (iflag(i),i=2,klon-1)
+c     endif
+      DO k = 1, klev
+        DO i = 1, klon
+           d_t(i,k) = dtime*d_t(i,k)
+           d_q(i,k) = dtime*d_q(i,k)
+           d_u(i,k) = dtime*d_u(i,k)
+           d_v(i,k) = dtime*d_v(i,k)
+        ENDDO
+        DO itra = 1,ntra
+          DO i = 1, klon
+            d_tra(i,k,itra) = 0.
+          ENDDO
+        ENDDO
+      ENDDO
+ 
+c
+c
+c
+      RETURN
+      END
+ 
Index: /LMDZ.3.3/trunk/libf/phylmd/convect1.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/convect1.F	(revision 254)
+++ /LMDZ.3.3/trunk/libf/phylmd/convect1.F	(revision 254)
@@ -0,0 +1,645 @@
+      subroutine convect1(len,nd,ndp1,noff,minorig,
+     &                   t,q,qs,u,v,
+     &                   p,ph,iflag,ft,
+     &                   fq,fu,fv,precip,cbmf,delt,Ma)
+C.............................START PROLOGUE............................
+C
+C SCCS IDENTIFICATION:  @(#)convect1.f	1.1 04/21/00
+C                       19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v
+C
+C CONFIGURATION IDENTIFICATION:  None
+C
+C MODULE NAME:  convect1
+C
+C DESCRIPTION:
+C
+C convect1     The Emanuel Cumulus Convection Scheme
+C
+C CONTRACT NUMBER AND TITLE:  None
+C
+C REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
+C
+C CLASSIFICATION:  Unclassified
+C
+C RESTRICTIONS: None
+C
+C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
+C
+C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
+C                  Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
+C
+C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
+C
+C USAGE: call convect1(len,nd,noff,minorig,
+C    &                   t,q,qs,u,v,
+C    &                   p,ph,iflag,ft,
+C    &                   fq,fu,fv,precip,cbmf,delt)
+C
+C PARAMETERS:
+C      Name            Type         Usage            Description
+C   ----------      ----------     -------  ----------------------------
+C
+C      len           Integer        Input        first (i) dimension
+C      nd            Integer        Input        vertical (k) dimension
+C      ndp1          Integer        Input        nd + 1
+C      noff          Integer        Input        integer limit for convection (nd-noff)
+C      minorig       Integer        Input        First level of convection
+C      t             Real           Input        temperature
+C      q             Real           Input        specific hum
+C      qs            Real           Input        sat specific hum
+C      u             Real           Input        u-wind
+C      v             Real           Input        v-wind
+C      p             Real           Input        full level pressure
+C      ph            Real           Input        half level pressure
+C      iflag         Integer        Output       iflag on latitude strip
+C      ft            Real           Output       temp tend
+C      fq            Real           Output       spec hum tend
+C      fu            Real           Output       u-wind tend
+C      fv            Real           Output       v-wind tend
+C      cbmf          Real           In/Out       cumulus mass flux
+C      delt          Real           Input        time step
+C      iflag         Integer        Output       integer flag for Emanuel conditions
+C
+C COMMON BLOCKS:
+C      Block      Name     Type    Usage              Notes
+C     --------  --------   ----    ------   ------------------------
+C
+C FILES: None
+C
+C DATA BASES: None
+C
+C NON-FILE INPUT/OUTPUT: None
+C
+C ERROR CONDITIONS: None
+C
+C ADDITIONAL COMMENTS: None
+C
+C.................MAINTENANCE SECTION................................
+C
+C MODULES CALLED:
+C         Name           Description
+C         convect2        Emanuel cumulus convection tendency calculations
+C        -------     ----------------------
+C LOCAL VARIABLES AND
+C          STRUCTURES:
+C Name     Type    Description
+C -------  ------  -----------
+C See Comments Below
+C
+C i        Integer loop index
+C k        Integer loop index
+c
+C METHOD:
+C
+C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
+C       convective scheme for use in climate models.
+C
+C FILES: None
+C
+C INCLUDE FILES: None
+C
+C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
+C
+C..............................END PROLOGUE.............................
+c
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphy.h"
+c
+      integer len
+      integer nd
+      integer ndp1
+      integer noff
+      real t(len,nd)
+      real q(len,nd)
+      real qs(len,nd)
+      real u(len,nd)
+      real v(len,nd)
+      real p(len,nd)
+      real ph(len,ndp1)
+      integer iflag(len)
+      real ft(len,nd)
+      real fq(len,nd)
+      real fu(len,nd)
+      real fv(len,nd)
+      real precip(len)
+      real cbmf(len)
+      real Ma(len,nd)
+      integer minorig
+      real delt,cpd,cpv,cl,rv,rd,lv0,g
+      real sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp
+      real alpha,entp,coeffs,coeffr,omtrain,cu
+c
+!-------------------------------------------------------------------
+! --- ARGUMENTS
+!-------------------------------------------------------------------
+! --- On input:
+!
+!  t:   Array of absolute temperature (K) of dimension ND, with first
+!       index corresponding to lowest model level. Note that this array
+!       will be altered by the subroutine if dry convective adjustment
+!       occurs and if IPBL is not equal to 0.
+!
+!  q:   Array of specific humidity (gm/gm) of dimension ND, with first
+!       index corresponding to lowest model level. Must be defined
+!       at same grid levels as T. Note that this array will be altered
+!       if dry convective adjustment occurs and if IPBL is not equal to 0.
+!
+!  qs:  Array of saturation specific humidity of dimension ND, with first
+!       index corresponding to lowest model level. Must be defined
+!       at same grid levels as T. Note that this array will be altered
+!       if dry convective adjustment occurs and if IPBL is not equal to 0.
+!
+!  u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
+!       index corresponding with the lowest model level. Defined at
+!       same levels as T. Note that this array will be altered if
+!       dry convective adjustment occurs and if IPBL is not equal to 0.
+!
+!  v:   Same as u but for meridional velocity.
+!
+!  tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
+!       where NTRA is the number of different tracers. If no
+!       convective tracer transport is needed, define a dummy
+!       input array of dimension (ND,1). Tracers are defined at
+!       same vertical levels as T. Note that this array will be altered
+!       if dry convective adjustment occurs and if IPBL is not equal to 0.
+!
+!  p:   Array of pressure (mb) of dimension ND, with first
+!       index corresponding to lowest model level. Must be defined
+!       at same grid levels as T.
+!
+!  ph:  Array of pressure (mb) of dimension ND+1, with first index
+!       corresponding to lowest level. These pressures are defined at
+!       levels intermediate between those of P, T, Q and QS. The first
+!       value of PH should be greater than (i.e. at a lower level than)
+!       the first value of the array P.
+!
+!  nl:  The maximum number of levels to which convection can penetrate, plus 1.
+!       NL MUST be less than or equal to ND-1.
+!
+!  delt: The model time step (sec) between calls to CONVECT
+!
+!----------------------------------------------------------------------------
+! ---   On Output:
+!
+!  iflag: An output integer whose value denotes the following:
+!       VALUE   INTERPRETATION
+!       -----   --------------
+!         0     Moist convection occurs.
+!         1     Moist convection occurs, but a CFL condition
+!               on the subsidence warming is violated. This
+!               does not cause the scheme to terminate.
+!         2     Moist convection, but no precip because ep(inb) lt 0.0001
+!         3     No moist convection because new cbmf is 0 and old cbmf is 0.
+!         4     No moist convection; atmosphere is not
+!               unstable
+!         6     No moist convection because ihmin le minorig.
+!         7     No moist convection because unreasonable
+!               parcel level temperature or specific humidity.
+!         8     No moist convection: lifted condensation
+!               level is above the 200 mb level.
+!         9     No moist convection: cloud base is higher
+!               then the level NL-1.
+!
+!  ft:   Array of temperature tendency (K/s) of dimension ND, defined at same
+!        grid levels as T, Q, QS and P.
+!
+!  fq:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
+!        defined at same grid levels as T, Q, QS and P.
+!
+!  fu:   Array of forcing of zonal velocity (m/s^2) of dimension ND,
+!        defined at same grid levels as T.
+!
+!  fv:   Same as FU, but for forcing of meridional velocity.
+!
+!  ftra: Array of forcing of tracer content, in tracer mixing ratio per
+!        second, defined at same levels as T. Dimensioned (ND,NTRA).
+!
+!  precip: Scalar convective precipitation rate (mm/day).
+!
+!  wd:   A convective downdraft velocity scale. For use in surface
+!        flux parameterizations. See convect.ps file for details.
+!
+!  tprime: A convective downdraft temperature perturbation scale (K).
+!          For use in surface flux parameterizations. See convect.ps
+!          file for details.
+!
+!  qprime: A convective downdraft specific humidity
+!          perturbation scale (gm/gm).
+!          For use in surface flux parameterizations. See convect.ps
+!          file for details.
+!
+!  cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
+!        BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
+!        ITS NEXT CALL. That is, the value of CBMF must be "remembered"
+!        by the calling program between calls to CONVECT.
+!
+!  det:   Array of detrainment mass flux of dimension ND.
+!
+!-------------------------------------------------------------------
+c
+c  Local arrays
+c
+      integer nl
+      integer nlp
+      integer nlm
+      integer i,k,n
+      real delti
+      real rowl
+      real clmcpv
+      real clmcpd
+      real cpdmcp
+      real cpvmcpd
+      real eps
+      real epsi
+      real epsim1
+      real ginv
+      real hrd
+      real prccon1
+      integer icbmax
+      real lv(klon,klev)
+      real cpn(klon,klev)
+      real cpx(klon,klev)
+      real tv(klon,klev)
+      real gz(klon,klev)
+      real hm(klon,klev)
+      real h(klon,klev)
+      real work(klon)
+      integer ihmin(klon)
+      integer nk(klon)
+      real rh(klon)
+      real chi(klon)
+      real plcl(klon)
+      integer icb(klon)
+      real tnk(klon)
+      real qnk(klon)
+      real gznk(klon)
+      real pnk(klon)
+      real qsnk(klon)
+      real ticb(klon)
+      real gzicb(klon)
+      real tp(klon,klev)
+      real tvp(klon,klev)
+      real clw(klon,klev)
+c
+      real ah0(klon),cpp(klon)
+      real tg,qg,s,alv,tc,ahg,denom,es,rg
+c
+      integer ncum
+      integer idcum(klon)
+c
+      cpd=1005.7
+      cpv=1870.0
+      cl=4190.0
+      rv=461.5
+      rd=287.04
+      lv0=2.501E6
+      g=9.8
+C
+C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
+C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
+C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
+C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
+C   ***               BETWEEN 0 C AND TLCRIT)                        ***
+C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
+C   ***                       FORMULATION                            ***
+C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
+C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
+C   ***                        OF CLOUD                              ***
+C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
+C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
+C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
+C   ***                          OF RAIN                             ***
+C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
+C   ***                          OF SNOW                             ***
+C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
+C   ***                         TRANSPORT                            ***
+C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
+C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
+C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
+C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
+C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
+C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
+c
+      sigs=0.12
+      sigd=0.05
+      elcrit=0.0011
+      tlcrit=-55.0
+      omtsnow=5.5
+      dtmax=0.9
+      damp=0.1
+      alpha=0.2
+      entp=1.5
+      coeffs=0.8
+      coeffr=1.0
+      omtrain=50.0
+c
+      cu=0.70
+      damp=0.1
+c
+c
+c Define nl, nlp, nlm, and delti
+c
+      nl=nd-noff
+      nlp=nl+1
+      nlm=nl-1
+      delti=1.0/delt
+!
+!-------------------------------------------------------------------
+! --- SET CONSTANTS
+!-------------------------------------------------------------------
+!
+      rowl=1000.0
+      clmcpv=cl-cpv
+      clmcpd=cl-cpd
+      cpdmcp=cpd-cpv
+      cpvmcpd=cpv-cpd
+      eps=rd/rv
+      epsi=1.0/eps
+      epsim1=epsi-1.0
+      ginv=1.0/g
+      hrd=0.5*rd
+      prccon1=86400.0*1000.0/(rowl*g)
+!
+! dtmax is the maximum negative temperature perturbation.
+!
+!=====================================================================
+! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
+!=====================================================================
+!
+      do 20 k=1,nd
+        do 10 i=1,len
+         ft(i,k)=0.0
+         fq(i,k)=0.0
+         fu(i,k)=0.0
+         fv(i,k)=0.0
+         tvp(i,k)=0.0
+         tp(i,k)=0.0
+         clw(i,k)=0.0
+         gz(i,k) = 0.
+ 10     continue
+ 20   continue
+      do 60 i=1,len
+        precip(i)=0.0
+        iflag(i)=0
+ 60   continue
+c
+!=====================================================================
+! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
+!=====================================================================
+      do 110 k=1,nl+1
+        do 100 i=1,len
+          lv(i,k)= lv0-clmcpv*(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)
+          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
+ 100    continue
+ 110  continue
+c
+c gz = phi at the full levels (same as p).
+c
+      do 120 i=1,len
+        gz(i,1)=0.0
+ 120  continue
+      do 140 k=2,nlp
+        do 130 i=1,len
+          gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
+     &         *(p(i,k-1)-p(i,k))/ph(i,k)
+ 130    continue
+ 140  continue
+c
+c h  = phi + cpT (dry static energy).
+c hm = phi + cp(T-Tbase)+Lq
+c
+      do 170 k=1,nlp
+        do 160 i=1,len
+          h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
+          hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k)
+ 160    continue
+ 170  continue
+c
+!-------------------------------------------------------------------
+! --- Find level of minimum moist static energy
+! --- If level of minimum moist static energy coincides with
+! --- or is lower than minimum allowable parcel origin level,
+! --- set iflag to 6.
+!-------------------------------------------------------------------
+      do 180 i=1,len
+       work(i)=1.0e12
+       ihmin(i)=nl
+ 180  continue
+      do 200 k=2,nlp
+        do 190 i=1,len
+         if((hm(i,k).lt.work(i)).and.
+     &      (hm(i,k).lt.hm(i,k-1)))then
+           work(i)=hm(i,k)
+           ihmin(i)=k
+         endif
+ 190    continue
+ 200  continue
+      do 210 i=1,len
+        ihmin(i)=min(ihmin(i),nlm)
+        if(ihmin(i).le.minorig)then
+          iflag(i)=6
+        endif
+ 210  continue
+c
+!-------------------------------------------------------------------
+! --- Find that model level below the level of minimum moist static
+! --- energy that has the maximum value of moist static energy
+!-------------------------------------------------------------------
+ 
+      do 220 i=1,len
+       work(i)=hm(i,minorig)
+       nk(i)=minorig
+ 220  continue
+      do 240 k=minorig+1,nl
+        do 230 i=1,len
+         if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
+           work(i)=hm(i,k)
+           nk(i)=k
+         endif
+ 230     continue
+ 240  continue
+!-------------------------------------------------------------------
+! --- Check whether parcel level temperature and specific humidity
+! --- are reasonable
+!-------------------------------------------------------------------
+       do 250 i=1,len
+       if(((t(i,nk(i)).lt.250.0).or.
+     &      (q(i,nk(i)).le.0.0).or.
+     &      (p(i,ihmin(i)).lt.400.0)).and.
+     &      (iflag(i).eq.0))iflag(i)=7
+ 250   continue
+!-------------------------------------------------------------------
+! --- Calculate lifted condensation level of air at parcel origin level
+! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
+!-------------------------------------------------------------------
+       do 260 i=1,len
+        tnk(i)=t(i,nk(i))
+        qnk(i)=q(i,nk(i))
+        gznk(i)=gz(i,nk(i))
+        pnk(i)=p(i,nk(i))
+        qsnk(i)=qs(i,nk(i))
+c
+        rh(i)=qnk(i)/qsnk(i)
+        rh(i)=min(1.0,rh(i))
+        chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
+        plcl(i)=pnk(i)*(rh(i)**chi(i))
+        if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
+     &   .and.(iflag(i).eq.0))iflag(i)=8
+ 260   continue
+!-------------------------------------------------------------------
+! --- Calculate first level above lcl (=icb)
+!-------------------------------------------------------------------
+      do 270 i=1,len
+       icb(i)=nlm
+ 270  continue
+c
+      do 290 k=minorig,nl
+        do 280 i=1,len
+          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
+     &    icb(i)=min(icb(i),k)
+ 280    continue
+ 290  continue
+c
+      do 300 i=1,len
+        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
+ 300  continue
+c
+c Compute icbmax.
+c
+      icbmax=2
+      do 310 i=1,len
+        icbmax=max(icbmax,icb(i))
+ 310  continue
+!
+!-------------------------------------------------------------------
+! --- Calculates the lifted parcel virtual temperature at nk,
+! --- the actual temperature, and the adiabatic
+! --- liquid water content. The procedure is to solve the equation.
+!     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
+!-------------------------------------------------------------------
+!
+      do 320 i=1,len
+        tnk(i)=t(i,nk(i))
+        qnk(i)=q(i,nk(i))
+        gznk(i)=gz(i,nk(i))
+        ticb(i)=t(i,icb(i))
+        gzicb(i)=gz(i,icb(i))
+ 320  continue
+c
+c   ***  Calculate certain parcel quantities, including static energy   ***
+c
+      do 330 i=1,len
+        ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
+     &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
+        cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
+ 330  continue
+c
+c   ***   Calculate lifted parcel quantities below cloud base   ***
+c
+        do 350 k=minorig,icbmax-1
+          do 340 i=1,len
+           tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i)
+           tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi)
+  340     continue
+  350   continue
+c
+c    ***  Find lifted parcel quantities above cloud base    ***
+c
+        do 360 i=1,len
+         tg=ticb(i)
+         qg=qs(i,icb(i))
+         alv=lv0-clmcpv*(ticb(i)-273.15)
+c
+c First iteration.
+c
+          s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i))
+          s=1./s
+          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+          tg=tg+s*(ah0(i)-ahg)
+          tg=max(tg,35.0)
+          tc=tg-273.15
+          denom=243.5+tc
+          if(tc.ge.0.0)then
+           es=6.112*exp(17.67*tc/denom)
+          else
+           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+          endif
+          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+c
+c Second iteration.
+c
+          s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i))
+          s=1./s
+          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+          tg=tg+s*(ah0(i)-ahg)
+          tg=max(tg,35.0)
+          tc=tg-273.15
+          denom=243.5+tc
+          if(tc.ge.0.0)then
+           es=6.112*exp(17.67*tc/denom)
+          else
+           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+          end if
+          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+c
+         alv=lv0-clmcpv*(ticb(i)-273.15)
+         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
+     &   -gz(i,icb(i))-alv*qg)/cpd
+         clw(i,icb(i))=qnk(i)-qg
+         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
+         rg=qg/(1.-qnk(i))
+         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
+  360   continue
+c
+      do 380 k=minorig,icbmax
+       do 370 i=1,len
+         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
+ 370   continue
+ 380  continue
+c
+!-------------------------------------------------------------------
+! --- Test for instability.
+! --- If there was no convection at last time step and parcel
+! --- is stable at icb, then set iflag to 4.
+!-------------------------------------------------------------------
+ 
+      do 390 i=1,len
+        if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and.
+     &  (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4
+ 390  continue
+ 
+!=====================================================================
+! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
+!=====================================================================
+c
+      ncum=0
+      do 400 i=1,len
+        if(iflag(i).eq.0)then
+           ncum=ncum+1
+           idcum(ncum)=i
+        endif
+ 400  continue
+c
+c Call convect2, which compresses the points and computes the heating,
+c moistening, velocity mixing, and precipiation.
+c
+c     print*,'cpd avant convect2 ',cpd
+      if(ncum.gt.0)then
+      call convect2(ncum,idcum,len,nd,ndp1,nl,minorig,
+     &              nk,icb,
+     &              t,q,qs,u,v,gz,tv,tp,tvp,clw,h,
+     &              lv,cpn,p,ph,ft,fq,fu,fv,
+     &              tnk,qnk,gznk,plcl,
+     &              precip,cbmf,iflag,
+     &              delt,cpd,cpv,cl,rv,rd,lv0,g,
+     &              sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
+     &              alpha,entp,coeffs,coeffr,omtrain,cu,Ma)
+      endif
+c
+      return
+      end
Index: /LMDZ.3.3/trunk/libf/phylmd/convect2.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/convect2.F	(revision 254)
+++ /LMDZ.3.3/trunk/libf/phylmd/convect2.F	(revision 254)
@@ -0,0 +1,1391 @@
+      subroutine convect2(ncum,idcum,len,nd,ndp1,nl,minorig,
+     &                 nk1,icb1,
+     &                 t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
+     &                 lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
+     &                 tnk1,qnk1,gznk1,plcl1,
+     &                 precip1,cbmf1,iflag1,
+     &                 delt,cpd,cpv,cl,rv,rd,lv0,g,
+     &                 sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
+     &                 alpha,entp,coeffs,coeffr,omtrain,cu,Ma)
+C.............................START PROLOGUE............................
+C
+C SCCS IDENTIFICATION:  @(#)convect2.f	1.2 05/18/00
+C                       22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v
+C
+C CONFIGURATION IDENTIFICATION:  None
+C
+C MODULE NAME:  convect2
+C
+C DESCRIPTION:
+C
+C convect1     The Emanuel Cumulus Convection Scheme - compute tendencies
+C
+C CONTRACT NUMBER AND TITLE:  None
+C
+C REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
+C
+C CLASSIFICATION:  Unclassified
+C
+C RESTRICTIONS: None
+C
+C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
+C
+C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
+C                  Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
+C
+C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
+C
+C USAGE: call convect2(ncum,idcum,len,nd,nl,minorig,
+C    &                 nk1,icb1,
+C    &                 t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
+C    &                 lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
+C    &                 tnk1,qnk1,gznk1,plcl1,
+C    &                 precip1,cbmf1,iflag1,
+C    &                 delt,cpd,cpv,cl,rv,rd,lv0,g,
+C    &                 sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
+C    &                 alpha,entp,coeffs,coeffr,omtrain,cu)
+C
+C PARAMETERS:
+C      Name            Type         Usage            Description
+C   ----------      ----------     -------  ----------------------------
+C
+C      ncum          Integer        Input        number of cumulus points
+C      idcum         Integer        Input        index of cumulus point
+C      len           Integer        Input        first dimension
+C      nd            Integer        Input        total vertical dimension
+C      ndp1          Integer        Input        nd + 1
+C      nl            Integer        Input        vertical dimension for convection
+C      minorig       Integer        Input        First level where convection is allow to begin
+C      nk1           Integer        Input        First level of convection
+C      ncb1          Integer        Input        Level of free convection
+C      t1            Real           Input        temperature
+C      q1            Real           Input        specific hum
+C      qs1           Real           Input        sat specific hum
+C      u1            Real           Input        u-wind
+C      v1            Real           Input        v-wind
+C      gz1           Real           Inout        geop
+C      tv1           Real           Input        virtual temp
+C      tp1           Real           Input
+C      clw1          Real           Inout        cloud liquid water
+C      h1            Real           Inout
+C      lv1           Real           Inout
+C      cpn1          Real           Inout
+C      p1            Real           Input        full level pressure
+C      ph1           Real           Input        half level pressure
+C      ft1           Real           Output       temp tend
+C      fq1           Real           Output       spec hum tend
+C      fu1           Real           Output       u-wind tend
+C      fv1           Real           Output       v-wind tend
+C      precip1       Real           Output       prec
+C      cbmf1         Real           In/Out       cumulus mass flux
+C      iflag1        Integer        Output       iflag on latitude strip
+C      delt          Real           Input        time step
+C      cpd           Integer        Input        See description below
+C      cpv           Integer        Input        See description below
+C      cl            Integer        Input        See description below
+C      rv            Integer        Input        See description below
+C      rd            Integer        Input        See description below
+C      lv0           Integer        Input        See description below
+C      g             Integer        Input        See description below
+C      sigs          Integer        Input        See description below
+C      sigd          Integer        Input        See description below
+C      elcrit        Integer        Input        See description below
+C      tlcrit        Integer        Input        See description below
+C      omtsnow       Integer        Input        See description below
+C      dtmax         Integer        Input        See description below
+C      damp          Integer        Input        See description below
+C      alpha         Integer        Input        See description below
+C      ent           Integer        Input        See description below
+C      coeffs        Integer        Input        See description below
+C      coeffr        Integer        Input        See description below
+C      omtrain       Integer        Input        See description below
+C      cu            Integer        Input        See description below
+C
+C COMMON BLOCKS:
+C      Block      Name     Type    Usage              Notes
+C     --------  --------   ----    ------   ------------------------
+C
+C FILES: None
+C
+C DATA BASES: None
+C
+C NON-FILE INPUT/OUTPUT: None
+C
+C ERROR CONDITIONS: None
+C
+C ADDITIONAL COMMENTS: None
+C
+C.................MAINTENANCE SECTION................................
+C
+C MODULES CALLED:
+C         Name           Description
+C         zilch        Zero out an array
+C        -------     ----------------------
+C LOCAL VARIABLES AND
+C          STRUCTURES:
+C Name     Type    Description
+C -------  ------  -----------
+C See Comments Below
+C
+C i        Integer loop index
+C k        Integer loop index
+c
+C METHOD:
+C
+C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
+C       convective scheme for use in climate models.
+C
+C FILES: None
+C
+C INCLUDE FILES: None
+C
+C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
+C
+C..............................END PROLOGUE.............................
+c
+c
+      implicit none
+c
+#include "dimensions.h"
+#include "dimphy.h"
+c
+      integer kmax2,imax2,kmin2,imin2
+      real ftmax2,ftmin2
+      integer kmax,imax,kmin,imin
+      real ftmax,ftmin
+c
+      integer ncum
+      integer idcum(len)
+      integer len
+      integer nd
+      integer ndp1
+      integer nl
+      integer minorig
+      integer nk1(len)
+      integer icb1(len)
+      real t1(len,nd)
+      real q1(len,nd)
+      real qs1(len,nd)
+      real u1(len,nd)
+      real v1(len,nd)
+      real gz1(len,nd)
+      real tv1(len,nd)
+      real tp1(len,nd)
+      real tvp1(len,nd)
+      real clw1(len,nd)
+      real h1(len,nd)
+      real lv1(len,nd)
+      real cpn1(len,nd)
+      real p1(len,nd)
+      real ph1(len,ndp1)
+      real ft1(len,nd)
+      real fq1(len,nd)
+      real fu1(len,nd)
+      real fv1(len,nd)
+      real tnk1(len)
+      real qnk1(len)
+      real gznk1(len)
+      real precip1(len)
+      real cbmf1(len)
+      real plcl1(len)
+      integer iflag1(len)
+      real delt
+      real cpd
+      real cpv
+      real cl
+      real rv
+      real rd
+      real lv0
+      real g
+      real sigs    ! SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE
+      real sigd    ! SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT
+      real elcrit  ! ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm)
+      real tlcrit  ! TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-
+c                     CONVERSION THRESHOLD IS ASSUMED TO BE ZERO
+      real omtsnow ! OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW
+      real dtmax   ! DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION
+c                    A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC.
+      real damp
+      real alpha
+      real entp    ! ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT FORMULATION
+      real coeffs  ! COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF SNOW
+      real coeffr  ! COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF RAIN
+      real omtrain ! OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN
+      real cu      ! CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM TRANSPORT
+c
+      real Ma(len,nd)
+c
+C
+C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
+C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
+C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
+C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
+C   ***               BETWEEN 0 C AND TLCRIT)                        ***
+C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
+C   ***                       FORMULATION                            ***
+C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
+C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
+C   ***                        OF CLOUD                              ***
+C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
+C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
+C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
+C   ***                          OF RAIN                             ***
+C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
+C   ***                          OF SNOW                             ***
+C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
+C   ***                         TRANSPORT                            ***
+C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
+C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
+C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
+C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
+C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
+C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
+c
+c Local arrays.
+c
+      real work(ncum)
+      real t(ncum,klev)
+      real q(ncum,klev)
+      real qs(ncum,klev)
+      real u(ncum,klev)
+      real v(ncum,klev)
+      real gz(ncum,klev)
+      real h(ncum,klev)
+      real lv(ncum,klev)
+      real cpn(ncum,klev)
+      real p(ncum,klev)
+      real ph(ncum,klev)
+      real ft(ncum,klev)
+      real fq(ncum,klev)
+      real fu(ncum,klev)
+      real fv(ncum,klev)
+      real precip(ncum)
+      real cbmf(ncum)
+      real plcl(ncum)
+      real tnk(ncum)
+      real qnk(ncum)
+      real gznk(ncum)
+      real tv(ncum,klev)
+      real tp(ncum,klev)
+      real tvp(ncum,klev)
+      real clw(ncum,klev)
+c     real det(ncum,klev)
+      real dph(ncum,klev)
+c     real wd(ncum)
+c     real tprime(ncum)
+c     real qprime(ncum)
+      real ah0(ncum)
+      real ep(ncum,klev)
+      real sigp(ncum,klev)
+      integer nent(ncum,klev)
+      real water(ncum,klev)
+      real evap(ncum,klev)
+      real mp(ncum,klev)
+      real m(ncum,klev)
+      real qti
+      real wt(ncum,klev)
+      real hp(ncum,klev)
+      real lvcp(ncum,klev)
+      real elij(ncum,klev,klev)
+      real ment(ncum,klev,klev)
+      real sij(ncum,klev,klev)
+      real qent(ncum,klev,klev)
+      real uent(ncum,klev,klev)
+      real vent(ncum,klev,klev)
+      real qp(ncum,klev)
+      real up(ncum,klev)
+      real vp(ncum,klev)
+      real cape(ncum)
+      real capem(ncum)
+      real frac(ncum)
+      real dtpbl(ncum)
+      real tvpplcl(ncum)
+      real tvaplcl(ncum)
+      real dtmin(ncum)
+      real w3d(ncum,klev)
+      real am(ncum)
+      real ents(ncum)
+      real uav(ncum)
+      real vav(ncum)
+c
+      integer iflag(ncum)
+      integer nk(ncum)
+      integer icb(ncum)
+      integer inb(ncum)
+      integer inb1(ncum)
+      integer jtt(ncum)
+c
+      integer nn,i,k,n,icbmax,nlp,j
+      integer ij
+      integer nn2,nn3
+      real clmcpv
+      real clmcpd
+      real cpdmcp
+      real cpvmcpd
+      real eps
+      real epsi
+      real epsim1
+      real tg,qg,s,alv,tc,ahg,denom,es,rg,ginv,rowl
+      real delti
+      real tca,elacrit
+      real by,defrac
+c     real byp
+      real byp(ncum)
+      logical lcape(ncum)
+      real dbo
+      real bf2,anum,dei,altem,cwat,stemp
+      real alt,qp1,smid,sjmax,sjmin
+      real delp,delm
+      real awat,coeff,afac,revap,dhdp,fac,qstm,rat
+      real qsm,sigt,b6,c6
+      real dpinv,cpinv
+      real fqold,ftold,fuold,fvold
+      real wdtrain(ncum),xxx
+      real bsum(ncum,klev)
+      real asij(ncum)
+      real smin(ncum)
+      real scrit(ncum)
+c     real amp1,ad
+      real amp1(ncum),ad(ncum)
+      logical lwork(ncum)
+      integer num1,num2
+c
+c     print*,'cpd en entree de convect2 ',cpd
+      nlp=nl+1
+c
+      rowl=1000.0
+      ginv=1.0/g
+      delti=1.0/delt
+c
+c Define some thermodynamic variables.
+c
+      clmcpv=cl-cpv
+      clmcpd=cl-cpd
+      cpdmcp=cpd-cpv
+      cpvmcpd=cpv-cpd
+      eps=rd/rv
+      epsi=1.0/eps
+      epsim1=epsi-1.0
+c
+c Compress the fields.
+c
+      do 110 k=1,nl+1
+       nn=0
+	do 100 i=1,len
+	  if(iflag1(i).eq.0)then
+	    nn=nn+1
+	    t(nn,k)=t1(i,k)
+	    q(nn,k)=q1(i,k)
+	    qs(nn,k)=qs1(i,k)
+	    u(nn,k)=u1(i,k)
+	    v(nn,k)=v1(i,k)
+	    gz(nn,k)=gz1(i,k)
+	    h(nn,k)=h1(i,k)
+	    lv(nn,k)=lv1(i,k)
+	    cpn(nn,k)=cpn1(i,k)
+	    p(nn,k)=p1(i,k)
+	    ph(nn,k)=ph1(i,k)
+	    tv(nn,k)=tv1(i,k)
+	    tp(nn,k)=tp1(i,k)
+	    tvp(nn,k)=tvp1(i,k)
+	    clw(nn,k)=clw1(i,k)
+	  endif
+ 100    continue
+c       print*,'100 ncum,nn',ncum,nn
+ 110  continue
+      nn=0
+      do 150 i=1,len
+	if(iflag1(i).eq.0)then
+	  nn=nn+1
+	  cbmf(nn)=cbmf1(i)
+	  plcl(nn)=plcl1(i)
+	  tnk(nn)=tnk1(i)
+	  qnk(nn)=qnk1(i)
+	  gznk(nn)=gznk1(i)
+	  nk(nn)=nk1(i)
+	  icb(nn)=icb1(i)
+	  iflag(nn)=iflag1(i)
+	endif
+ 150  continue
+c       print*,'150 ncum,nn',ncum,nn
+c
+c Initialize the tendencies, det, wd, tprime, qprime.
+c
+      do 170 k=1,nl
+	do 160 i=1,ncum
+c         det(i,k)=0.0
+	  ft(i,k)=0.0
+	  fu(i,k)=0.0
+	  fv(i,k)=0.0
+	  fq(i,k)=0.0
+	  dph(i,k)=ph(i,k)-ph(i,k+1)
+	  ep(i,k)=0.0
+	  sigp(i,k)=sigs
+ 160    continue
+ 170  continue
+      do 180 i=1,ncum
+c      wd(i)=0.0
+c      tprime(i)=0.0
+c      qprime(i)=0.0
+       precip(i)=0.0
+       ft(i,nl+1)=0.0
+       fu(i,nl+1)=0.0
+       fv(i,nl+1)=0.0
+       fq(i,nl+1)=0.0
+ 180  continue
+c
+c Compute icbmax.
+c
+      icbmax=2
+      do 230 i=1,ncum
+	icbmax=max(icbmax,icb(i))
+ 230  continue
+c
+c
+!=====================================================================
+! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
+!=====================================================================
+c
+c ---       The procedure is to solve the equation.
+c              cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
+c
+c   ***  Calculate certain parcel quantities, including static energy   ***
+c
+c
+      do 240 i=1,ncum
+	ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
+     &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
+ 240  continue
+c
+c
+c    ***  Find lifted parcel quantities above cloud base    ***
+c
+c
+	do 300 k=minorig+1,nl
+	  do 290 i=1,ncum
+	    if(k.ge.(icb(i)+1))then
+	      tg=t(i,k)
+	      qg=qs(i,k)
+	      alv=lv0-clmcpv*(t(i,k)-273.15)
+c
+c First iteration.
+c
+	       s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k))
+	       s=1./s
+	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
+	       tg=tg+s*(ah0(i)-ahg)
+	       tg=max(tg,35.0)
+	       tc=tg-273.15
+	       denom=243.5+tc
+	       if(tc.ge.0.0)then
+		es=6.112*exp(17.67*tc/denom)
+	       else
+		es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+	       endif
+		qg=eps*es/(p(i,k)-es*(1.-eps))
+c
+c Second iteration.
+c
+	       s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k))
+	       s=1./s
+	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
+	       tg=tg+s*(ah0(i)-ahg)
+	       tg=max(tg,35.0)
+	       tc=tg-273.15
+	       denom=243.5+tc
+	       if(tc.ge.0.0)then
+		es=6.112*exp(17.67*tc/denom)
+	       else
+		es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+	       endif
+		qg=eps*es/(p(i,k)-es*(1.-eps))
+c
+	       alv=lv0-clmcpv*(t(i,k)-273.15)
+c      print*,'cpd dans convect2 ',cpd
+c      print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
+c      print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
+	       tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)
+     &                  /cpd
+c              if (.not.cpd.gt.1000.) then
+c                  print*,'CPD=',cpd
+c                  stop
+c              endif
+               clw(i,k)=qnk(i)-qg
+               clw(i,k)=max(0.0,clw(i,k))
+               rg=qg/(1.-qnk(i))
+               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
+            endif
+  290     continue
+  300   continue
+c
+!=====================================================================
+! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
+! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
+! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
+!=====================================================================
+c
+      do 320 k=minorig+1,nl
+        do 310 i=1,ncum
+          if(k.ge.(nk(i)+1))then
+            tca=tp(i,k)-273.15
+            if(tca.ge.0.0)then
+              elacrit=elcrit
+            else
+              elacrit=elcrit*(1.0-tca/tlcrit)
+            endif
+            elacrit=max(elacrit,0.0)
+            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
+            ep(i,k)=max(ep(i,k),0.0 )
+            ep(i,k)=min(ep(i,k),1.0 )
+            sigp(i,k)=sigs
+          endif
+ 310    continue
+ 320  continue
+c
+!=====================================================================
+! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
+! --- VIRTUAL TEMPERATURE
+!=====================================================================
+c
+      do 340 k=minorig+1,nl
+        do 330 i=1,ncum
+        if(k.ge.(icb(i)+1))then
+          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
+c         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
+c         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
+        endif
+ 330    continue
+ 340  continue
+      do 350 i=1,ncum
+       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
+ 350  continue
+c
+c
+c=====================================================================
+c --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
+c=====================================================================
+c
+        do 360 i=1,ncum*nlp
+          nent(i,1)=0
+          water(i,1)=0.0
+          evap(i,1)=0.0
+          mp(i,1)=0.0
+          m(i,1)=0.0
+          wt(i,1)=omtsnow
+          hp(i,1)=h(i,1)
+c         if(.not.cpn(i,1).gt.900.) then
+c         print*,'i,lv(i,1),cpn(i,1)'
+c         print*, i,lv(i,1),cpn(i,1)
+c         k=(i-1)/ncum+1
+c         print*,'i,k',mod(i,ncum),k,'  cpn',cpn(mod(i,ncum),k)
+c         stop
+c         endif
+          lvcp(i,1)=lv(i,1)/cpn(i,1)
+ 360    continue
+c
+      do 380 i=1,ncum*nlp*nlp
+        elij(i,1,1)=0.0
+        ment(i,1,1)=0.0
+        sij(i,1,1)=0.0
+ 380  continue
+c
+      do 400 k=1,nlp
+       do 390 j=1,nlp
+          do 385 i=1,ncum
+            qent(i,k,j)=q(i,j)
+            uent(i,k,j)=u(i,j)
+            vent(i,k,j)=v(i,j)
+ 385      continue
+ 390    continue
+ 400  continue
+c
+      do 420 i=1,ncum
+        qp(i,1)=q(i,1)
+        up(i,1)=u(i,1)
+        vp(i,1)=v(i,1)
+ 420  continue
+      do 440 k=2,nlp
+        do 430 i=1,ncum
+          qp(i,k)=q(i,k-1)
+          up(i,k)=u(i,k-1)
+          vp(i,k)=v(i,k-1)
+ 430    continue
+ 440  continue
+c
+c=====================================================================
+c  --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
+c  --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
+c  --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
+c=====================================================================
+c
+      do 510 i=1,ncum
+        cape(i)=0.0
+        capem(i)=0.0
+        inb(i)=icb(i)+1
+        inb1(i)=inb(i)
+ 510  continue
+c
+c Originial Code
+c
+c     do 530 k=minorig+1,nl-1
+c       do 520 i=1,ncum
+c         if(k.ge.(icb(i)+1))then
+c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+c           byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+c           cape(i)=cape(i)+by
+c           if(by.ge.0.0)inb1(i)=k+1
+c           if(cape(i).gt.0.0)then
+c             inb(i)=k+1
+c             capem(i)=cape(i)
+c           endif
+c         endif
+c520    continue
+c530  continue
+c     do 540 i=1,ncum
+c         byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
+c         cape(i)=capem(i)+byp
+c         defrac=capem(i)-cape(i)
+c         defrac=max(defrac,0.001)
+c         frac(i)=-cape(i)/defrac
+c         frac(i)=min(frac(i),1.0)
+c         frac(i)=max(frac(i),0.0)
+c540   continue
+c
+c K Emanuel fix
+c
+c     call zilch(byp,ncum)
+c     do 530 k=minorig+1,nl-1
+c       do 520 i=1,ncum
+c         if(k.ge.(icb(i)+1))then
+c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+c           cape(i)=cape(i)+by
+c           if(by.ge.0.0)inb1(i)=k+1
+c           if(cape(i).gt.0.0)then
+c             inb(i)=k+1
+c             capem(i)=cape(i)
+c             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+c           endif
+c         endif
+c520    continue
+c530  continue
+c     do 540 i=1,ncum
+c         inb(i)=max(inb(i),inb1(i))
+c         cape(i)=capem(i)+byp(i)
+c         defrac=capem(i)-cape(i)
+c         defrac=max(defrac,0.001)
+c         frac(i)=-cape(i)/defrac
+c         frac(i)=min(frac(i),1.0)
+c         frac(i)=max(frac(i),0.0)
+c540   continue
+c
+c J Teixeira fix
+c
+      call zilch(byp,ncum)
+      do 515 i=1,ncum
+        lcape(i)=.true.
+ 515  continue
+      do 530 k=minorig+1,nl-1
+        do 520 i=1,ncum
+          if(cape(i).lt.0.0)lcape(i)=.false.
+          if((k.ge.(icb(i)+1)).and.lcape(i))then
+            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+            cape(i)=cape(i)+by
+            if(by.ge.0.0)inb1(i)=k+1
+            if(cape(i).gt.0.0)then
+              inb(i)=k+1
+              capem(i)=cape(i)
+            endif
+          endif
+ 520    continue
+ 530  continue
+      do 540 i=1,ncum
+          cape(i)=capem(i)+byp(i)
+          defrac=capem(i)-cape(i)
+          defrac=max(defrac,0.001)
+          frac(i)=-cape(i)/defrac
+          frac(i)=min(frac(i),1.0)
+          frac(i)=max(frac(i),0.0)
+ 540  continue
+c
+c=====================================================================
+c ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
+c=====================================================================
+c
+      do 600 k=minorig+1,nl
+        do 590 i=1,ncum
+        if((k.ge.icb(i)).and.(k.le.inb(i)))then
+          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
+        endif
+ 590    continue
+ 600  continue
+c
+c=====================================================================
+c ---  CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),
+c --- AT EACH MODEL LEVEL
+c=====================================================================
+c
+c tvpplcl = parcel temperature lifted adiabatically from level
+c           icb-1 to the LCL.
+c tvaplcl = virtual temperature at the LCL.
+c
+      do 610 i=1,ncum
+        dtpbl(i)=0.0
+        tvpplcl(i)=tvp(i,icb(i)-1)
+     &  -rd*tvp(i,icb(i)-1)*(p(i,icb(i)-1)-plcl(i))
+     &  /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
+        tvaplcl(i)=tv(i,icb(i))
+     &  +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
+     &  /(p(i,icb(i))-p(i,icb(i)+1))
+ 610  continue
+c
+c-------------------------------------------------------------------
+c --- Interpolate difference between lifted parcel and
+c --- environmental temperatures to lifted condensation level
+c-------------------------------------------------------------------
+c
+c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
+c
+      do 630 k=minorig,icbmax
+        do 620 i=1,ncum
+        if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
+          dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
+        endif
+ 620    continue
+ 630  continue
+      do 640 i=1,ncum
+        dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
+        dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
+ 640  continue
+c
+c-------------------------------------------------------------------
+c --- Adjust cloud base mass flux
+c-------------------------------------------------------------------
+c
+      do 650 i=1,ncum
+       work(i)=cbmf(i)
+       cbmf(i)=max(0.0,(1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
+       if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
+         iflag(i)=3
+       endif
+ 650  continue
+c
+c-------------------------------------------------------------------
+c --- Calculate rates of mixing,  m(i)
+c-------------------------------------------------------------------
+c
+      call zilch(work,ncum)
+c
+      do 670 j=minorig+1,nl
+        do 660 i=1,ncum
+          if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
+             k=min(j,inb1(i))
+             dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
+     &       +entp*0.04*(ph(i,k)-ph(i,k+1))
+             work(i)=work(i)+dbo
+             m(i,j)=cbmf(i)*dbo
+          endif
+ 660    continue
+ 670  continue
+      do 690 k=minorig+1,nl
+        do 680 i=1,ncum
+          if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
+            m(i,k)=m(i,k)/work(i)
+          endif
+ 680    continue
+ 690  continue
+c
+c
+c=====================================================================
+c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
+c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
+c --- FRACTION (sij)
+c=====================================================================
+c
+c
+       do 750 i=minorig+1,nl
+         do 710 j=minorig+1,nl
+           do 700 ij=1,ncum
+             if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
+     &         .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
+               qti=qnk(ij)-ep(ij,i)*clw(ij,i)
+               bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j)
+     &         /(rv*t(ij,j)*t(ij,j)*cpd)
+               anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j))
+               denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j)
+               dei=denom
+               if(abs(dei).lt.0.01)dei=0.01
+               sij(ij,i,j)=anum/dei
+               sij(ij,i,i)=1.0
+               altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
+               altem=altem/bf2
+               cwat=clw(ij,j)*(1.-ep(ij,j))
+               stemp=sij(ij,i,j)
+               if((stemp.lt.0.0.or.stemp.gt.1.0.or.
+     1           altem.gt.cwat).and.j.gt.i)then
+                 anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2)
+                 denom=denom+lv(ij,j)*(q(ij,i)-qti)
+                 if(abs(denom).lt.0.01)denom=0.01
+                 sij(ij,i,j)=anum/denom
+                 altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
+                 altem=altem-(bf2-1.)*cwat
+               endif
+               if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
+                 qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
+     &                        +(1.-sij(ij,i,j))*qti
+                 uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
+     &                        +(1.-sij(ij,i,j))*u(ij,nk(ij))
+                 vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
+     &                        +(1.-sij(ij,i,j))*v(ij,nk(ij))
+                 elij(ij,i,j)=altem
+                 elij(ij,i,j)=max(0.0,elij(ij,i,j))
+                 ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
+                 nent(ij,i)=nent(ij,i)+1
+               endif
+             sij(ij,i,j)=max(0.0,sij(ij,i,j))
+             sij(ij,i,j)=min(1.0,sij(ij,i,j))
+             endif
+  700      continue
+  710    continue
+c
+c   ***   If no air can entrain at level i assume that updraft detrains  ***
+c   ***   at that level and calculate detrained air flux and properties  ***
+c
+           do 740 ij=1,ncum
+             if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
+     &       .and.(nent(ij,i).eq.0))then
+               ment(ij,i,i)=m(ij,i)
+               qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
+               uent(ij,i,i)=u(ij,nk(ij))
+               vent(ij,i,i)=v(ij,nk(ij))
+               elij(ij,i,i)=clw(ij,i)
+               sij(ij,i,i)=1.0
+             endif
+ 740       continue
+ 750   continue
+c
+      do 770 i=1,ncum
+        sij(i,inb(i),inb(i))=1.0
+ 770  continue
+c
+c=====================================================================
+c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
+c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
+c=====================================================================
+c
+c
+       call zilch(bsum,ncum*nlp)
+       do 780 ij=1,ncum
+         lwork(ij)=.false.
+ 780   continue
+       do 789 i=minorig+1,nl
+c
+         num1=0
+         do 779 ij=1,ncum
+           if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
+ 779     continue
+         if(num1.le.0)go to 789
+c
+           do 781 ij=1,ncum
+             if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
+                lwork(ij)=(nent(ij,i).ne.0)
+                qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
+                anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
+                denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
+                if(abs(denom).lt.0.01)denom=0.01
+                scrit(ij)=anum/denom
+                alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
+                if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
+                asij(ij)=0.0
+                smin(ij)=1.0
+             endif
+ 781       continue
+         do 783 j=minorig,nl
+c
+         num2=0
+         do 778 ij=1,ncum
+             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
+     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
+     &       .and.lwork(ij))num2=num2+1
+ 778     continue
+         if(num2.le.0)go to 783
+c
+           do 782 ij=1,ncum
+             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
+     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
+                  if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
+                    if(j.gt.i)then
+                      smid=min(sij(ij,i,j),scrit(ij))
+                      sjmax=smid
+                      sjmin=smid
+                        if(smid.lt.smin(ij)
+     &                  .and.sij(ij,i,j+1).lt.smid)then
+                          smin(ij)=smid
+                          sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
+                          sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
+                          sjmin=min(sjmin,scrit(ij))
+                        endif
+                    else
+                      sjmax=max(sij(ij,i,j+1),scrit(ij))
+                      smid=max(sij(ij,i,j),scrit(ij))
+                      sjmin=0.0
+                      if(j.gt.1)sjmin=sij(ij,i,j-1)
+                      sjmin=max(sjmin,scrit(ij))
+                    endif
+                    delp=abs(sjmax-smid)
+                    delm=abs(sjmin-smid)
+                    asij(ij)=asij(ij)+(delp+delm)
+     &                           *(ph(ij,j)-ph(ij,j+1))
+                    ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
+     &                           *(ph(ij,j)-ph(ij,j+1))
+                  endif
+              endif
+  782    continue
+  783    continue
+            do 784 ij=1,ncum
+            if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
+               asij(ij)=max(1.0e-21,asij(ij))
+               asij(ij)=1.0/asij(ij)
+               bsum(ij,i)=0.0
+            endif
+ 784        continue
+            do 786 j=minorig,nl+1
+              do 785 ij=1,ncum
+                if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
+     &          .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
+     &          .and.lwork(ij))then
+                   ment(ij,i,j)=ment(ij,i,j)*asij(ij)
+                   bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
+                endif
+ 785     continue
+ 786     continue
+             do 787 ij=1,ncum
+               if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
+     &         .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
+                 nent(ij,i)=0
+                 ment(ij,i,i)=m(ij,i)
+                 qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
+                 uent(ij,i,i)=u(ij,nk(ij))
+                 vent(ij,i,i)=v(ij,nk(ij))
+                 elij(ij,i,i)=clw(ij,i)
+                 sij(ij,i,i)=1.0
+               endif
+  787        continue
+  789  continue
+c
+c=====================================================================
+c --- PRECIPITATING DOWNDRAFT CALCULATION
+c=====================================================================
+c
+c   ***  Check whether ep(inb)=0, if so, skip precipitating    ***
+c   ***             downdraft calculation                      ***
+c
+c
+c   ***  Integrate liquid water equation to find condensed water   ***
+c   ***                and condensed water flux                    ***
+c
+c
+      do 890 i=1,ncum
+        jtt(i)=2
+        if(ep(i,inb(i)).le.0.0001)iflag(i)=2
+        if(iflag(i).eq.0)then
+          lwork(i)=.true.
+        else
+          lwork(i)=.false.
+        endif
+ 890  continue
+c
+c    ***                    Begin downdraft loop                    ***
+c
+c
+        call zilch(wdtrain,ncum)
+        do 899 i=nl+1,1,-1
+c
+          num1=0
+          do 879 ij=1,ncum
+            if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
+ 879      continue
+          if(num1.le.0)go to 899
+c
+c
+c    ***        Calculate detrained precipitation             ***
+c
+          do 891 ij=1,ncum
+            if((i.le.inb(ij)).and.(lwork(ij)))then
+            wdtrain(ij)=g*ep(ij,i)*m(ij,i)*clw(ij,i)
+            endif
+ 891      continue
+c
+          if(i.gt.1)then
+            do 893 j=1,i-1
+              do 892 ij=1,ncum
+                if((i.le.inb(ij)).and.(lwork(ij)))then
+                  awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
+                  awat=max(0.0,awat)
+                  wdtrain(ij)=wdtrain(ij)+g*awat*ment(ij,j,i)
+                endif
+ 892          continue
+ 893      continue
+          endif
+c
+c    ***    Find rain water and evaporation using provisional   ***
+c    ***              estimates of qp(i)and qp(i-1)             ***
+c
+c
+c  ***  Value of terminal velocity and coeffecient of evaporation for snow   ***
+c
+          do 894 ij=1,ncum
+            if((i.le.inb(ij)).and.(lwork(ij)))then
+            coeff=coeffs
+            wt(ij,i)=omtsnow
+c
+c  ***  Value of terminal velocity and coeffecient of evaporation for rain   ***
+c
+            if(t(ij,i).gt.273.0)then
+              coeff=coeffr
+              wt(ij,i)=omtrain
+            endif
+            qsm=0.5*(q(ij,i)+qp(ij,i+1))
+            afac=coeff*ph(ij,i)*(qs(ij,i)-qsm)
+     &       /(1.0e4+2.0e3*ph(ij,i)*qs(ij,i))
+            afac=max(afac,0.0)
+            sigt=sigp(ij,i)
+            sigt=max(0.0,sigt)
+            sigt=min(1.0,sigt)
+            b6=100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij,i)
+            c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
+            revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
+            evap(ij,i)=sigt*afac*revap
+            water(ij,i)=revap*revap
+c
+c    ***  Calculate precipitating downdraft mass flux under     ***
+c    ***              hydrostatic approximation                 ***
+c
+            if(i.gt.1)then
+              dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
+              dhdp=max(dhdp,10.0)
+              mp(ij,i)=100.*ginv*lv(ij,i)*sigd*evap(ij,i)/dhdp
+              mp(ij,i)=max(mp(ij,i),0.0)
+c
+c   ***   Add small amount of inertia to downdraft              ***
+c
+              fac=20.0/(ph(ij,i-1)-ph(ij,i))
+              mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
+c
+c    ***      Force mp to decrease linearly to zero                 ***
+c    ***      between about 950 mb and the surface                  ***
+c
+              if(p(ij,i).gt.(0.949*p(ij,1)))then
+                 jtt(ij)=max(jtt(ij),i)
+                 mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
+     &           /(p(ij,1)-p(ij,jtt(ij)))
+              endif
+            endif
+c
+c    ***       Find mixing ratio of precipitating downdraft     ***
+c
+            if(i.ne.inb(ij))then
+              if(i.eq.1)then
+                qstm=qs(ij,1)
+              else
+                qstm=qs(ij,i-1)
+              endif
+              if(mp(ij,i).gt.mp(ij,i+1))then
+                 rat=mp(ij,i+1)/mp(ij,i)
+                 qp(ij,i)=qp(ij,i+1)*rat+q(ij,i)*(1.0-rat)+100.*ginv*
+     &             sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
+                 up(ij,i)=up(ij,i+1)*rat+u(ij,i)*(1.-rat)
+                 vp(ij,i)=vp(ij,i+1)*rat+v(ij,i)*(1.-rat)
+               else
+                 if(mp(ij,i+1).gt.0.0)then
+                   qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
+     &               +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
+     &               *(cl-cpd))+cpd*(t(ij,i+1)-t(ij,i)))
+     &               /(lv(ij,i)+t(ij,i)*(cl-cpd))
+                   up(ij,i)=up(ij,i+1)
+                   vp(ij,i)=vp(ij,i+1)
+                 endif
+              endif
+              qp(ij,i)=min(qp(ij,i),qstm)
+              qp(ij,i)=max(qp(ij,i),0.0)
+            endif
+            endif
+ 894      continue
+ 899    continue
+c
+c   ***  Calculate surface precipitation in mm/day     ***
+c
+        do 1190 i=1,ncum
+          if(iflag(i).le.1)then
+cc            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
+cc     &                /(rowl*g)
+cc            precip(i)=precip(i)*delt/86400.
+            precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
+          endif
+ 1190   continue
+c
+c
+c   ***  Calculate downdraft velocity scale and surface temperature and  ***
+c   ***                    water vapor fluctuations                      ***
+c
+c     wd=beta*abs(mp(icb))*0.01*rd*t(icb)/(sigd*p(icb))
+c     qprime=0.5*(qp(1)-q(1))
+c     tprime=lv0*qprime/cpd
+c
+c   ***  Calculate tendencies of lowest level potential temperature  ***
+c   ***                      and mixing ratio                        ***
+c
+        do 1200 i=1,ncum
+          work(i)=0.01/(ph(i,1)-ph(i,2))
+          am(i)=0.0
+ 1200   continue
+        do 1220 k=2,nl
+          do 1210 i=1,ncum
+            if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
+              am(i)=am(i)+m(i,k)
+            endif
+ 1210     continue
+ 1220   continue
+        do 1240 i=1,ncum
+          if((g*work(i)*am(i)).ge.delti)iflag(i)=1
+          ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1)
+     &    +(gz(i,2)-gz(i,1))/cpn(i,1))
+          ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1)
+          ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2)
+     &     -t(i,1))*work(i)/cpn(i,1)
+          fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))*
+     &    work(i)+sigd*evap(i,1)
+          fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i)
+          fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))
+     &    +am(i)*(u(i,2)-u(i,1)))
+          fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))
+     &    +am(i)*(v(i,2)-v(i,1)))
+ 1240   continue
+        do 1260 j=2,nl
+           do 1250 i=1,ncum
+             if(j.le.inb(i))then
+               fq(i,1)=fq(i,1)
+     &                 +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1))
+               fu(i,1)=fu(i,1)
+     &                 +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1))
+               fv(i,1)=fv(i,1)
+     &                 +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1))
+             endif
+ 1250      continue
+ 1260   continue
+c
+c   ***  Calculate tendencies of potential temperature and mixing ratio  ***
+c   ***               at levels above the lowest level                   ***
+c
+c   ***  First find the net saturated updraft and downdraft mass fluxes  ***
+c   ***                      through each level                          ***
+c
+        do 1500 i=2,nl+1
+c
+          num1=0
+          do 1265 ij=1,ncum
+            if(i.le.inb(ij))num1=num1+1
+ 1265     continue
+          if(num1.le.0)go to 1500
+c
+          call zilch(amp1,ncum)
+          call zilch(ad,ncum)
+c
+          do 1280 k=i+1,nl+1
+            do 1270 ij=1,ncum
+              if((i.ge.nk(ij)).and.(i.le.inb(ij))
+     &            .and.(k.le.(inb(ij)+1)))then
+                amp1(ij)=amp1(ij)+m(ij,k)
+              endif
+ 1270         continue
+ 1280     continue
+c
+          do 1310 k=1,i
+            do 1300 j=i+1,nl+1
+               do 1290 ij=1,ncum
+                 if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
+                   amp1(ij)=amp1(ij)+ment(ij,k,j)
+                 endif
+ 1290          continue
+ 1300       continue
+ 1310     continue
+          do 1340 k=1,i-1
+            do 1330 j=i,nl+1
+              do 1320 ij=1,ncum
+                if((i.le.inb(ij)).and.(j.le.inb(ij)))then
+                   ad(ij)=ad(ij)+ment(ij,j,k)
+                endif
+ 1320         continue
+ 1330       continue
+ 1340     continue
+c
+          do 1350 ij=1,ncum
+          if(i.le.inb(ij))then
+            dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
+            cpinv=1.0/cpn(ij,i)
+c
+            ft(ij,i)=ft(ij,i)
+     &       +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i)
+     &       +(gz(ij,i+1)-gz(ij,i))*cpinv)
+     &       -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv))
+     &       -sigd*lvcp(ij,i)*evap(ij,i)
+            ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+
+     &        t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
+            ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)*
+     &        (t(ij,i+1)-t(ij,i))*dpinv*cpinv
+            fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))-
+     &        ad(ij)*(q(ij,i)-q(ij,i-1)))
+            fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))-
+     &        ad(ij)*(u(ij,i)-u(ij,i-1)))
+            fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))-
+     &        ad(ij)*(v(ij,i)-v(ij,i-1)))
+         endif
+ 1350    continue
+         do 1370 k=1,i-1
+           do 1360 ij=1,ncum
+             if(i.le.inb(ij))then
+               awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
+               awat=max(awat,0.0)
+               fq(ij,i)=fq(ij,i)
+     &         +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
+               fu(ij,i)=fu(ij,i)
+     &         +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
+               fv(ij,i)=fv(ij,i)
+     &         +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
+             endif
+ 1360      continue
+ 1370    continue
+         do 1390 k=i,nl+1
+           do 1380 ij=1,ncum
+             if((i.le.inb(ij)).and.(k.le.inb(ij)))then
+               fq(ij,i)=fq(ij,i)
+     &                  +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
+               fu(ij,i)=fu(ij,i)
+     &                  +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
+               fv(ij,i)=fv(ij,i)
+     &                  +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
+             endif
+ 1380      continue
+ 1390    continue
+          do 1400 ij=1,ncum
+           if(i.le.inb(ij))then
+             fq(ij,i)=fq(ij,i)
+     &                +sigd*evap(ij,i)+g*(mp(ij,i+1)*
+     &                (qp(ij,i+1)-q(ij,i))
+     &                -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
+             fu(ij,i)=fu(ij,i)
+     &                +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)*
+     &                (up(ij,i)-u(ij,i-1)))*dpinv
+             fv(ij,i)=fv(ij,i)
+     &               +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
+     &               (vp(ij,i)-v(ij,i-1)))*dpinv
+           endif
+ 1400     continue
+ 1500   continue
+c
+c   *** Adjust tendencies at top of convection layer to reflect  ***
+c   ***       actual position of the level zero cape             ***
+c
+        do 503 ij=1,ncum
+        fqold=fq(ij,inb(ij))
+        fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
+        fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
+     &   +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
+     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
+     &   /lv(ij,inb(ij)-1)
+        ftold=ft(ij,inb(ij))
+        ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
+        ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
+     &   +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
+     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
+     &   /cpn(ij,inb(ij)-1)
+        fuold=fu(ij,inb(ij))
+        fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
+        fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
+     &   +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
+     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
+        fvold=fv(ij,inb(ij))
+        fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
+        fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
+     &  +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
+     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
+ 503    continue
+c
+c   ***   Very slightly adjust tendencies to force exact   ***
+c   ***     enthalpy, momentum and tracer conservation     ***
+c
+        do 682 ij=1,ncum
+        ents(ij)=0.0
+        uav(ij)=0.0
+        vav(ij)=0.0
+        do 681 i=1,inb(ij)
+         ents(ij)=ents(ij)
+     &  +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1))	
+         uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
+         vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
+  681	continue
+  682   continue
+        do 683 ij=1,ncum
+        ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
+        uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
+        vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
+ 683    continue
+        do 642 ij=1,ncum
+        do 641 i=1,inb(ij)
+         ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
+         fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
+         fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
+  641	continue
+ 642    continue
+c
+        do 1810 k=1,nl+1
+          do 1800 i=1,ncum
+            if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
+ 1800     continue
+ 1810   continue
+c
+c
+        do 1900 i=1,ncum
+          if(iflag(i).gt.2)then
+          precip(i)=0.0
+          cbmf(i)=0.0
+          endif
+ 1900   continue
+        do 1920 k=1,nl
+         do 1910 i=1,ncum
+           if(iflag(i).gt.2)then
+             ft(i,k)=0.0
+             fq(i,k)=0.0
+             fu(i,k)=0.0
+             fv(i,k)=0.0
+           endif
+ 1910    continue
+ 1920   continue
+        do 2000 i=1,ncum
+         precip1(idcum(i))=precip(i)
+         cbmf1(idcum(i))=cbmf(i)
+         iflag1(idcum(i))=iflag(i)
+ 2000   continue
+        do 2020 k=1,nl
+          do 2010 i=1,ncum
+            ft1(idcum(i),k)=ft(i,k)
+            fq1(idcum(i),k)=fq(i,k)
+            fu1(idcum(i),k)=fu(i,k)
+            fv1(idcum(i),k)=fv(i,k)
+ 2010     continue
+ 2020   continue
+c
+      DO k=1,nd
+        DO i=1,len
+         Ma(i,k) = 0.
+        ENDDO
+      ENDDO
+      DO k=nl,1,-1
+        DO i=1,ncum
+          Ma(i,k) = Ma(i,k+1)+m(i,k)
+        ENDDO
+      ENDDO
+c
+        return
+        end
+
Index: /LMDZ.3.3/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/physiq.F	(revision 253)
+++ /LMDZ.3.3/trunk/libf/phylmd/physiq.F	(revision 254)
@@ -472,4 +472,5 @@
       SAVE bbase
       REAL rflag(klon)          ! flag fonctionnement de convect
+      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
 c -- convect43:
       INTEGER ntra              ! nb traceurs pour convect4.3
@@ -691,5 +692,5 @@
 cKE43
 c Initialisation pour la convection de K.E. (sb):
-         IF (iflag_con.EQ.4) THEN
+         IF (iflag_con.GE.3) THEN
 
          PRINT*, "*** Convection de Kerry Emanuel 4.3  "
@@ -1086,5 +1087,5 @@
      .                "ave(X)", zsto,zout)
 cKE43
-      IF (iflag_con .EQ. 4) THEN ! sb
+      IF (iflag_con .GE. 3) THEN ! sb
 c
          CALL histdef(nid_mth, "cape", "Conv avlbl pot ener", "J/Kg",
@@ -1271,5 +1272,5 @@
 c
 cKE43
-      IF (iflag_con.EQ.4) THEN ! (sb)
+      IF (iflag_con.GE.3) THEN ! (sb)
 c
          CALL histdef(nid_mth, "upwd", "saturated updraft", "Kg/m2/s",
@@ -1713,11 +1714,5 @@
          itop_con(i) = klev+1 - kctop(i)
       ENDDO
-      ELSE IF (iflag_con.EQ.3) THEN
-          stop'reactiver le call conlmd dans physiq.F'
-c     CALL conccm (dtime,paprs,pplay,t_seri,q_seri,conv_q,
-c    s             d_t_con, d_q_con, 
-c    s             rain_con, snow_con, ibas_con, itop_con)
-cKE43
-      ELSE IF (iflag_con.EQ.4) THEN
+      ELSE IF (iflag_con.GE.3) THEN
 c nb of tracers for the KE convection:
           if (nqmax .GE. 4) then
@@ -1726,25 +1721,30 @@
               ntra = 1 
           endif
-cke43 (arguments inutiles enleves => des SAVE dans conema43?)
-c$$$          CALL conema43(dtime,paprs,pplay,t_seri,q_seri,
-c$$$     $        u_seri,v_seri,tr_seri,nbtr,
-c$$$     .        ema_workcbmf,
-c$$$     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
-c$$$     .        wdn, tdn, qdn,
-c$$$     .        rain_con, snow_con, ibas_con, itop_con,
-c$$$     .        upwd,dnwd,dnwd0,bas,top,Ma,cape,tvp,rflag,
-c$$$     .        pbase
-c$$$     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
-c$$$     .        pori_con,plcl_con,dtma_con,dtlcl_con)
-          CALL conema (dtime,paprs,pplay,t_seri,q_seri,
-     $        u_seri,v_seri,tr_seri,nbtr,
+          if (iflag_con.eq.4) then ! vectorise
+          CALL conemav (dtime,paprs,pplay,t_seri,q_seri,
+     .        u_seri,v_seri,tr_seri,nbtr,
      .        ema_work1,ema_work2,
      .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
-c$$$     .        wdn, tdn, qdn,
      .        rain_con, snow_con, ibas_con, itop_con,
-     .        upwd,dnwd,dnwd0,bas,top,Ma,cape,tvp,rflag,
-     .        pbase
+     .        upwd,dnwd,dnwd0,
+c    .        Ma,cape,tvp,(/(nint(rflag(i)),i=1,size(rflag))/),
+     .        Ma,cape,tvp,iflagctrl,
+     .       pbase
      .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
-c$$$     .        pori_con,plcl_con,dtma_con,dtlcl_con)
+
+          else
+
+          CALL conema (dtime,paprs,pplay,t_seri,q_seri,
+     .        u_seri,v_seri,tr_seri,nbtr,
+     .        ema_work1,ema_work2,
+     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
+     .        rain_con, snow_con, ibas_con, itop_con,
+     .        upwd,dnwd,dnwd0,bas,top,
+     .        Ma,cape,tvp,rflag,
+     .       pbase
+     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
+          endif
+
+
           DO i = 1, klon
             ema_pcb(i)  = pbase(i)
@@ -1809,5 +1809,5 @@
       IF (nqmax.GT.2) THEN !--melange convectif de traceurs
 c
-          IF (iflag_con .NE. 2 .AND.  iflag_con .NE. 4 ) THEN 
+          IF (iflag_con .NE. 2 .AND.  iflag_con .LT. 3 ) THEN 
               PRINT*, 'Pour l instant, seul conflx fonctionne ',
      $            'avec traceurs', iflag_con
@@ -2101,5 +2101,5 @@
 C Pascale : il faut quand meme apeller phytrac car il gere les sorties
 cKE43       des traceurs => il faut donc mettre des flags a .false.
-      IF (iflag_con.EQ.4) THEN
+      IF (iflag_con.GE.3) THEN
 c           on ajoute les tendances calculees par KE43
         DO iq=1, nqmax-2 ! Sandrine a -3 ???
@@ -2440,5 +2440,5 @@
       CALL histwrite(nid_mth,"vq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 cKE43
-      IF (iflag_con .EQ. 4) THEN ! sb
+      IF (iflag_con .GE. 3) THEN ! sb
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, cape,zx_tmp_2d)
@@ -2542,9 +2542,9 @@
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zpt_conv, zx_tmp_3d)
       CALL histwrite(nid_mth,"ptconv",itap,zx_tmp_3d,
-     .                                   iim*jjmp1*klev,ndex3d)
+     .                                   iim*(jjmp1)*klev,ndex3d)
 c
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, ratqs, zx_tmp_3d)
       CALL histwrite(nid_mth,"ratqs",itap,zx_tmp_3d,
-     .                                   iim*jjmp1*klev,ndex3d)
+     .                                   iim*(jjmp1)*klev,ndex3d)
 c
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
@@ -2619,5 +2619,5 @@
       ENDIF
 cKE43
-      IF (iflag_con.EQ.4) THEN ! (sb)
+      IF (iflag_con.GE.3) THEN ! (sb)
 c
       CALL gr_fi_ecrit(klev,klon,iim,jjmp1, upwd, zx_tmp_3d)
Index: /LMDZ.3.3/trunk/libf/phylmd/printflag.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/printflag.F	(revision 253)
+++ /LMDZ.3.3/trunk/libf/phylmd/printflag.F	(revision 254)
@@ -35,5 +35,8 @@
      ,          ******'
        ELSE IF ( iflag_con.EQ. 3 )   THEN
-           PRINT *,' *****           Shema  convection    CCM      
+           PRINT *,' *****   Shema  convection Emanuel version LMD
+     ,          ******'
+       ELSE IF ( iflag_con.EQ. 4 )   THEN
+           PRINT *,' *****   Shema  convection Emanuel version 4vec
      ,          ******'
        ENDIF
Index: /LMDZ.3.3/trunk/libf/phylmd/zilch.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/zilch.F	(revision 254)
+++ /LMDZ.3.3/trunk/libf/phylmd/zilch.F	(revision 254)
@@ -0,0 +1,13 @@
+      subroutine zilch(x,m)
+c
+c Zero the real array x dimensioned m.
+c
+      implicit none
+c
+      integer m,i
+      real x(m)
+      do 1 i=1,m
+      x(i)= 0.0  
+    1 continue
+      return
+      end
