! $Header$

    cape, tvp, pbase, buoybase, &  ! ccc     *

                                   ! DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)

    dtvpdt1, dtvpdq1, dplcldt, dplcldr, & ! sbl

    ft2, fr2, fu2, fv2, wd, qcond, qcondc) ! sbl

  ! ***  THE PARAMETER NA SHOULD IN GENERAL EQUAL ND   ***

  ! #################################################################

  ! Fleur       Introduction des traceurs dans convect3 le 6 juin 200

  ! #################################################################

  USE dimphy

  USE infotrac_phy, ONLY: nbtr

  IMPLICIT NONE

  INTEGER na

  PARAMETER (na=60)

  REAL deltac ! cld

  PARAMETER (deltac=0.01) ! cld

  INTEGER nent(na)

  INTEGER nd, ndp1, nl, ntra, iflag, icb, inb

  REAL dtime, epmax, delt, precip, cape

  REAL dplcldt, dplcldr

  REAL t1(nd), r1(nd), rs(nd), u(nd), v(nd), tra(nd, ntra)

  REAL p(nd), ph(ndp1)

  REAL ft(nd), fr(nd), fu(nd), fv(nd), ftra(nd, ntra)

  REAL sig(nd), w0(nd)

  REAL m(na), mp(na), ment(na, na), qent(na, na), elij(na, na)

  REAL sij(na, na), tvp(na), tv(na), water(na)

  REAL rp(na), ep(na), th(na), wt(na), evap(na), clw(na)

  REAL sigp(na), b(na), tp(na), cpn(na)

  REAL lv(na), lvcp(na), h(na), hp(na), gz(na)

  REAL t(na), rr(na)

  REAL ft2(nd), fr2(nd), fu2(nd), fv2(nd) ! added sbl

  REAL buoy(na) !  Lifted parcel buoyancy

  REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual

  ! temperature wrt T1 and Q1

  REAL clw_new(na), qi(na)

  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

  LOGICAL ice_conv, ok_adj

  PARAMETER (ice_conv=.TRUE.)

  ! ccccccccccccccccccccccccccccccccccccccccccccc

  ! declaration des variables a sortir

  ! cccccccccccccccccccccccccccccccccccccccccccc

  REAL mke(nd)

  REAL mike(nd)

  REAL ma(nd)

  REAL tps(nd) !temperature dans les ascendances non diluees

  REAL tls(nd) !temperature potentielle

  REAL ments(nd, nd)

  REAL qents(nd, nd)

  REAL sigij(klev, klev)

  REAL pbase ! pressure at the cloud base level

  REAL buoybase ! buoyancy at the cloud base level

  ! ccccccccccccccccccccccccccccccccccccccccccccc

  REAL :: cpv,cl,cpvmcl,eps,alv0,rdcp,pbcrit,ptcrit,sigd,spfac

  REAL :: tau,beta,alpha,dtcrit,dtovsh,ahm,rm,um,vm,dphinv

  REAL :: a2,x,tvx,tvy,plcl,pden,dpbase,tvpbase,tvbase,tdif

  REAL :: ath1,ath,delti,deltap,dcape,dlnp,sigold,dtmin,fac,w

  REAL :: amu,rti,cpd,bf2,anum,denom,dei,altem,cwat,stemp,qp

  REAL :: scrit,alt,smax,asij,wgh,sjmax,sjmin,smid,delp,delm

  REAL :: asum,bsum,csum,wflux,tinv,wdtrain,awat,afac,afac1,afac2

  REAL :: bfac,pr1,pr2,sigt,b6,c6,revap,tevap,delth,amfac,amp2

  REAL :: xf,tf,af,bf,fac2,ur,sru,d,ampmax,dpinv,am,amde,cpinv

  REAL :: amp1,ad,rat,ax,bx,cx,dx,ex,dsum

  INTEGER :: nk,i,j,nopt,jn,k,im,jm,n

  REAL dnwd0(nd) !  precipitation driven unsaturated downdraft flux

  REAL dnwd(nd), dn1 ! in-cloud saturated downdraft mass flux

  REAL upwd(nd), up1 ! in-cloud saturated updraft mass flux

  ! ***         ASSIGN VALUES OF THERMODYNAMIC CONSTANTS        ***

  ! ***             THESE SHOULD BE CONSISTENT WITH             ***

  ! ***              THOSE USED IN CALLING PROGRAM              ***

  ! ***     NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT  ***

  ! sb      CPD=1005.7

  ! sb      CPV=1870.0

  ! sb      CL=4190.0

  ! sb      CPVMCL=CL-CPV

  ! sb      RV=461.5

  ! sb      RD=287.04

  ! sb      EPS=RD/RV

  ! sb      ALV0=2.501E6

  ! cccccccccccccccccccccc

  ! constantes coherentes avec le modele du Centre Europeen

  ! sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644

  ! sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153

  ! sb      CPD = 3.5 * RD

  ! sb      CPV = 4.0 * RV

  ! sb      CL = 4218.0

  ! sb      CPVMCL=CL-CPV

  ! sb      EPS=RD/RV

  ! sb      ALV0=2.5008E+06

  ! ccccccccccccccccccccc

  ! on utilise les constantes thermo du Centre Europeen: (SB)

  include "YOMCST.h"

  ! ccccccccccccccccccccc

  ! ***  INITIALIZE OUTPUT ARRAYS AND PARAMETERS  ***

    END DO

  END DO

  END DO

  ! ************************************************************

  ! *    CALCUL DES TEMPERATURES POTENTIELLES A SORTIR

  ! ************************************************************

  END DO

  ! ***********************************************************

  ! ***                    SPECIFY PARAMETERS                        ***

  ! ***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE   ***

  ! ***       PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO         ***

  ! ***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP.      ***

  ! ***            EFFICIENCY IS ASSUMED TO BE UNITY                 ***

  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***

  ! ***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE      ***

  ! ***                        OF CLOUD                              ***

  ! ***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF    ***

  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***

  ! ***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY)    ***

  ! ***           (BETA MUST BE LESS THAN OR EQUAL TO 1)             ***

  ! ***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE    ***

  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***

  ! ***                     IT MUST BE LESS THAN 0                   ***

  pbcrit = 150.0

  ptcrit = 500.0

  sigd = 0.01

  spfac = 0.15

  ! sb:

  ! EPMAX=0.993 ! precip efficiency less than unity

  ! EPMAX=1. ! precip efficiency less than unity

  ! jyg

  ! CC      BETA=0.96

  ! Beta is now expressed as a function of the characteristic time

  ! of the convective process.

  ! CC        Old value : TAU = 15000.   !(for dtime = 600.s)

  ! CC        Other value (inducing little change) :TAU = 8000.

  tau = 8000.

  ! jyg

  ! CC      ALPHA=1.0

  ! Increase alpha in order to compensate W decrease

  ! jyg (voir CONVECT 3)

  ! CC      DTCRIT=-0.2

  dtcrit = -2.

  ! gf&jyg

  ! CC     DT pour l'overshoot.

  dtovsh = -0.2

  ! ***        INCREMENT THE COUNTER       ***

  ! ***    IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT     ***

  ! ***     RETURNS ARRAYS T AND R THAT MAY HAVE BEEN      ***

  ! ***  ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE    ***

  ! ***        THE RETURNED ARRAYS ARE UNALTERED.          ***

  nopt = 0

  ! !      NOPT=1 ! sbl

  ! ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***

  ! ***  DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A  ***

  ! ***                BOUNDARY LAYER SCHEME !!!               ***

      jn = 0

      END DO

      ahm = 0.0

      rm = 0.0

      um = 0.0

      vm = 0.0

      END DO

        END DO

      END DO

      END DO

      a2 = 0.0

        END DO

      END DO

      END DO

  END IF ! added sbl

  ! ***   RESET INPUT ARRAYS IF ok_adj 0   ***

      ! !            T1(I)=T(I)      ! commente sbl

      ! !            R1(I)=RR(I)     ! commente sbl

    END DO

  END IF

  ! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY

  END DO

  ! ***  CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL ***

  ! ***       (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)     ***

    ! sb3d         print*,'je suis passe par 366'

  END IF

  ! jyg1 Utilisation de la subroutine CLIFT

  ! C      RH=RR(1)/RS(1)

  ! C      CHI=T(1)/(1669.0-122.0*RH-T(1))

  ! C      PLCL=P(1)*(RH**CHI)

  ! jyg2

  ! sb3d      PRINT *,' em_plcl,p1,t1,r1,rs1,rh '

  ! sb3d     $        ,PLCL,P(1),T(1),RR(1),RS(1),RH

  END IF

  ! jyg1

  ! Essais de modification de ICB

  ! ***  CALCULATE FIRST LEVEL ABOVE LCL (=ICB)  ***

  ! C      ICB=NL-1

  ! C      DO 50 I=2,NL-1

  ! C         IF(P(I).LT.PLCL)THEN

  ! C            ICB=MIN(ICB,I)   ! ICB sup ou egal a 2

  ! C         END IF

  ! C   50 CONTINUE

  ! C      IF(ICB.EQ.(NL-1))THEN

  ! C         IFLAG=3

  ! C         RETURN

  ! C      END IF

  ! *** CALCULATE LAYER CONTAINING LCL (=ICB)   ***

  ! sb      DO 50 I=2,NL-1

    ! la modification consiste a comparer PLCL a PH et non a P:

    ! ICB est defini par :  PH(ICB)<PLCL<PH(ICB-!)

    END IF

  END DO

  END IF

  ! jyg2

  ! *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL

  ! ***

  ! ***  TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC

  ! ***

  ! ***                   LIQUID WATER CONTENT

  ! ***

  ! jyg1

  ! make sure that "Cloud base" seen by TLIFT is actually the

  ! fisrt level where adiabatic ascent is saturated

    ! sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL)

    CALL tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tp, clw, nd, nl, &

  ELSE

    ! sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL)

    CALL tlift(p, t, rr, rs, gz, plcl, icb+1, nk, tvp, tp, clw, nd, nl, &

  END IF

  ! jyg2

  ! *****************************************************************************

  ! ***     SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE

  ! *****************************************************************************

  END DO

  ! *****************************************************************************

  ! ***  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)     ***

    pden = ptcrit - pbcrit

    ! jyg

    ! cc         EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN

    ! sb         EP(I)=(PLCL-P(I)-PBCRIT)/PDEN

    ! sb         EP(I)=AMIN1(EP(I),1.0)

    ! ***       CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL     ***

    ! ***                    VIRTUAL TEMPERATURE                    ***

    ! cd1

    ! . Keep all liquid water in lifted parcel (-> adiabatic CAPE)

    ! cc    TVP(I)=TVP(I)-TP(I)*(RR(1)-EP(I)*CLW(I))

    ! !!! sb         TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift

    ! cd2

    ! ***       Calculate first estimate of buoyancy

  END DO

  ! ***   Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy

  dpbase = -40. !That is 400m above LCL

  tvpbase = tvp(icb)*(pbase-p(icb+1))/(p(icb)-p(icb+1)) + &

  tvbase = tv(icb)*(pbase-p(icb+1))/(p(icb)-p(icb+1)) + &

  ! test sb:

  ! @      write(*,*) '++++++++++++++++++++++++++++++++++++++++'

  ! @      write(*,*) 'plcl,dpbas,tvpbas,tvbas,tvp(icb),tvp(icb1)

  ! @     :             ,tv(icb),tv(icb1)'

  ! @      write(*,*) plcl,dpbase,tvpbase,tvbase,tvp(icb)

  ! @     L          ,tvp(icb+1),tv(icb),tv(icb+1)

  ! @      write(*,*) '++++++++++++++++++++++++++++++++++++++++'

  ! fin test sb

  ! C       Set buoyancy = BUOYBASE for all levels below BASE.

  ! C       For safety, set : BUOY(ICB) = BUOYBASE

    END IF

  END DO

  ! sb3d      print *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl'

  ! sb3d     $,        buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl

  ! sb3d      print *,'TVP ',(tvp(i),i=1,nl)

  ! sb3d      print *,'TV ',(tv(i),i=1,nl)

  ! sb3d      print *, 'P ',(p(i),i=1,nl)

  ! sb3d      print *,'ICB ',icb

  ! test sb:

  ! @      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'

  ! @      write(*,*) 'icb,icbs,inb,buoybase:'

  ! @      write(*,*) icb,icb+1,inb,buoybase

  ! @      write(*,*) 'k,tvp,tv,tp,buoy,ep: '

  ! @      do k=1,nl

  ! @      write(*,*) k,tvp(k),tv(k),tp(k),buoy(k),ep(k)

  ! @      enddo

  ! @      write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'

  ! fin test sb

  ! ***   MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE  ***

  ! ***    AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT  ***

  ! ***                         AT CLOUD BASE                         ***

  ! ***       IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING       ***

  ! ***                        SIG(I) AND W0(I)                       ***

  ! jyg

  ! CC      TDIF=TVP(ICB)-TV(ICB)

  tdif = buoy(icb)

  ! jyg

  ! CC      ATH=TH(ICB-1)-1.0

  ! ATH=0.                          ! ajout sb

  ! IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb

    END DO

  END IF

  ! ***  IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY

  ! ***

  ! ***        NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS

  ! ***

    END DO

      END DO

    END DO

  END DO

  ! ***  FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S       ***

  ! ***                LEVEL OF NEUTRAL BUOYANCY                   ***

    ! jyg

    ! CC         IF((TVP(I)-TV(I)).LT.DTCRIT)THEN

    END IF

  END DO

  ! ***          RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB       ***

      ! jyg

      ! CC            SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(INB)-TVP(INB))*

      ! CC     1              ABS(TV(INB)-TVP(INB))

    END DO

  END IF

    ! jyg

    ! CC         SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(ICB)-TVP(ICB))*

    ! CC     1           (TV(ICB)-TVP(ICB))

  END DO

  ! ***    RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME    ***

  ! ***           AND AFTER 10 TIME STEPS OF NO CONVECTION              ***

    END DO

  END IF

  ! ***   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL   ***

  END DO

  ! ***  CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE),  ***

  ! ***     VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY     ***

  ! ***     UNDILUTE UPDRAFT (SIG),  AND UPDRAFT MASS FLUX (M)    ***

    ! jyg1

    ! CC         CAPE=CAPE+RD*(TVP(I-1)-TV(I-1))*(PH(I-1)-PH(I))/P(I-1)

    ! CC         DCAPE=RD*BUOY(I-1)*(PH(I-1)-PH(I))/P(I-1)

    ! CC         DLNP=(PH(I-1)-PH(I))/P(I-1)

    ! The interval on which CAPE is computed starts at PBASE :

    dcape = rd*buoy(i-1)*deltap/p(i-1)

    dlnp = deltap/p(i-1)

    ! jyg2

    ! sb3d         print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape

    ! test sb:

    ! @       write(*,*) '############################################'

    ! @         write(*,*) 'cape,rrd,buoy,deltap,p,pbase,ph:'

    ! @     :     ,cape,rd,buoy(i-1),deltap,p(i-1),pbase,ph(i)

    ! @       write(*,*) '############################################'

    ! fin test sb

    dtmin = 100.0

      ! jyg

      ! CC            DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J)))

    END DO

    ! sb3d     print *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I)

    ! jyg

    ! C    Essais de reduction de la vitesse

    ! C         FAC = FAC*.5

    ! --------- test sb:

    ! write(*,*) '############################################'

    ! write(*,*) 'k,amu,buoy(k-1),deltap,w,beta,fac,cape,w0(k)'

    ! write(*,*) i,amu,buoy(i-1),deltap

    ! :           ,w,beta,fac,cape,w0(i)

    ! write(*,*) '############################################'

    ! ---------

  END DO

  ! jyg1

  ! sb3d      print *, 'Cloud base, c. top, CAPE',ICB,INB,cape

  ! sb3d      print *, 'SIG ',(sig(i),i=1,inb)

  ! sb3d      print *, 'W ',(w0(i),i=1,inb)

  ! sb3d      print *, 'M ',(m(i), i=1,inb)

  ! sb3d      print *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb)

  ! sb3d      print *, 'Dt_vrai ',(buoy(i),i=1,inb)

  ! jyg2

  ! ***  CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING  ***

  ! ***     RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING     ***

  ! ***                        FRACTION (SIJ)                          ***

      dei = denom

      stemp = sij(i, j)

      END IF

        END DO

        elij(i, j) = altem

      END IF

    END DO

    ! ***   IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS

    ! ***

    ! ***   AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES

    ! ***

      END DO

    END IF

    END DO

  END DO

  ! ***  NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL  ***

  ! ***              PROBABILITIES OF MIXING                     ***

      smax = 0.0

      asij = 0.0

          wgh = 1.0

          ELSE

            sjmin = 0.0

          END IF

        END IF

      END DO

      END DO

      asum = 0.0

      bsum = 0.0

      END DO

      END DO

      csum = 0.0

      END DO

        END DO

      END IF

    END IF

  END DO

  ! **************************************************************

  ! *       CALCUL DES MENTS(I,J) ET DES QENTS(I,J)

  ! *************************************************************

    END DO

  END DO

  ! **********************************************************

  ! --- test sb:

  ! @       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'

  ! @       write(*,*) 'inb,m(inb),ment(inb,inb),sigij(inb,inb):'

  ! @       write(*,*) inb,m(inb),ment(inb,inb),sigij(inb,inb)

  ! @       write(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'

  ! ---

  ! ***  CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING    ***

  ! ***             DOWNDRAFT CALCULATION                      ***

  ! ***  INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER   ***

  ! ***                AND CONDENSED WATER FLUX                    ***

  wflux = 0.0

  tinv = 1./3.

  ! ***                    BEGIN DOWNDRAFT LOOP                    ***

    ! ***              CALCULATE DETRAINED PRECIPITATION             ***

      END DO

    END IF

    ! ***    FIND RAIN WATER AND EVAPORATION USING PROVISIONAL   ***

    ! ***              ESTIMATES OF RP(I)AND RP(I-1)             ***

    END IF

    ELSE

    END IF

    bfac = 1./(sigd*wt(i))

    ! jyg1

    ! CC        SIGT=1.0

    ! CC        IF(I.GE.ICB)SIGT=SIGP(I)

    ! Prise en compte de la variation progressive de SIGT dans

    ! les couches ICB et ICB-1:

    ! Pour PLCL<PH(I+1), PR1=0 & PR2=1

    ! Pour PLCL>PH(I),   PR1=1 & PR2=0

    ! Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval

    ! sur le nuage, et PR2 est la proportion sous la base du

    ! nuage.

    ! sb3d         print *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2

    ! jyg2

    ELSE

      evap(i) = -evap(i+1) + 0.02*(wdtrain+sigd*wt(i)*water(i+1))/(sigd*(ph(i &

    END IF

    ! ***  CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER     ***

    ! ***              HYDROSTATIC APPROXIMATION                 ***

    ! ***           IF HYDROSTATIC ASSUMPTION FAILS,             ***

    ! ***   SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA  ***

    ! ***  AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS ***

      bf = 2.*(tinv*mp(i+1))**3 + tinv*mp(i+1)*xf*tf + &

      fac2 = 1.0

        fac = 1.0

        mp(i) = mp(i+1)*tinv + (0.5*bf+sru)**tinv + &

      ELSE

      END IF

      b(i-1) = b(i) + 100.0*(p(i-1)-p(i))*tevap/(mp(i)+sigd*0.1) - &

    END IF

    ! ***         LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION      ***

    ! ***      FORCE MP TO DECREASE LINEARLY TO ZERO                 ***

    ! ***       BETWEEN CLOUD BASE AND THE SURFACE                   ***