SUBROUTINE MF1_mf(                       &
                 & kts,kte,dt,zw,dz,p,      &
                 & momentum_opt,            &
                 & tke_opt,                 &
                 & u,v,w,th,thl,thv,tk,     &
                 & qt,qv,qc,qke,            &
                 & exner,vt,vq,sgm,         &
                 & ust,flt,flq,flqv,flqc,   &
                 & pblh,kpbl,DX,landsea,ts, &
            ! outputs - tendencies
            !  &dth,dqv,dqc,du,dv,&
            ! outputs - updraft properties   
                 & edmf_a,edmf_w,           &
                 & edmf_qt,edmf_thl,        &
                 & edmf_ent,edmf_qc,        &
            ! outputs - variables needed for solver 
                 & s_aw,s_awthl,s_awqt,     &
                 & s_awqv,s_awqc,           &
                 & s_awu,s_awv,s_awqke,     &
#if (WRF_CHEM == 1)
                 & nchem,chem,s_awchem,      &
#endif
            ! in/outputs - subgrid scale clouds
                 & qc_bl1d,cldfra_bl1d,     &
            ! inputs - flags for moist arrays
                 &F_QC,F_QI,                &
                 &Psig_shcu,                &
            ! output info
                 &nup2,ktop,maxmf,          &
                 ) 

  ! inputs:
     INTEGER, INTENT(IN) :: KTS,KTE,momentum_opt,tke_opt,kpbl

#ifdef HARDCODE_VERTICAL
# define kts 1
# define kte HARDCODE_VERTICAL
#endif

     REAL,DIMENSION(KTS:KTE), INTENT(IN) :: U,V,W,TH,THL,TK,QT,QV,QC,&
                                            THV,P,qke,exner,dz
     REAL,DIMENSION(KTS:KTE+1), INTENT(IN) :: ZW  !height at full-sigma
     REAL, INTENT(IN) :: DT,UST,FLT,FLQ,FLQV,FLQC,PBLH,&
                         DX,Psig_shcu,landsea,ts
     LOGICAL, OPTIONAL :: F_QC,F_QI

  ! outputs - tendencies
  !   REAL,DIMENSION(KTS:KTE), INTENT(OUT) :: DTH,DQV,DQC,DU,DV
  ! outputs - updraft properties
     REAL,DIMENSION(KTS:KTE), INTENT(OUT) :: edmf_a,edmf_w,        &
                      & edmf_qt,edmf_thl, edmf_ent,edmf_qc
  ! output
     INTEGER, INTENT(OUT) :: nup2,ktop
     REAL, INTENT(OUT) :: maxmf
  ! outputs - variables needed for solver
     REAL,DIMENSION(KTS:KTE+1) :: s_aw,      & !sum ai*wis_awphi
                               s_awthl,      & !sum ai*wi*phii
                                s_awqt,      &
                                s_awqv,      &
                                s_awqc,      &
                                 s_awu,      &
                                 s_awv,      &
                               s_awqke

     REAL,DIMENSION(KTS:KTE), INTENT(INOUT) :: qc_bl1d,cldfra_bl1d

    INTEGER, PARAMETER :: NUP=10, debug_mf=0
  ! local variables
  ! updraft properties
     REAL,DIMENSION(KTS:KTE+1,1:NUP) :: UPW,UPTHL,UPQT,UPQC,UPQV,  &
                                        UPA,UPU,UPV,UPTHV,UPQKE
  ! entrainment variables
     REAL,DIMENSION(KTS:KTE,1:NUP) :: ENT,ENTf
     INTEGER,DIMENSION(KTS:KTE,1:NUP) :: ENTi
  ! internal variables
     INTEGER :: K,I,k50
     REAL :: fltv,wstar,qstar,thstar,sigmaW,sigmaQT,sigmaTH,z0,    &
             pwmin,pwmax,wmin,wmax,wlv,wtv,Psig_w,maxw,maxqc,wpbl
     REAL :: B,QTn,THLn,THVn,QCn,Un,Vn,QKEn,Wn2,Wn,EntEXP,EntW,BCOEFF

  ! w parameters
     REAL,PARAMETER :: &
          &Wa=2./3., &
          &Wb=0.002,&
          &Wc=1.5 
        
  ! Lateral entrainment parameters ( L0=100 and ENT0=0.1) were taken from
  ! Suselj et al (2013, jas). Note that Suselj et al (2014,waf) use L0=200 and ENT0=0.2.
     REAL,PARAMETER :: &
         & L0=100.,&
         & ENT0=0.1

  ! Implement ideas from Neggers (2016, JAMES):
     REAL, PARAMETER :: Atot = 0.10 ! Maximum total fractional area of all updrafts
     REAL, PARAMETER :: lmax = 1000.! diameter of largest plume
     REAL, PARAMETER :: dl   = 100. ! diff size of each plume - the differential multiplied by the integrand
     REAL, PARAMETER :: dcut = 1.0  ! max diameter of plume to parameterize relative to dx (km)
     REAL ::  d            != -2.3 to -1.7  ;=-1.9 in Neggers paper; power law exponent for number density (N=Cl^d).
          ! Note that changing d to -2.0 makes each size plume equally contribute to the total coverage of all plumes.
          ! Note that changing d to -1.7 doubles the area coverage of the largest plumes relative to the smallest plumes.
     REAL :: cn,c,l,n,an2,hux,maxwidth,wspd_pbl,cloud_base,width_flx

#if (WRF_CHEM == 1)
     INTEGER, INTENT(IN) :: nchem
     REAL,DIMENSION(kts:kte, nchem) :: chem
     REAL,DIMENSION(kts:kte+1, nchem) :: s_awchem
     REAL,DIMENSION(nchem) :: chemn
     REAL,DIMENSION(KTS:KTE+1,1:NUP, nchem) :: UPCHEM
     INTEGER :: ic
     REAL,DIMENSION(KTS:KTE+1, nchem) :: edmf_chem
#endif

  !JOE: add declaration of ERF
   REAL :: ERF

   LOGICAL :: superadiabatic

  ! VARIABLES FOR CHABOUREAU-BECHTOLD CLOUD FRACTION
   REAL,DIMENSION(KTS:KTE), INTENT(INOUT) :: vt, vq, sgm
   REAL :: sigq,xl,tlk,qsat_tl,rsl,cpm,a,qmq,mf_cf,diffqt,&
           Fng,qww,alpha,beta,bb,f,pt,t,q2p,b9,satvp,rhgrid

   ! WA TEST 11/9/15 for consistent reduction of updraft params
   REAL :: csigma,acfac,EntThrottle

   !JOE- plume overshoot
   INTEGER :: overshoot
   REAL :: bvf, Frz

! Initialize individual updraft properties
  UPW=0.
  UPTHL=0.
  UPTHV=0.
  UPQT=0.
  UPA=0.
  UPU=0.
  UPV=0.
  UPQC=0.
  UPQV=0.
  UPQKE=0.
#if (WRF_CHEM == 1)
  UPCHEM(KTS:KTE+1,1:NUP,1:nchem)=0.0
#endif
  ENT=0.001
! Initialize mean updraft properties
  edmf_a  =0.
  edmf_w  =0.
  edmf_qt =0.
  edmf_thl=0.
  edmf_ent=0.
  edmf_qc =0.
#if (WRF_CHEM == 1)
  edmf_chem(kts:kte+1,1:nchem) = 0.0
#endif
! Initialize the variables needed for implicit solver
  s_aw=0.
  s_awthl=0.
  s_awqt=0.
  s_awqv=0.
  s_awqc=0.
  s_awu=0.
  s_awv=0.
  s_awqke=0.
#if (WRF_CHEM == 1)
  s_awchem(kts:kte+1,1:nchem) = 0.0
#endif


  ! Taper off MF scheme when significant resolved-scale motions
  ! are present This function needs to be asymetric...
  k      = 1
  maxw   = 0.0
  cloud_base  = 9000.0
  DO k=1,kte
     IF(ZW(k) > pblh + 500.) exit

     wpbl = w(k)
     IF(w(k) < 0.)wpbl = 2.*w(k)
     maxw = MAX(maxw,ABS(wpbl))

     !Find highest k-level below 50m AGL
     IF(ZW(k)<=50.)k50=k

     !Search for cloud base
     IF(qc(k)>1E-5 .AND. cloud_base == 9000.0)THEN
       cloud_base = 0.5*(ZW(k)+ZW(k+1))
     ENDIF
  ENDDO
  maxw = MAX(0.,maxw - 0.5)         ! do nothing for small w, but
  Psig_w = MAX(0.0, 1.0 - maxw/0.5) ! linearly taper off for w > 0.5 m/s
  Psig_w = MIN(Psig_w, Psig_shcu)

  fltv = flt + svp1*flq

  !Completely shut off MF scheme for strong resolved-scale vertical velocities.
  IF(Psig_w == 0.0 .and. fltv > 0.0) fltv = -1.*fltv

! if surface buoyancy is positive we do integration, otherwise not, and make sure that 
! PBLH > twice the height of the surface layer (set at z0 = 50m)
! Also, ensure that it is at least slightly superadiabatic up through 50 m
      superadiabatic = .false.
  IF((landsea-1.5).GE.0)THEN
     hux = -0.002   ! WATER  ! dT/dz must be < - 0.2 K per 100 m.
  ELSE
     hux = -0.005  ! LAND    ! dT/dz must be < - 0.5 K per 100 m.
  ENDIF
  DO k=1,MAX(1,k50-1)
    IF (k == 1) then
      IF ((th(k)-ts)/(0.5*dz(k)) < hux) THEN
        superadiabatic = .true.
      ELSE
        superadiabatic = .false.
        exit
      ENDIF
    ELSE
      IF ((th(k)-th(k-1))/(0.5*(dz(k)+dz(k-1))) < hux) THEN
        superadiabatic = .true.
      ELSE
        superadiabatic = .false.
        exit
      ENDIF
    ENDIF
  ENDDO

  ! Determine the numer of updrafts/plumes in the grid column:
  ! Some of these criteria may be a little redundant but useful for bullet-proofing.
  !   (1) largest plume = 1.0 * dx.
  !   (2) Apply a scale-break, assuming no plumes with diameter larger than PBLH can exist.
  !   (3) max plume size beneath clouds deck approx = 0.5 * cloud_base.
  !   (4) add shear-dependent limit, when plume model breaks down. (taken out)
  !   (5) land-only limit to reduce plume sizes in weakly forced conditions
  ! Criteria (1)
  NUP2 = max(1,min(NUP,INT(dx*dcut/dl)))
  ! Criteria (2) and (4)
  !wspd_pbl=SQRT(MAX(u(kpbl)**2 + v(kpbl)**2, 0.01))
  maxwidth = 1.0*PBLH !- MIN(15.*MAX(wspd_pbl - 7.5, 0.), 0.3*PBLH)
  ! Criteria (3)
  maxwidth = MIN(maxwidth,0.5*cloud_base)
  ! Criteria (5)
  IF((landsea-1.5).LT.0)THEN
    IF (cloud_base .LT. 2000.) THEN
      width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.120)/0.03) + .5),1000.), 0.)
    ELSE
      width_flx = MAX(MIN(1000.*(0.6*tanh((flt - 0.050)/0.03) + .5),1000.), 0.)
    ENDIF
    maxwidth = MIN(maxwidth,width_flx)
  ENDIF
  ! Convert maxwidth to number of plumes
  NUP2 = MIN(MAX(INT((maxwidth - MOD(maxwidth,100.))/100), 0), NUP2)

  !Initialize values:
  ktop = 0
  maxmf= 0.0

  IF ( fltv > 0.002 .AND. NUP2 .GE. 1 .AND. superadiabatic) then
    !PRINT*," Conditions met to run mass-flux scheme",fltv,pblh

    ! Find coef C for number size density N
    cn = 0.
    d=-1.9  !set d to value suggested by Neggers 2015 (JAMES).
    do I=1,NUP !NUP2
       IF(I > NUP2) exit
       l  = dl*I                            ! diameter of plume
       cn = cn + l**d * (l*l)/(dx*dx) * dl  ! sum fractional area of each plume
    enddo
    C = Atot/cn   !Normalize C according to the defined total fraction (Atot)

    ! Find the portion of the total fraction (Atot) of each plume size:
    An2 = 0.
    do I=1,NUP !NUP2
       IF(I > NUP2) exit
       l  = dl*I                            ! diameter of plume
       N = C*l**d                           ! number density of plume n
       UPA(1,I) = N*l*l/(dx*dx) * dl        ! fractional area of plume n
       ! Make updraft area (UPA) a function of the buoyancy flux
       acfac = .5*tanh((fltv - 0.03)/0.09) + .5
       UPA(1,I)=UPA(1,I)*acfac
       An2 = An2 + UPA(1,I)                 ! total fractional area of all plumes
       !print*," plume size=",l,"; area=",An,"; total=",An2
    end do

    ! set initial conditions for updrafts
    z0=50.
    pwmin=0.1       ! was 0.5
    pwmax=0.5       ! was 3.0

    wstar=max(1.E-2,(g/thv(1)*fltv*pblh)**(1./3.))
    qstar=max(flq,1.0E-5)/wstar
    thstar=flt/wstar

    IF((landsea-1.5).GE.0)THEN
       csigma = 1.34   ! WATER
    ELSE
       csigma = 1.34   ! LAND
    ENDIF
    sigmaW =1.34*wstar*(z0/pblh)**(1./3.)*(1 - 0.8*z0/pblh)
    sigmaQT=csigma*qstar*(z0/pblh)**(-1./3.)
    sigmaTH=csigma*thstar*(z0/pblh)**(-1./3.)

    wmin=sigmaW*pwmin
    wmax=sigmaW*pwmax

    !recompute acfac for plume excess
    acfac = .5*tanh((fltv - 0.08)/0.07) + .5

    !SPECIFY SURFACE UPDRAFT PROPERTIES
    DO I=1,NUP !NUP2
       IF(I > NUP2) exit
       wlv=wmin+(wmax-wmin)/NUP2*(i-1)
       wtv=wmin+(wmax-wmin)/NUP2*i

       !SURFACE UPDRAFT VERTICAL VELOCITY
       UPW(1,I)=wmin + REAL(i)/REAL(NUP)*(wmax-wmin)

       UPU(1,I)=U(1)
       UPV(1,I)=V(1)
       UPQC(1,I)=0
       UPQT(1,I) = 0.
       UPTHV(1,I)= 0.
       UPTHL(1,I)= 0.
       k50=1 !for now, keep at lowest model layer...
       DO k=1,k50
         UPQT(1,I) = UPQT(1,I) +QT(k) +0.58*UPW(1,I)*sigmaQT/sigmaW *acfac
         UPTHV(1,I)= UPTHV(1,I)+THV(k)+0.58*UPW(1,I)*sigmaTH/sigmaW *acfac
         UPTHL(1,I)= UPTHL(1,I)+THL(k)+0.58*UPW(1,I)*sigmaTH/sigmaW *acfac
       ENDDO
       UPQT(1,I) = UPQT(1,I)/REAL(k50)
       UPTHV(1,I)= UPTHV(1,I)/REAL(k50)
       UPTHL(1,I)= UPTHL(1,I)/REAL(k50)             ! now, if the lowest layer is saturated, it will be counted for.
       UPQKE(1,I)= QKE(1)
#if (WRF_CHEM == 1)
       do ic = 1,nchem
          UPCHEM(1,I,ic)= CHEM(1,ic)
       enddo
#endif

    ENDDO

  EntThrottle = 0.001
  ! do integration  updraft
    DO I=1,NUP !NUP2
       IF(I > NUP2) exit
       QCn = 0.
       overshoot = 0
       l  = dl*I                            ! diameter of plume
       DO k=KTS+1,KTE
          !w-dependency for entrainment a la Tian and Kuang (2016)
! AEPS in Angevine et al. paper for sensitivity tests is the leading coeff
! of this equation for ENT
          ENT(k,i) = 0.5/(MIN(MAX(UPW(K-1,I),0.75),1.5)*l)
          !JOE - implement minimum background entrainment 
          ENT(k,i) = max(ENT(k,i),0.0003)
          !JOE - increase entrainment for plumes extending very high.
          IF(ZW(k) >= pblh+1500.) ENT(k,i) =ENT(k,i) + (ZW(k)-(pblh+1500.)) * 5.0E-6
          IF(UPW(K-1,I) > 2.0) ENT(k,i) = ENT(k,i) + EntThrottle*(UPW(K-1,I) - 2.0)
          ENT(k,i) = min(ENT(k,i),0.9/(ZW(k)-ZW(k-1)))

          ! Linear entrainment:
          EntExp= ENT(K,I)*(ZW(k)-ZW(k-1))
          QTn =UPQT(k-1,I) *(1.-EntExp) + QT(k-1)*EntExp
          THLn=UPTHL(k-1,I)*(1.-EntExp) + THL(k-1)*EntExp
          Un  =UPU(k-1,I)  *(1.-EntExp) + U(k-1)*EntExp
          Vn  =UPV(k-1,I)  *(1.-EntExp) + V(k-1)*EntExp
          QKEn=UPQKE(k-1,I)*(1.-EntExp) + QKE(k-1)*EntExp

#if (WRF_CHEM == 1)
          do ic = 1,nchem
             ! Exponential Entrainment:
             !chemn(ic) = chem(k,ic)*(1-EntExp)+UPCHEM(K-1,I,ic)*EntExp
             ! Linear entrainment:
             chemn(ic)=UPCHEM(k-1,I,ic)*(1.-EntExp) + chem(k-1,ic)*EntExp
         enddo
#endif

          ! get thvn,qcn
          call condensation_edmf(QTn,THLn,(P(K)+P(K-1))/2.,ZW(k),THVn,QCn)

          B=g*(0.5*(THVn+UPTHV(k-1,I))/THV(k-1) - 1.0)
! BCOEFF below is perturbed for sensitivity tests in Angevine et al. paper
          IF(B>0.)THEN
            BCOEFF = 0.15        !0.11 gives maxmf similar to Kay's scheme, w typically stays < 2.0, so doesn't hit the limits nearly as much ! orig
          ELSE
            BCOEFF = 0.2 !0.33
          ENDIF

          !WA: TEMF form. 
          !JOE-add dz max of 250, so accelerations aren't unphysically large for coarse 
          !vertical grid spacing.
          IF (UPW(K-1,I) < 0.2 ) THEN
             Wn = UPW(K-1,I) + (-2. * ENT(K,I) * UPW(K-1,I) + BCOEFF*B / MAX(UPW(K-1,I),0.2)) * MIN(ZW(k)-ZW(k-1), 250.)
          ELSE
             Wn = UPW(K-1,I) + (-2. * ENT(K,I) * UPW(K-1,I) + BCOEFF*B / UPW(K-1,I)) * MIN(ZW(k)-ZW(k-1), 250.)
          ENDIF
          !Do not allow a parcel to accelerate more than 1.25 m/s over 200 m. 
          !Add max increase of 2.0 m/s for coarse vertical resolution.
          IF(Wn > UPW(K-1,I) + MIN(1.25*(ZW(k)-ZW(k-1))/200., 2.0) ) THEN
             Wn = UPW(K-1,I) + MIN(1.25*(ZW(k)-ZW(k-1))/200., 2.0)
          ENDIF
          Wn = MIN(MAX(Wn,0.0), 3.0) ! orig

          IF (debug_mf == 1) THEN
            IF (Wn .GE. 3.0) THEN
              ! surface values
              print *," **** SUSPICIOUSLY LARGE W:"
              print *,' QCn:',QCn,' ENT=',ENT(k,i),' Nup2=',Nup2
              print *,'pblh:',pblh,' Wn:',Wn,' UPW(k-1)=',UPW(K-1,I)
              print *,'K=',k,' B=',B,' dz=',ZW(k)-ZW(k-1) 
            ENDIF
          ENDIF

!          !Allow strongly forced plumes to overshoot if KE is sufficient
          IF (fltv > 0.05 .AND. Wn <= 0 .AND. overshoot == 0) THEN
             overshoot = 1
             IF ( THV(k)-THV(k-1) .GT. 0.0 ) THEN
                bvf = SQRT( gtr*(THV(k)-THV(k-1))/(0.5*(dz(k)+dz(k-1))) )
                !vertical Froude number
                Frz = UPW(K-1,I)/(bvf*0.5*(dz(k)+dz(k-1)))
                IF ( Frz >= 0.5 ) Wn =  MIN(Frz,1.0)*UPW(K-1,I)
             ENDIF
          ELSEIF (fltv > 0.05 .AND. overshoot == 1) THEN
             !Do not let overshooting parcel go more than 1 layer up
             Wn = 0.0
          ENDIF

          !Limit very tall plumes
          Wn=Wn*EXP(-MAX(ZW(k) - MIN(pblh+2000.,3000.),0.0)/1000.)
          IF(ZW(k) >= MIN(pblh+3000.,4000.))Wn=0.

          !JOE- minimize the plume penetratration in stratocu-topped PBL
          IF (fltv < 0.06) THEN
             IF(ZW(k) >= pblh-200. .AND. qc(k) > 1e-5 .AND. I > 4) Wn=0.
          ENDIF

          IF (Wn > 0.) THEN
             UPW(K,I)=Wn  !Wn !sqrt(Wn2)
             UPTHV(K,I)=THVn
             UPTHL(K,I)=THLn
             UPQT(K,I)=QTn
             UPQC(K,I)=QCn
             UPU(K,I)=Un
             UPV(K,I)=Vn
             UPQKE(K,I)=QKEn
             UPA(K,I)=UPA(K-1,I)
#if (WRF_CHEM == 1)
             do ic = 1,nchem
                UPCHEM(k,I,ic) = chemn(ic)
             enddo
#endif
             ktop = MAX(ktop,k)
          ELSE
             exit  !exit k-loop
          END IF
       ENDDO
       IF (debug_mf == 1) THEN
          IF (MAXVAL(UPW(:,I)) .GE. 3. .OR. MINVAL(UPA(:,I)) < 0.0 .OR. &
              MAXVAL(UPA(:,I)) > Atot .OR. NUP2 > 10) THEN
             ! surface values
             print *,'flq:',flq,' fltv:',fltv,' Nup2=',Nup2
             print *,'pblh:',pblh,' wstar:',wstar,' ktop=',ktop
             print *,'sigmaW=',sigmaW,' sigmaTH=',sigmaTH,' sigmaQT=',sigmaQT
             ! means
             print *,'u:',u
             print *,'v:',v
             print *,'thl:',thl
             print *,'UPA:',UPA(:,I)
             print *,'UPW:',UPW(:,I)
             print *,'UPTHL:',UPTHL(:,I)
             print *,'UPQT:',UPQT(:,I)
             print *,'ENT:',ENT(:,I)
          ENDIF
       ENDIF
    ENDDO
  ELSE
    !At least one of the conditions was not met for activating the MF scheme.
    NUP2=0. 
  END IF !end criteria for mass-flux scheme

  ktop=MIN(ktop,KTE-1)  !  Just to be safe...
!        
!  get updraft properties, for saving
!
  IF(nup2 > 0) THEN

    DO k=KTS,KTE-1
      IF(k > KTOP) exit
      DO I=1,NUP !NUP2
        IF(I > NUP2) exit
        edmf_a(K)=edmf_a(K)+UPA(K+1,I)
        edmf_w(K)=edmf_w(K)+UPA(K+1,I)*UPW(K+1,I)
        edmf_qt(K)=edmf_qt(K)+UPA(K+1,I)*UPQT(K+1,I)
        edmf_thl(K)=edmf_thl(K)+UPA(K+1,I)*UPTHL(K+1,I)
        edmf_ent(K)=edmf_ent(K)+UPA(K+1,I)*ENT(K+1,I)
        edmf_qc(K)=edmf_qc(K)+UPA(K+1,I)*UPQC(K+1,I)
#if (WRF_CHEM == 1)
        do ic = 1,nchem
          edmf_chem(k,ic) = edmf_chem(k,ic) + UPA(K+1,I)*UPCHEM(k,I,ic)
        enddo
#endif
      ENDDO

      IF (edmf_a(k)>0.) THEN
        edmf_w(k)=edmf_w(k)/edmf_a(k)
        edmf_qt(k)=edmf_qt(k)/edmf_a(k)
        edmf_thl(k)=edmf_thl(k)/edmf_a(k)
        edmf_ent(k)=edmf_ent(k)/edmf_a(k)
        edmf_qc(k)=edmf_qc(k)/edmf_a(k)
#if (WRF_CHEM == 1)
        do ic = 1,nchem
          edmf_chem(k,ic) = edmf_chem(k,ic)/edmf_a(k)
        enddo
#endif
        edmf_a(k)=edmf_a(k)*Psig_w

        !FIND MAXIMUM MASS-FLUX IN THE COLUMN:
        IF(edmf_a(k)*edmf_w(k) > maxmf) maxmf = edmf_a(k)*edmf_w(k)
      ENDIF
    ENDDO

    DO k=KTS,KTE
      IF(k > KTOP) exit
      DO I=1,NUP !NUP2
        IF(I > NUP2) exit
        s_aw(k)   = s_aw(K)    + UPA(K,I)*UPW(K,I)*Psig_w * (1.0+rstoch_col(k))
        s_awthl(k)= s_awthl(K) + UPA(K,i)*UPW(K,I)*UPTHL(K,I)*Psig_w * (1.0+rstoch_col(k))
        s_awqt(k) = s_awqt(K)  + UPA(K,i)*UPW(K,I)*UPQT(K,I)*Psig_w * (1.0+rstoch_col(k))
        s_awqc(k) = s_awqc(K)  + UPA(K,i)*UPW(K,I)*UPQC(K,I)*Psig_w * (1.0+rstoch_col(k))
        IF (momentum_opt > 0) THEN
          s_awu(k)  = s_awu(K)   + UPA(K,i)*UPW(K,I)*UPU(K,I)*Psig_w * (1.0+rstoch_col(k))
          s_awv(k)  = s_awv(K)   + UPA(K,i)*UPW(K,I)*UPV(K,I)*Psig_w * (1.0+rstoch_col(k))
        ENDIF
        IF (tke_opt > 0) THEN
          s_awqke(k)= s_awqke(K) + UPA(K,i)*UPW(K,I)*UPQKE(K,I)*Psig_w * (1.0+rstoch_col(k))
        ENDIF
#if (WRF_CHEM == 1)
        do ic = 1,nchem
          s_awchem(k,ic) = s_awchem(k,ic) + UPA(K,i)*UPW(K,I)*UPCHEM(K,I,ic)*Psig_w * (1.0+rstoch_col(k))
        enddo
#endif
      ENDDO
      s_awqv(k) = s_awqt(k)  - s_awqc(k)
    ENDDO

!JOE: ADD CLDFRA_bl1d, qc_bl1d. Note that they have already been defined in
!     mym_condensation. Here, a shallow-cu component is added.  
     DO K=KTS,KTE
        IF(k > KTOP) exit
        IF(edmf_qc(k)>0.0)THEN
            satvp = 3.80*exp(17.27*(th(k)-273.)/ &
                   (th(k)-36.))/(.01*p(k))
            rhgrid = max(.01,MIN( 1., qv(k) /satvp))

            !COMPUTE CLDFRA & QC_BL FROM MASS-FLUX SCHEME and recompute vt & vq

            xl = xl_blend(tk(k))                ! obtain blended heat capacity 
            tlk = thl(k)*(p(k)/p1000mb)**rcp    ! recover liquid temp (tl) from thl
            qsat_tl = qsat_blend(tlk,p(k))      ! get saturation water vapor mixing ratio
                                                !   at tl and p
            rsl = xl*qsat_tl / (r_v*tlk**2)     ! slope of C-C curve at t = tl
                                                ! CB02, Eqn. 4
            cpm = cp + qt(k)*cpv                ! CB02, sec. 2, para. 1
            a   = 1./(1. + xl*rsl/cpm)          ! CB02 variable "a"
            b9  = a*rsl                         ! CB02 variable "b" 

            q2p  = xlvcp/exner(k)
            pt = thl(k) +q2p*edmf_qc(k) ! potential temp
            bb = b9*tk(k)/pt ! bb is "b9" in BCMT95.  Their "b9" differs from
                           ! "b9" in CB02 by a factor
                           ! of T/theta.  Strictly, b9 above is formulated in
                           ! terms of sat. mixing ratio, but bb in BCMT95 is
                           ! cast in terms of sat. specific humidity.  The
                           ! conversion is neglected here.
            qww   = 1.+0.61*qt(k)
            alpha = 0.61*pt
            t     = th(k)*exner(k)
            beta  = pt*xl/(t*cp) - 1.61*pt
            !Buoyancy flux terms have been moved to the end of this section...

            !Now calculate convective component of the cloud fraction:
            if (a > 0.0) then
               f = MIN(1.0/a, 4.0)              ! f is vertical profile scaling function (CB2005)
            else
               f = 1.0
            endif
            sigq = 9.E-3 * edmf_a(k) * edmf_w(k) * f ! convective component of sigma (CB2005)
            sigq = SQRT(sigq**2 + sgm(k)**2)    ! combined conv + stratus components

            qmq = a * (qt(k) - qsat_tl)         ! saturation deficit/excess;
                                                !   the numerator of Q1
            mf_cf = min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.02),0.6)
            IF (rhgrid >= .93) THEN
               !IN high RH, defer to stratus component if > convective component
               cldfra_bl1d(k) = MAX(mf_cf, cldfra_bl1d(k))
               IF (cldfra_bl1d(k) > edmf_a(k)) THEN
                  qc_bl1d(k) = edmf_qc(k)*edmf_a(k)/cldfra_bl1d(k)
               ELSE
                 cldfra_bl1d(k)=edmf_a(k)
                 qc_bl1d(k) = edmf_qc(k)
               ENDIF
            ELSE
               IF (mf_cf > edmf_a(k)) THEN
                  cldfra_bl1d(k) = mf_cf
                  qc_bl1d(k) = edmf_qc(k)*edmf_a(k)/mf_cf
               ELSE
                  cldfra_bl1d(k)=edmf_a(k)
                  qc_bl1d(k) = edmf_qc(k)
               ENDIF
            ENDIF
            !Now recalculate the terms for the buoyancy flux for mass-flux clouds:
            !See mym_condensation for details on these formulations.  The
            !cloud-fraction bounding was added to improve cloud retention,
            !following RAP and HRRR testing.
            Fng = 2.05 ! the non-Gaussian transport factor (assumed constant)
            vt(k) = qww   - MIN(0.20,cldfra_bl1D(k))*beta*bb*Fng - 1.
            vq(k) = alpha + MIN(0.20,cldfra_bl1D(k))*beta*a*Fng  - tv0
         ENDIF

      ENDDO

    ENDIF  !end nup2 > 0

    !modify output (negative: dry plume, positive: moist plume)
    IF (ktop > 0) THEN
      maxqc = maxval(edmf_qc(1:ktop)) 
      IF ( maxqc < 1.E-8) maxmf = -1.0*maxmf
    ENDIF
       
IF (edmf_w(1) > 4.0) THEN 
! surface values
    print *,'flq:',flq,' fltv:',fltv
    print *,'pblh:',pblh,' wstar:',wstar
    print *,'sigmaW=',sigmaW,' sigmaTH=',sigmaTH,' sigmaQT=',sigmaQT
 
! mean updrafts
   print *,' edmf_a',edmf_a(1:14)
   print *,' edmf_w',edmf_w(1:14)
   print *,' edmf_qt:',edmf_qt(1:14)
   print *,' edmf_thl:',edmf_thl(1:14)
ENDIF !END Debugging


#ifdef HARDCODE_VERTICAL
# undef kts
# undef kte
#endif

END SUBROUTINE MF1_MF