
module modal_aero_convproc
!---------------------------------------------------------------------------------
! Purpose:
!
! CAM interface to aerosol/trace-gas convective cloud processing scheme
!
! currently these routines assume stratiform and convective clouds only interact 
! through the detrainment of convective cloudborne material into stratiform clouds
!
! thus the stratiform-cloudborne aerosols (in the qqcw array) are not processed 
! by the convective up/downdrafts, but are affected by the detrainment
!
! Author: R. C. Easter
!
!---------------------------------------------------------------------------------

   use shr_kind_mod, only: r8=>shr_kind_r8
   use physconst,    only: gravit                              
   use ppgrid,       only: pver, pcols, pverp, begchunk, endchunk
   use cam_history,  only: outfld, addfld, horiz_only, add_default
   use cam_logfile,  only: iulog
   use cam_abortutils, only: endrun
   use physconst,    only: spec_class_aerosol, spec_class_gas

   implicit none

   save
   private                         ! Make default type private to the module

! Public methods

   public :: &
      ma_convproc_register,         &!
      ma_convproc_init,             &!
      ma_convproc_intr               !


!
! module data
!
   !deepconv_wetdep_history variable controls history output of additonal deep-convection wet deposition fields 
   !(in addition to the normal fields for total-convection wet deposition)
   logical, parameter, public :: deepconv_wetdep_history = .true.
   logical, parameter :: use_cwaer_for_activate_maxsat = .false.
   logical, parameter :: apply_convproc_tend_to_ptend = .true.

   real(r8) :: hund_ovr_g ! = 100.0_r8/gravit
!  used with zm_conv mass fluxes and delta-p
!     for mu = [mbar/s],   mu*hund_ovr_g = [kg/m2/s]
!     for dp = [mbar] and q = [kg/kg],   q*dp*hund_ovr_g = [kg/m2]

!  method1_activate_nlayers = number of layers (including cloud base) where activation is applied
   integer, parameter  :: method1_activate_nlayers = 2
!  method2_activate_smaxmax = the uniform or peak supersat value (as 0-1 fraction = percent*0.01)
   real(r8), parameter :: method2_activate_smaxmax = 0.003_r8

!  method_reduce_actfrac = 1 -- multiply activation fractions by factor_reduce_actfrac
!                               (this works ok with convproc_method_activate = 1 but not for ... = 2)
!                        = 2 -- do 2 iterations to get an overall reduction by factor_reduce_actfrac
!                               (this works ok with convproc_method_activate = 1 or 2)
!                        = other -- do nothing involving reduce_actfrac
   integer, parameter  :: method_reduce_actfrac = 0
   real(r8), parameter :: factor_reduce_actfrac = 0.5_r8

!
! Private module data
!

   logical, private :: convproc_do_gas, convproc_do_aer
   logical, private :: convproc_prevap_resusp_fixaa = .false. ! REASTER 08/05/2015
   !  convproc_method_fixaa - see explanation in subr. ma_convproc_tend(                                           &
   integer, private :: convproc_method_activate
   !  convproc_method_activate - 1=apply abdulrazzak-ghan to entrained aerosols for lowest nlayers
   !                             2=do secondary activation with prescribed supersat

!=========================================================================================
  contains


!=========================================================================================
subroutine ma_convproc_register

!----------------------------------------
! Purpose: register fields with the physics buffer
!----------------------------------------

  use physics_buffer, only:  pbuf_add_field

  implicit none

  integer idx

end subroutine ma_convproc_register



!=========================================================================================
subroutine ma_convproc_init

!----------------------------------------
! Purpose:  declare output fields, initialize variables needed by convection
!----------------------------------------

  use cam_history,    only: outfld, addfld, horiz_only, add_default
  use physics_buffer, only: pbuf_add_field
  use phys_control,   only: phys_getopts
  use ppgrid,         only: pcols, pver
  use spmd_utils,     only: masterproc
  use error_messages, only: alloc_err	

  implicit none

  integer :: npass_calc_updraft
  logical :: history_aerosol      ! Output the MAM aerosol tendencies
  logical :: resus_fix
  character(len=100) :: msg

! 
! Add history fields
!
    call phys_getopts( history_aerosol_out=history_aerosol, &
        convproc_do_aer_out = convproc_do_aer, &
        convproc_do_gas_out = convproc_do_gas, &
        convproc_method_activate_out = convproc_method_activate )

    call addfld(      'SH_MFUP_MAX', horiz_only, 'A', 'kg/m2', &
                      'Shallow conv. column-max updraft mass flux' )
    call addfld(      'SH_WCLDBASE', horiz_only, 'A', 'm/s', &
                      'Shallow conv. cloudbase vertical velocity' )
    call addfld(      'SH_KCLDBASE', horiz_only, 'A', '1', &
                      'Shallow conv. cloudbase level index' )

    call addfld(      'DP_MFUP_MAX', horiz_only, 'A', 'kg/m2', &
                      'Deep conv. column-max updraft mass flux' )
    call addfld(      'DP_WCLDBASE', horiz_only, 'A', 'm/s', &
                      'Deep conv. cloudbase vertical velocity' )
    call addfld(      'DP_KCLDBASE', horiz_only, 'A', '1', &
                      'Deep conv. cloudbase level index' )

    if ( history_aerosol .and. convproc_do_aer ) then
       call add_default( 'SH_MFUP_MAX', 1, ' ' )
       call add_default( 'SH_WCLDBASE', 1, ' ' )
       call add_default( 'SH_KCLDBASE', 1, ' ' )
       call add_default( 'DP_MFUP_MAX', 1, ' ' )
       call add_default( 'DP_WCLDBASE', 1, ' ' )
       call add_default( 'DP_KCLDBASE', 1, ' ' )
    end if

!
! Print control variable settings
!
   if ( .not. masterproc ) return

   write(*,'(a,l12)')     'ma_convproc_init - convproc_do_aer               = ', &
      convproc_do_aer
   write(*,'(a,l12)')     'ma_convproc_init - convproc_do_gas               = ', &
      convproc_do_gas
   write(*,'(a,l12)')     'ma_convproc_init - use_cwaer_for_activate_maxsat = ', &
      use_cwaer_for_activate_maxsat
   write(*,'(a,l12)')     'ma_convproc_init - apply_convproc_tend_to_ptend  = ', &
      apply_convproc_tend_to_ptend
   write(*,'(a,i12)')     'ma_convproc_init - convproc_method_activate      = ', &
      convproc_method_activate
   write(*,'(a,i12)')     'ma_convproc_init - method1_activate_nlayers      = ', &
      method1_activate_nlayers
   write(*,'(a,1pe12.4)') 'ma_convproc_init - method2_activate_smaxmax      = ', &
      method2_activate_smaxmax
   write(*,'(a,i12)')     'ma_convproc_init - method_reduce_actfrac         = ', &
      method_reduce_actfrac
   write(*,'(a,1pe12.4)') 'ma_convproc_init - factor_reduce_actfrac         = ', &
      factor_reduce_actfrac

   npass_calc_updraft = 1
   if ( (method_reduce_actfrac == 2)      .and. &
        (factor_reduce_actfrac >= 0.0_r8) .and. &
        (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2
   write(*,'(a,i12)')     'ma_convproc_init - npass_calc_updraft            = ', &
      npass_calc_updraft

   return
end subroutine ma_convproc_init



!=========================================================================================
subroutine ma_convproc_intr( state, ptend, pbuf, ztodt,             &
                           dp_frac, icwmrdp, rprddp, evapcdp,       &
                           sh_frac, icwmrsh, rprdsh, evapcsh,       &
                           dlf, dlfsh, cmfmcsh, sh_e_ed_ratio,      &
                           nsrflx_mzaer2cnvpr, qsrflx_mzaer2cnvpr,  &
                           aerdepwetis,                             &
                           mu, md, du, eu,                          &
                           ed, dp, dsubcld,                         &
                           jt, maxg, ideep, lengath, species_class, &
                           mam_prevap_resusp_optaa,                 &
                           history_aero_prevap_resusp               )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Convective cloud processing (transport, activation/resuspension,
!    wet removal) of aerosols and trace gases.
!    (Currently no aqueous chemistry and no trace-gas wet removal)
! Does aerosols    when convproc_do_aer is .true.
! Does trace gases when convproc_do_gas is .true.
!
! Does deep and shallow convection
! Uses mass fluxes, cloud water, precip production from the
!    convective cloud routines
! 
! Author: R. Easter
! 
!-----------------------------------------------------------------------

   use physics_types, only: physics_state, physics_ptend, physics_ptend_init
   use time_manager,  only: get_nstep
   use physics_buffer, only: physics_buffer_desc, pbuf_get_index
   use constituents,  only: pcnst, cnst_name
   use error_messages, only: alloc_err	

   use modal_aero_data, only: lmassptr_amode, nspec_amode, ntot_amode, numptr_amode
 
! Arguments
   type(physics_state), intent(in ) :: state          ! Physics state variables
   type(physics_ptend), intent(inout) :: ptend          ! indivdual parameterization tendencies
   type(physics_buffer_desc), pointer :: pbuf(:)
   real(r8), intent(in) :: ztodt                          ! 2 delta t (model time increment)

   real(r8), intent(in)    :: dp_frac(pcols,pver) ! Deep conv cloud frac (0-1)
   real(r8), intent(in)    :: icwmrdp(pcols,pver) ! Deep conv cloud condensate (kg/kg - in cloud)
   real(r8), intent(in)    :: rprddp(pcols,pver)  ! Deep conv precip production (kg/kg/s - grid avg)
   real(r8), intent(in)    :: evapcdp(pcols,pver) ! Deep conv precip evaporation (kg/kg/s - grid avg)
   real(r8), intent(in)    :: sh_frac(pcols,pver) ! Shal conv cloud frac (0-1)
   real(r8), intent(in)    :: icwmrsh(pcols,pver) ! Shal conv cloud condensate (kg/kg - in cloud)
   real(r8), intent(in)    :: rprdsh(pcols,pver)  ! Shal conv precip production (kg/kg/s - grid avg)
   real(r8), intent(in)    :: evapcsh(pcols,pver) ! Shal conv precip evaporation (kg/kg/s - grid avg)
   real(r8), intent(in)    :: dlf(pcols,pver)     ! Tot  conv cldwtr detrainment (kg/kg/s - grid avg)
   real(r8), intent(in)    :: dlfsh(pcols,pver)   ! Shal conv cldwtr detrainment (kg/kg/s - grid avg)
   real(r8), intent(in)    :: cmfmcsh(pcols,pverp) ! Shal conv mass flux (kg/m2/s)
   real(r8), intent(in)    :: sh_e_ed_ratio(pcols,pver)  ! shallow conv [ent/(ent+det)] ratio
   integer,  intent(in)    :: nsrflx_mzaer2cnvpr
   real(r8), intent(in)    :: qsrflx_mzaer2cnvpr(pcols,pcnst,nsrflx_mzaer2cnvpr)
   real(r8), intent(inout) :: aerdepwetis(pcols,pcnst)  ! aerosol wet deposition (interstitial)

                                               ! mu, md, ..., ideep, lengath are all deep conv variables
                                               ! *** AND ARE GATHERED ***
   real(r8), intent(in)    :: mu(pcols,pver)   ! Updraft mass flux (positive)
   real(r8), intent(in)    :: md(pcols,pver)   ! Downdraft mass flux (negative)
   real(r8), intent(in)    :: du(pcols,pver)   ! Mass detrain rate from updraft
   real(r8), intent(in)    :: eu(pcols,pver)   ! Mass entrain rate into updraft
   real(r8), intent(in)    :: ed(pcols,pver)   ! Mass entrain rate into downdraft
                           ! eu, ed, du are "d(massflux)/dp" and are all positive
   real(r8), intent(in)    :: dp(pcols,pver)   ! Delta pressure between interfaces
   real(r8), intent(in)    :: dsubcld(pcols)   ! Delta pressure from cloud base to sfc

   integer,  intent(in)    :: jt(pcols)         ! Index of cloud top for each column
   integer,  intent(in)    :: maxg(pcols)       ! Index of cloud top for each column
   integer,  intent(in)    :: ideep(pcols)      ! Gathering array
   integer,  intent(in)    :: lengath           ! Gathered min lon indices over which to operate
   integer,  intent(in)    :: species_class(:)
   integer,  intent(in)    :: mam_prevap_resusp_optaa
   logical,  intent(in)    :: history_aero_prevap_resusp


! Local variables
   integer, parameter :: nsrflx = 6        ! last dimension of qsrflx ! REASTER 08/05/2015
   integer  :: i, ii, itmpa
   integer  :: k
   integer  :: l, ll, lchnk
   integer  :: n, ncol, nstep

   real(r8) :: dlfdp(pcols,pver)
   real(r8) :: dpdry(pcols,pver)
   real(r8) :: dqdt(pcols,pver,pcnst)
   real(r8) :: dt
   real(r8) :: qa(pcols,pver,pcnst), qb(pcols,pver,pcnst)
   real(r8) :: qsrflx(pcols,pcnst,nsrflx)
   real(r8) :: sflxic(pcols,pcnst)
   real(r8) :: sflxid(pcols,pcnst)
   real(r8) :: sflxec(pcols,pcnst)
   real(r8) :: sflxed(pcols,pcnst)
   real(r8) :: tmpa, tmpb, tmpg

   logical  :: dotend(pcnst)

! physics buffer fields 
   integer itim, ifld
   real(r8), pointer, dimension(:,:,:) :: fracis  ! fraction of transported species that are insoluble

!
! Initialize
!

! apply this minor fix when doing resuspend to coarse mode
   if (mam_prevap_resusp_optaa >= 30) convproc_prevap_resusp_fixaa = .true.

   lchnk = state%lchnk
   ncol  = state%ncol
   nstep = get_nstep()
   dt = ztodt

   hund_ovr_g = 100.0_r8/gravit
!  used with zm_conv mass fluxes and delta-p
!     for mu = [mbar/s],   mu*hund_ovr_g = [kg/m2/s]
!     for dp = [mbar] and q = [kg/kg],   q*dp*hund_ovr_g = [kg/m2]

   sflxic(:,:) = 0.0_r8
   sflxid(:,:) = 0.0_r8
   sflxec(:,:) = 0.0_r8
   sflxed(:,:) = 0.0_r8
   do l = 1, pcnst
      if ( (species_class(l) == spec_class_aerosol) .and. ptend%lq(l) ) then
         sflxec(1:ncol,l) = qsrflx_mzaer2cnvpr(1:ncol,l,1) 
         sflxed(1:ncol,l) = qsrflx_mzaer2cnvpr(1:ncol,l,2) 
      end if
   end do

!
! Associate pointers with physics buffer fields
!
!  ifld = pbuf_get_fld_idx('FRACIS')
!  fracis  => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,state%lchnk,1:pcnst)


!
! prepare for deep conv processing
!
  do l = 1, pcnst
     if ( ptend%lq(l) ) then
        ! calc new q (after calcaersize and mz_aero_wet_intr)
        qa(1:ncol,:,l) = state%q(1:ncol,:,l) + dt*ptend%q(1:ncol,:,l)
        qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) ) 

!    skip this -- if code generates negative q, 
!    then you need to see the messages and fix it
!       if ( apply_convproc_tend_to_ptend ) then
!          ! adjust ptend%q when qa < 0.0
!          ptend%q(1:ncol,:,l) = ptend%q(1:ncol,:,l) &
!             + ( qb(1:ncol,:,l) - qa(1:ncol,:,l) )/dt
!       end if

     else
        ! use old q
        qb(1:ncol,:,l) = state%q(1:ncol,:,l)
     end if
  end do
  dqdt(:,:,:) = 0.0_r8
  qsrflx(:,:,:) = 0.0_r8
  dotend(:) = .false.


!
! do deep conv processing
!
  if (convproc_do_aer .or. convproc_do_gas) then

  dlfdp(1:ncol,:) = max( (dlf(1:ncol,:) - dlfsh(1:ncol,:)), 0.0_r8 )

  call ma_convproc_dp_intr(                    &
     state, pbuf, dt,                          &
     dp_frac, icwmrdp, rprddp, evapcdp, dlfdp, &
     mu, md, du, eu,                           &
     ed, dp, dsubcld,                          &
     jt, maxg, ideep, lengath,                 &
     qb, dqdt, dotend, nsrflx, qsrflx,         &
     species_class, mam_prevap_resusp_optaa    )


! apply deep conv processing tendency and prepare for shallow conv processing
  do l = 1, pcnst
     if ( .not. dotend(l) ) cycle

     ! calc new q (after ma_convproc_dp_intr)
     qa(1:ncol,:,l) = qb(1:ncol,:,l) + dt*dqdt(1:ncol,:,l)
     qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) ) 

! skip this -- if code generates negative q, 
! then you need to see the messages and fix it
!    ! adjust dqdt when qa < 0.0
!    dqdt(1:ncol,:,l) = dqdt(1:ncol,:,l) &
!       + ( qb(1:ncol,:,l) - qa(1:ncol,:,l) )/dt

     if ( apply_convproc_tend_to_ptend ) then
        ! add dqdt onto ptend%q and set ptend%lq
        ptend%q(1:ncol,:,l) = ptend%q(1:ncol,:,l) + dqdt(1:ncol,:,l)
        ptend%lq(l) = .true.
     end if

     if ((species_class(l) == spec_class_aerosol) .or. &
         (species_class(l) == spec_class_gas    )) then
        ! these used for history file wetdep diagnostics
        sflxic(1:ncol,l) = sflxic(1:ncol,l) + qsrflx(1:ncol,l,4) 
        sflxid(1:ncol,l) = sflxid(1:ncol,l) + qsrflx(1:ncol,l,4) 
        sflxec(1:ncol,l) = sflxec(1:ncol,l) + qsrflx(1:ncol,l,6) ! REASTER 08/05/2015
        sflxed(1:ncol,l) = sflxed(1:ncol,l) + qsrflx(1:ncol,l,6) ! REASTER 08/05/2015
     end if

     if (species_class(l) == spec_class_aerosol) then
        ! this used for surface coupling
        aerdepwetis(1:ncol,l) = aerdepwetis(1:ncol,l) &
           + qsrflx(1:ncol,l,4) + qsrflx(1:ncol,l,5) 
     end if

  end do ! l

  dqdt(:,:,:) = 0.0_r8
  qsrflx(:,:,:) = 0.0_r8
  dotend(:) = .false.

  end if ! (convproc_do_aer  .or. convproc_do_gas ) then


!
! do shallow conv processing
!
  if (convproc_do_aer .or. convproc_do_gas ) then

  call ma_convproc_sh_intr(                    &
     state, pbuf, dt,                          &
     sh_frac, icwmrsh, rprdsh, evapcsh, dlfsh, &
     cmfmcsh, sh_e_ed_ratio,                   &
     qb, dqdt, dotend, nsrflx, qsrflx,         &
     species_class, mam_prevap_resusp_optaa    )


! apply shallow conv processing tendency
  do l = 1, pcnst
     if ( .not. dotend(l) ) cycle

     ! calc new q (after ma_convproc_sh_intr)
     qa(1:ncol,:,l) = qb(1:ncol,:,l) + dt*dqdt(1:ncol,:,l)
     qb(1:ncol,:,l) = max( 0.0_r8, qa(1:ncol,:,l) ) 

! skip this -- if code generates negative q, 
! then you need to see the messages and fix it
!    ! adjust dqdt when qa < 0.0
!    dqdt(1:ncol,:,l) = dqdt(1:ncol,:,l) &
!       + ( qb(1:ncol,:,l) - qa(1:ncol,:,l) )/dt

     if ( apply_convproc_tend_to_ptend ) then
        ! add dqdt onto ptend%q and set ptend%lq
        ptend%q(1:ncol,:,l) = ptend%q(1:ncol,:,l) + dqdt(1:ncol,:,l)
        ptend%lq(l) = .true.
     end if

     if ((species_class(l) == spec_class_aerosol) .or. &
         (species_class(l) == spec_class_gas    )) then
        sflxic(1:ncol,l) = sflxic(1:ncol,l) + qsrflx(1:ncol,l,4) 
        sflxec(1:ncol,l) = sflxec(1:ncol,l) + qsrflx(1:ncol,l,6) ! REASTER 08/05/2015
     end if

     if (species_class(l) == spec_class_aerosol) then
        aerdepwetis(1:ncol,l) = aerdepwetis(1:ncol,l) &
           + qsrflx(1:ncol,l,4) + qsrflx(1:ncol,l,5) 
     end if

  end do ! l

  end if ! (convproc_do_aer  .or. convproc_do_gas) then


! output wet deposition fields to history
!    I = in-cloud removal;     E = precip-evap resuspension
!    C = convective (total);   D = deep convective
! note that the precip-evap resuspension includes that resulting from
!    below-cloud removal, calculated in mz_aero_wet_intr
  if (convproc_do_aer .and. apply_convproc_tend_to_ptend ) then
     do n = 1, ntot_amode
     do ll = 0, nspec_amode(n)
        if (ll == 0) then
           l = numptr_amode(n)
        else
           l = lmassptr_amode(ll,n)
        end if

        call outfld( trim(cnst_name(l))//'SFWET', aerdepwetis(:,l), pcols, lchnk )
        call outfld( trim(cnst_name(l))//'SFSIC', sflxic(:,l), pcols, lchnk )
        if ( history_aero_prevap_resusp ) &
        call outfld( trim(cnst_name(l))//'SFSEC', sflxec(:,l), pcols, lchnk )

        if ( deepconv_wetdep_history ) then
           call outfld( trim(cnst_name(l))//'SFSID', sflxid(:,l), pcols, lchnk )
           if ( history_aero_prevap_resusp ) &
           call outfld( trim(cnst_name(l))//'SFSED', sflxed(:,l), pcols, lchnk )
        end if
     end do ! ll
     end do ! n
  end if


end subroutine ma_convproc_intr



!=========================================================================================
subroutine ma_convproc_dp_intr(                &
     state, pbuf, dt,                          &
     dp_frac, icwmrdp, rprddp, evapcdp, dlfdp, &
     mu, md, du, eu,                           &
     ed, dp, dsubcld,                          &
     jt, maxg, ideep, lengath,                 &
     q, dqdt, dotend, nsrflx, qsrflx,          &
     species_class, mam_prevap_resusp_optaa    )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Convective cloud processing (transport, activation/resuspension,
!    wet removal) of aerosols and trace gases.
!    (Currently no aqueous chemistry and no trace-gas wet removal)
! Does aerosols    when convproc_do_aer is .true.
! Does trace gases when convproc_do_gas is .true.
!
! This routine does deep convection
! Uses mass fluxes, cloud water, precip production from the
!    convective cloud routines
! 
! Author: R. Easter
! 
!-----------------------------------------------------------------------

   use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
   use time_manager,   only: get_nstep
   use physics_buffer, only: pbuf_get_index, physics_buffer_desc, pbuf_get_field
   use constituents,   only: pcnst, cnst_get_ind, cnst_name
   use error_messages, only: alloc_err	

   use physconst,      only: gravit, rair
   use phys_grid,      only: get_lat_all_p, get_lon_all_p, get_rlat_all_p, get_rlon_all_p

   use modal_aero_data, only: lmassptr_amode, nspec_amode, ntot_amode, numptr_amode
 
! Arguments
   type(physics_state), intent(in ) :: state          ! Physics state variables
   type(physics_buffer_desc), pointer :: pbuf(:)

   real(r8), intent(in) :: dt                         ! delta t (model time increment)

   real(r8), intent(in)    :: q(pcols,pver,pcnst)
   real(r8), intent(inout) :: dqdt(pcols,pver,pcnst)
   logical,  intent(inout) :: dotend(pcnst)
   integer,  intent(in)    :: nsrflx
   real(r8), intent(inout) :: qsrflx(pcols,pcnst,nsrflx)

   real(r8), intent(in)    :: dp_frac(pcols,pver) ! Deep conv cloud frac (0-1)
   real(r8), intent(in)    :: icwmrdp(pcols,pver) ! Deep conv cloud condensate (kg/kg - in cloud)
   real(r8), intent(in)    :: rprddp(pcols,pver)  ! Deep conv precip production (kg/kg/s - grid avg)
   real(r8), intent(in)    :: evapcdp(pcols,pver) ! Deep conv precip evaporation (kg/kg/s - grid avg)
   real(r8), intent(in)    :: dlfdp(pcols,pver)   ! Deep conv cldwtr detrainment (kg/kg/s - grid avg)

                                               ! mu, md, ..., ideep, lengath are all deep conv variables
   real(r8), intent(in)    :: mu(pcols,pver)   ! Updraft mass flux (positive)
   real(r8), intent(in)    :: md(pcols,pver)   ! Downdraft mass flux (negative)
   real(r8), intent(in)    :: du(pcols,pver)   ! Mass detrain rate from updraft
   real(r8), intent(in)    :: eu(pcols,pver)   ! Mass entrain rate into updraft
   real(r8), intent(in)    :: ed(pcols,pver)   ! Mass entrain rate into downdraft
                           ! eu, ed, du are "d(massflux)/dp" and are all positive
   real(r8), intent(in)    :: dp(pcols,pver)   ! Delta pressure between interfaces
   real(r8), intent(in)    :: dsubcld(pcols)   ! Delta pressure from cloud base to sfc

   integer,  intent(in)    :: jt(pcols)         ! Index of cloud top for each column
   integer,  intent(in)    :: maxg(pcols)       ! Index of cloud top for each column
   integer,  intent(in)    :: ideep(pcols)      ! Gathering array
   integer,  intent(in)    :: lengath           ! Gathered min lon indices over which to operate
   integer,  intent(in)    :: species_class(:)
   integer,  intent(in)    :: mam_prevap_resusp_optaa

!  real(r8), intent(in)    :: concld(pcols,pver) ! Convective cloud cover

! Local variables
   integer, parameter :: idiag_aa = -1  ! controls diagnostic output at 2 selected grid points
   integer, parameter :: idiag_gg = -1  ! controls using special profiles for gases at these grid points

   integer :: i, ii, itmpa
   integer :: ixcldice, ixcldliq              ! constituent indices for cloud liquid and ice water.
   integer :: ixh2o2, ixbc_a1, ixbc_a5, ixso4_a1
   integer :: itmpveca(pcols), itmpvecb(pcols)
   integer :: k, kaa, kbb, kk
   integer :: l, ll, lchnk, lun
   integer :: lat_ndx(pcols), lon_ndx(pcols)
   integer :: n, ncol, nstep

   real(r8) :: dpdry(pcols,pver)     ! layer delta-p-dry (mb)
   real(r8) :: fracice(pcols,pver)   ! Ice fraction of cloud droplets
   real(r8) :: lat_deg(pcols), lon_deg(pcols)
   real(r8) :: qaa(pcols,pver,pcnst), qbb(pcols,pver,pcnst)
   real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg
   real(r8) :: tmpveca(300), tmpvecb(300), tmpvecc(300)
   real(r8) :: xx_mfup_max(pcols), xx_wcldbase(pcols), xx_kcldbase(pcols)

! physics buffer fields 
   integer itim, ifld
   real(r8), pointer, dimension(:,:,:) :: fracis  ! fraction of transported species that are insoluble

!
! Initialize
!
   lun = iulog

! call physics_ptend_init(ptend)

!
! Associate pointers with physics buffer fields
!
   ifld = pbuf_get_index('FRACIS') 
   call pbuf_get_field(pbuf, ifld, fracis)

   fracice(:,:) = 0.0_r8


!
! Transport all constituents except cloud water and ice
!

   lchnk = state%lchnk
   ncol  = state%ncol
   nstep = get_nstep()

!
!     Convective transport of all trace species except cloud liquid 
!     and cloud ice done here because we need to do the scavenging first
!     to determine the interstitial fraction.
!
!  call cnst_get_ind('CLDLIQ', ixcldliq)
!  call cnst_get_ind('CLDICE', ixcldice)

! initialize dpdry (units=mb), which is used for tracers of dry mixing ratio type
   dpdry = 0._r8
   do i = 1, lengath
      dpdry(i,:) = state%pdeldry(ideep(i),:)/100._r8
   end do

! qaa hold tracer mixing ratios
   qaa = q

! turn on/off calculations for aerosols and trace gases
   do l = 1, pcnst
      dotend(l) = .false.
      if (species_class(l) == spec_class_aerosol) then
         if (convproc_do_aer) dotend(l) = .true.
      else if (species_class(l) == spec_class_gas) then
         if (convproc_do_gas) dotend(l) = .true.
      end if
   end do


! set itmpveca=1 for the 2 "strong convection" grid points, -1 otherwise
! set itmpveca=associated index to the gathered arrays
   call get_lat_all_p(  lchnk, pcols, lat_ndx )
   call get_lon_all_p(  lchnk, pcols, lon_ndx )
   call get_rlat_all_p( lchnk, pcols, lat_deg )
   call get_rlon_all_p( lchnk, pcols, lon_deg )
   lat_deg(1:ncol) = lat_deg(1:ncol) *180.0/3.1415926536
   lon_deg(1:ncol) = lon_deg(1:ncol) *180.0/3.1415926536

   itmpa = 0
   itmpveca(:) = -1
   itmpvecb(:) = -1

   if (idiag_aa > 0) then
   do i = 1, ncol
!     if (lat_ndx(i) /= 48) cycle
!     if ( (lon_ndx(i) /=  40) .and. &
!          (lon_ndx(i) /= 116) ) cycle
      if ( (lchnk == 146 .and. i == 10) .or. &
           (lchnk == 765 .and. i == 14) ) then
         continue
      else
         cycle
      end if

      itmpa = itmpa + 1
      itmpveca(i) = 1
      do ii = 1, lengath
         if (ideep(ii) == i) itmpvecb(i) = ii
      end do
      write(lun,'(a)')
      write(lun,'(a,i9,7i5,2f6.1)') &
         'qakn lchnk, ncol, i, lat_ndx, lon_ndx, itmpveca/b, nstep, lat, lon', &
         lchnk, ncol, i, lat_ndx(i), lon_ndx(i), itmpveca(i), itmpvecb(i), nstep, &
         lat_deg(i), lon_deg(i)
   end do ! i
   end if ! (idiag_aa > 0)


! change profiles of first 4 gases
   call cnst_get_ind('H2O2',   ixh2o2)
   call cnst_get_ind('so4_a1', ixso4_a1)
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
   call cnst_get_ind('bc_a5',  ixbc_a5)
#else
   call cnst_get_ind('bc_a1',  ixbc_a1)
#endif

   if (idiag_aa > 0 .and. idiag_gg > 0) then
   if (itmpa .gt. 0) then
      write(lun,'(a,2i4,1p,e12.4,i9,i5)') 'qako ixh2o2, ixso4_a1, dt, lchnk, nstep ', &
         ixh2o2, ixso4_a1, dt, lchnk, nstep
   end if
   if (ixh2o2   < 6 .or. ixh2o2   > pcnst-4) &
      call endrun( "*** ma_convproc_dp_intr -- bad ixh2o2" )
   if (ixso4_a1 < 6 .or. ixso4_a1 > pcnst) &
      call endrun( "*** ma_convproc_dp_intr -- bad ixso4_a1" )
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
   if (ixbc_a5  < 6 .or. ixbc_a5  > pcnst) &
      call endrun( "*** ma_convproc_dp_intr -- bad ixbc_a5" )
#else
   if (ixbc_a1  < 6 .or. ixbc_a1  > pcnst) &
      call endrun( "*** ma_convproc_dp_intr -- bad ixbc_a1" )
#endif

   do ll = 1, 4
      l = ixh2o2 + (ll-1)
      kbb = pver - 5*(ll-1) ; kbb = max( kbb, 1 )
      kaa = kbb -4          ; kaa = max( kaa, 1 )
      qaa(1:ncol,:,l) = 0.0_r8
!     qaa(1:ncol,kaa:kbb,l) = 1.0e-9_r8  ! qaa holds "modified q" before convtran
      qaa(1:ncol,kaa:kbb,l) = 1.0_r8  ! qaa holds "modified q" before convtran
   end do
   dotend(ixh2o2:ixh2o2+3) = .true.  ! for initial testing
   end if ! (idiag_aa > 0 .and. idiag_gg > 0)


! output profiles of some conv cloud variables
! (when idiag_aa <= 0, itmpa == 0, and this is inactive)
   if (itmpa > 0) then
   do i = 1, ncol
      if (itmpveca(i) <= 0) cycle

      ii = itmpvecb(i) ; kaa = 11 ; kbb = pver
      write(lun,'(a,i9,4i5,1p,2e12.4)') &
         'qakq lchnk, i, lat_ndx, lon_ndx, nstep                 ', &
         lchnk, i, lat_ndx(i), lon_ndx(i), nstep
      write(lun,'(2a)') &
         'qakq zi and fracis(so4_a1); mu and md; sum(eu+du), sum(ed); ', &
         'eu, du, ed; wup, dp_frac, icwmrdp (in-cld, grid-av), rprddp*dt'

      tmpveca(:) = 0.0 ; tmpvecb(:) = 0.0
      if (ii > 0) then
         k = pver
         tmpveca(k) = hund_ovr_g*(eu(ii,k)-du(ii,k))*dp(ii,k)
         do k = pver-1, 1, -1
            tmpveca(k) = tmpveca(k+1) + hund_ovr_g*(eu(ii,k)-du(ii,k))*dp(ii,k)
         end do
         do k = 2, pver
            tmpvecb(k) = tmpvecb(k-1) - hund_ovr_g*ed(ii,k-1)*dp(ii,k-1)
         end do
      end if

      tmpvecc(1:pver) = dp_frac(i,1:pver)
      do k = pver, 3, -1
         if (tmpvecc(k) < 1.0e-4_r8) then
            do kk = k-1, 1, -1
               if (tmpvecc(kk) >= 1.0e-4_r8) then
                  tmpvecc(k) = tmpvecc(kk) 
                  exit
               end if
            end do
         end if
         tmpvecc(k) = max( tmpvecc(k), 1.0e-4_r8 ) 
      end do

      tmpa = 0.0 ; tmpb = 0.0
      tmpc = 0.0 ; tmpd = 0.0
      tmpe = 0.0 ; tmpf = 0.0
      do k = 10, pver
         if (ii > 0) then
            tmpa = hund_ovr_g*mu(ii,k) ; tmpb = hund_ovr_g*md(ii,k)
            tmpc = hund_ovr_g*eu(ii,k) ; tmpd = hund_ovr_g*du(ii,k)
            tmpe = hund_ovr_g*ed(ii,k)
            tmpf = tmpa / ( (state%pmid(i,k)/(state%t(i,k)*rair)) * tmpvecc(k) )
         end if
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i4,1p, 4e11.3, 2x,2e11.3, 2x,3e11.3, &
              &2x,2e10.3, 2x,3e10.3)') 'qakq', k, &
            state%zi(i,k+1), fracis(i,k,ixso4_a1), tmpa, tmpb, &
            tmpveca(k), tmpvecb(k), tmpc, tmpd, tmpe, &
            tmpf, dp_frac(i,k), icwmrdp(i,k), &
            icwmrdp(i,k)*dp_frac(i,k), rprddp(i,k)*dt
      end do

      tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0
      if (ii > 0) then
         do k = pver-5, pver
            tmpd = hund_ovr_g*eu(ii,k)*dp(ii,k)
            tmpa = tmpa + tmpd*fracis(i,k,ixso4_a1)
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
            tmpb = tmpb + tmpd*fracis(i,k,ixbc_a5)
#else
            tmpb = tmpb + tmpd*fracis(i,k,ixbc_a1)
#endif
            tmpc = tmpc + tmpd
         end do
         tmpa = tmpa/ max(tmpc,1.0e-35_r8)
         tmpb = tmpb/ max(tmpc,1.0e-35_r8)
      end if
      write(lun,'(a,3i10,1p,3e10.2,2x,a)') 'qakq111222', &
         nstep, lchnk, i, tmpa, tmpb, tmpc, &
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
         'nstep, lchnk, i, k25:30-avg fracis(so4_a1), fracis(bc_a5), eu*dp'
#else
         'nstep, lchnk, i, k25:30-avg fracis(so4_a1), fracis(bc_a1), eu*dp'
#endif

      write(lun,'(/2a)') &
         'qakq --- numb and mass for modes 1-3'
      do k = 10, pver
         do n = 1, ntot_amode
            tmpveca(n) = qaa(i,k,numptr_amode(n))
            tmpvecb(n) = 0.0_r8
            do ll = 1, nspec_amode(n)
               tmpvecb(n) = tmpvecb(n) + qaa(i,k,lmassptr_amode(ll,n))
            end do
         end do
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i4,1p, 7(2x,2e10.2))' ) 'qakq ---', k, &
            (tmpveca(l), tmpvecb(l), l=1,min(ntot_amode,7))
      end do

   end do ! i
   end if ! (itmpa > 0)


!
! do ma_convproc_tend call
!
! question/issue - when computing first-order removal rate for convective cloud water,
!    should dlf be included as is done in wetdepa?
! detrainment does not change the in-cloud (= in updraft) cldwtr mixing ratio
! when you have detrainment, the updraft air mass flux is decreasing with height,
!    and the cldwtr flux may be decreasing also, 
!    but the in-cloud cldwtr mixing ratio is not changed by detrainment itself
! this suggests that wetdepa is incorrect, and dlf should not be included
!
! if dlf should be included, then you want to calculate
!    rprddp / (dp_frac*icwmrdp + dt*(rprddp + dlfdp)]
! so need to pass both rprddp and dlfdp to ma_convproc_tend
!

!  tine ma_convproc_tend(                                            &
!                    convtype,                                       &
!                    lchnk,      ncnst,      nstep,      dt,         &
!                    t,          pmid,       pdel,       q,          &   
!                    mu,         md,         du,         eu,         &   
!                    ed,         dp,         dsubcld,    jt,         &   
!                    mx,         ideep,      il1g,       il2g,       &       
!                    concld,     icwmr1,     cmfdqrzh,   fracice,    &
!                    dqdt,       doconvproc, nsrflx,     qsrflx      )

   call ma_convproc_tend(                                            &
                     'deep',                                         &
                     lchnk,      pcnst,      nstep,      dt,         &
                     state%t,    state%pmid, state%pdel, qaa,        &   
                     mu,         md,         du,         eu,         &   
                     ed,         dp,         dpdry,      jt,         &   
                     maxg,       ideep,      1,          lengath,    &       
                     dp_frac,    icwmrdp,    rprddp,     evapcdp,    &
                     fracice,                                        &
                     dqdt,       dotend,     nsrflx,     qsrflx,     &
                     species_class, mam_prevap_resusp_optaa,         & ! REASTER 08/05/2015
                     xx_mfup_max, xx_wcldbase, xx_kcldbase,          &
                     lun,        itmpveca                            )
!                    ed,         dp,         dsubcld,    jt,         &   



! set qbb = "modified q" with convtran tendency applied,
! and set tendencies to zero
   if (idiag_aa > 0 .and. idiag_gg > 0) then
   do ll = 1, 4
      l = ixh2o2 + (ll-1)
      qbb(1:ncol,:,l) = qaa(1:ncol,:,l) + dqdt(1:ncol,:,l)*dt
      dqdt(1:ncol,:,l) = 0.0  ! set tendency to zero
   end do
   end if ! (idiag_aa > 0 .and. idiag_gg > 0)

! output profiles of the "modified" gases
! (when idiag_aa <= 0, itmpa == 0, and this is inactive)
   if (itmpa > 0) then
   do i = 1, ncol
      if (itmpveca(i) <= 0) cycle

      ii = itmpvecb(i) ; kaa = 11 ; kbb = pver
      tmpa = 0.0
      if (ii > 0) then
         tmpa = 0.0
         tmpa = hund_ovr_g * sum(mu(ii,kaa:kbb)) / (kbb-kaa+1)
      end if
      write(lun,'(a,i9,4i5,1p,2e12.4)') &
         'qakp lchnk, i, lat_ndx, lon_ndx, nstep, mu-avg, dsubcld', &
         lchnk, i, lat_ndx(i), lon_ndx(i), nstep, tmpa, dsubcld(ii)
      write(lun,'(a,i9,4i5,1p,e12.4)') &
         'qakp zi and fracis(so4_a1); mu and md; gases 1-4'

      tmpa = 0.0 ; tmpb = 0.0
      do k = 10, pver
         if (ii > 0) then
            tmpa = 0.0
            tmpa = hund_ovr_g*mu(ii,k) ; tmpb = hund_ovr_g*md(ii,k)
         end if
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i4,1p,4e11.3,5(2x,2e10.3))') 'qakp', k, &
            state%zi(i,k+1), fracis(i,k,ixso4_a1), tmpa, tmpb, &
            (qaa(i,k,l), (qbb(i,k,l)-qaa(i,k,l)), l=ixh2o2,ixh2o2+3)
      end do

      tmpveca(:) = 0.0 ; tmpa = 0.0
      if (ii > 0) then
         do l = ixh2o2, ixh2o2+3
            tmpa = sum( dpdry(ii,:) )
            tmpveca(l) = sum( dpdry(ii,:) * (qbb(i,:,l)-qaa(i,:,l)) )
            tmpveca(l) = tmpveca(l) / max( tmpa, 1.0e-10_r8 )
         end do
      end if
      write(lun,'(a,i4,1p,e11.3,33x,5(12x,e10.3))') 'qakp', -1, &
         tmpa, (tmpveca(l), l=ixh2o2,ixh2o2+3)

      k = 26
      write(lun,'(a,i9,4i5,1p,2e12.4)') &
         'qakt lchnk, i, k, nstep // l, qold, qnew', &
         lchnk, i, k, nstep
      do l = ixh2o2, pcnst
         tmpa = qaa(i,k,l) + dt*dqdt(i,k,l) 
         if (tmpa .le. -1.0e-14_r8) write(lun,'(2a,i4,1p,2e12.4)') &
            'qakt ', cnst_name(l)(1:8), l, qaa(i,k,l), tmpa
      end do

   end do ! i
   end if ! (itmpa > 0)


   call outfld( 'DP_MFUP_MAX', xx_mfup_max, pcols, lchnk )
   call outfld( 'DP_WCLDBASE', xx_wcldbase, pcols, lchnk )
   call outfld( 'DP_KCLDBASE', xx_kcldbase, pcols, lchnk )


end subroutine ma_convproc_dp_intr



!=========================================================================================
subroutine ma_convproc_sh_intr(                 &
     state, pbuf, dt,                           &
     sh_frac, icwmrsh, rprdsh, evapcsh, dlfsh,  &
     cmfmcsh, sh_e_ed_ratio,                    &
     q, dqdt, dotend, nsrflx, qsrflx,           &
     species_class, mam_prevap_resusp_optaa     )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Convective cloud processing (transport, activation/resuspension,
!    wet removal) of aerosols and trace gases.
!    (Currently no aqueous chemistry and no trace-gas wet removal)
! Does aerosols    when convproc_do_aer is .true.
! Does trace gases when convproc_do_gas is .true.
!
! This routine does shallow convection
! Uses mass fluxes, cloud water, precip production from the
!    convective cloud routines
! 
! Author: R. Easter
! 
!-----------------------------------------------------------------------

   use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
   use time_manager,   only: get_nstep
   use physics_buffer, only: pbuf_get_index, physics_buffer_desc, pbuf_get_field
   use constituents,   only: pcnst, cnst_get_ind, cnst_name
   use error_messages, only: alloc_err	

   use physconst,      only: gravit, rair
   use phys_grid,      only: get_lat_all_p, get_lon_all_p, get_rlat_all_p, get_rlon_all_p

   use modal_aero_data, only: lmassptr_amode, nspec_amode, ntot_amode, numptr_amode
 
! Arguments
   type(physics_state), intent(in ) :: state          ! Physics state variables
   type(physics_buffer_desc), pointer :: pbuf(:)

   real(r8), intent(in) :: dt                         ! delta t (model time increment)

   real(r8), intent(in)    :: q(pcols,pver,pcnst)
   real(r8), intent(inout) :: dqdt(pcols,pver,pcnst)
   logical,  intent(inout) :: dotend(pcnst)
   integer,  intent(in)    :: nsrflx
   real(r8), intent(inout) :: qsrflx(pcols,pcnst,nsrflx)

   real(r8), intent(in)    :: sh_frac(pcols,pver) ! Shallow conv cloud frac (0-1)
   real(r8), intent(in)    :: icwmrsh(pcols,pver) ! Shallow conv cloud condensate (kg/kg - in cloud)
   real(r8), intent(in)    :: rprdsh(pcols,pver)  ! Shallow conv precip production (kg/kg/s - grid avg)
   real(r8), intent(in)    :: evapcsh(pcols,pver) ! Shallow conv precip evaporation (kg/kg/s - grid avg)
   real(r8), intent(in)    :: dlfsh(pcols,pver)   ! Shallow conv cldwtr detrainment (kg/kg/s - grid avg)
   real(r8), intent(in)    :: cmfmcsh(pcols,pverp) ! Shallow conv mass flux (kg/m2/s)
   real(r8), intent(in)    :: sh_e_ed_ratio(pcols,pver)  ! shallow conv [ent/(ent+det)] ratio
   integer,  intent(in)    :: species_class(:)
   integer,  intent(in)    :: mam_prevap_resusp_optaa

!  real(r8), intent(in)    :: concld(pcols,pver) ! Convective cloud cover

! Local variables
   integer, parameter :: idiag_aa = -1  ! controls diagnostic output at 2 selected grid points
   integer, parameter :: idiag_gg = -1  ! controls using special profiles for gases at these grid points

   integer :: i, ii, itmpa
   integer :: ido_20000_block_now
   integer :: ixcldice, ixcldliq              ! constituent indices for cloud liquid and ice water.
   integer :: ixh2o2, ixbc_a1, ixbc_a5, ixso4_a1
   integer :: itmpveca(pcols), itmpveca2(pcols)
   integer :: k, kaa, kbb, kcc, kk
   integer :: l, ll, lchnk, lun
   integer :: lat_ndx(pcols), lon_ndx(pcols)
   integer :: maxg_minval
   integer :: n, ncol, nstep

   real(r8) :: dpdry(pcols,pver)     ! layer delta-p-dry (mb)
   real(r8) :: fracice(pcols,pver)   ! Ice fraction of cloud droplets
   real(r8) :: lat_deg(pcols), lon_deg(pcols)
   real(r8) :: qaa(pcols,pver,pcnst), qbb(pcols,pver,pcnst)
   real(r8) :: tmpa, tmpb, tmpc, tmpd, tmpe, tmpf, tmpg
   real(r8) :: tmpveca(300), tmpvecb(300), tmpvecc(300)
   real(r8) :: xx_mfup_max(pcols), xx_wcldbase(pcols), xx_kcldbase(pcols)

! variables that mimic the zm-deep counterparts
                                               ! mu, md, ..., ideep, lengath are all deep conv variables
   real(r8)  :: mu(pcols,pver)   ! Updraft mass flux (positive)
   real(r8)  :: md(pcols,pver)   ! Downdraft mass flux (negative)
   real(r8)  :: du(pcols,pver)   ! Mass detrain rate from updraft
   real(r8)  :: eu(pcols,pver)   ! Mass entrain rate into updraft
   real(r8)  :: ed(pcols,pver)   ! Mass entrain rate into downdraft
                           ! eu, ed, du are "d(massflux)/dp" and are all positive
   real(r8)  :: dp(pcols,pver)   ! Delta pressure between interfaces

   integer   :: jt(pcols)         ! Index of cloud top for each column
   integer   :: maxg(pcols)       ! Index of cloud bot for each column
   integer   :: ideep(pcols)      ! Gathering array
   integer   :: lengath           ! Gathered min lon indices over which to operate

! physics buffer fields 
   integer itim, ifld
   real(r8), pointer, dimension(:,:,:) :: fracis  ! fraction of transported species that are insoluble

!
! Initialize
!
   lun = iulog

! call physics_ptend_init(ptend)

!
! Associate pointers with physics buffer fields
!
   ifld = pbuf_get_index('FRACIS')
   call pbuf_get_field(pbuf, ifld, fracis)
   
   fracice(:,:) = 0.0_r8


!
! Transport all constituents except cloud water and ice
!

   lchnk = state%lchnk
   ncol  = state%ncol
   nstep = get_nstep()

!
!     Convective transport of all trace species except cloud liquid 
!     and cloud ice done here because we need to do the scavenging first
!     to determine the interstitial fraction.
!
!  call cnst_get_ind('CLDLIQ', ixcldliq)
!  call cnst_get_ind('CLDICE', ixcldice)

!
! create mass flux, entrainment, detrainment, and delta-p arrays 
! with same units as the zm-deep
!
   mu(:,:) = 0.0_r8
   md(:,:) = 0.0_r8
   du(:,:) = 0.0_r8
   eu(:,:) = 0.0_r8
   ed(:,:) = 0.0_r8
   jt(:) = -1
   maxg(:) = -1
   ideep(:) = -1
   lengath = ncol
   maxg_minval = pver*2

! these dp and dpdry have units of mb
   dpdry(1:ncol,:) = state%pdeldry(1:ncol,:)/100._r8
   dp(   1:ncol,:) = state%pdel(   1:ncol,:)/100._r8

   do i = 1, ncol
      ideep(i) = i

! load updraft mass flux from cmfmcsh
      kk = 0
      do k = 2, pver
! if mass-flux < 1e-7 kg/m2/s ~= 1e-7 m/s ~= 1 cm/day, treat as zero
         if (cmfmcsh(i,k) >= 1.0e-7_r8) then
! mu has units of mb/s
            mu(i,k) = cmfmcsh(i,k) / hund_ovr_g
            kk = kk + 1
            if (kk == 1) jt(i) = k - 1
            maxg(i) = k
         end if
      end do ! k
      if (kk <= 0) cycle  ! current column has no convection
      
! extend below-cloud source region downwards (how far?)
      maxg_minval = min( maxg_minval, maxg(i) )
      kaa = maxg(i)
      kbb = min( kaa+4, pver )
!     kbb = pver
      if (kbb > kaa) then
         tmpa = sum( dpdry(i,kaa:kbb) )
         do k = kaa+1, kbb
            mu(i,k) = mu(i,kaa)*sum( dpdry(i,k:kbb) )/tmpa
         end do ! k
         maxg(i) = kbb
      end if

! calc ent / detrainment, using the [ent/(ent+det)] ratio from uw scheme
!    which is equal to [fer_out/(fer_out+fdr_out)]  (see uwshcu.F90)
!
! note that the ratio is set to -1.0 (invalid) when both fer and fdr are very small
!    and the ratio values are often strange (??) at topmost layer
!
! for initial testing, impose a limit of 
!    entrainment <= 4 * (net entrainment), OR
!    detrainment <= 4 * (net detrainment)
      do k = jt(i), maxg(i)
         if (k < pver) then
            tmpa = (mu(i,k) - mu(i,k+1))/dpdry(i,k)
         else
            tmpa = mu(i,k)/dpdry(i,k)
         end if
         tmpb = sh_e_ed_ratio(i,k)
!        tmpb = -1.0  ! force ent only or det only
         if (tmpb < -1.0e-5_r8) then
            ! do ent only or det only
            if (tmpa >= 0.0_r8) then
               ! net entrainment
               eu(i,k) = tmpa
            else
               ! net detrainment
               du(i,k) = -tmpa
            end if
         else
            if (tmpa >= 0.0_r8) then
               ! net entrainment
               if (k >= kaa .or. tmpb < 0.0_r8) then
                  ! layers at/below initial maxg, or sh_e_ed_ratio is invalid
                  eu(i,k) = tmpa
               else
                  tmpb = max( tmpb, 0.571_r8 )
                  eu(i,k) = tmpa*(tmpb/(2.0_r8*tmpb - 1.0_r8))
                  du(i,k) = eu(i,k) - tmpa
               end if
            else
               ! net detrainment
               tmpa = -tmpa
               if (k <= jt(i) .or. tmpb < 0.0_r8) then
                  ! layers at/above jt (where ratio is strange??), or sh_e_ed_ratio is invalid
                  du(i,k) = tmpa
               else
                  tmpb = min( tmpb, 0.429_r8 )
                  du(i,k) = tmpa*(1.0_r8 - tmpb)/(1.0_r8 - 2.0_r8*tmpb)
                  eu(i,k) = du(i,k) - tmpa
               end if
            end if
         end if
      end do ! k

   end do ! i


! qaa hold tracer mixing ratios
   qaa = q


! turn on/off calculations for aerosols and trace gases
   do l = 1, pcnst
      dotend(l) = .false.
      if (species_class(l) == spec_class_aerosol) then
         if (convproc_do_aer) dotend(l) = .true.
      else if (species_class(l) == spec_class_gas) then
         if (convproc_do_gas) dotend(l) = .true.
      end if
   end do


! set itmpveca=1 for the 2 "strong convection" grid points, -1 otherwise
! set itmpveca=associated index to the gathered arrays
   call get_lat_all_p(  lchnk, pcols, lat_ndx )
   call get_lon_all_p(  lchnk, pcols, lon_ndx )
   call get_rlat_all_p( lchnk, pcols, lat_deg )
   call get_rlon_all_p( lchnk, pcols, lon_deg )
   lat_deg(1:ncol) = lat_deg(1:ncol) *180.0/3.1415926536
   lon_deg(1:ncol) = lon_deg(1:ncol) *180.0/3.1415926536

   itmpa = 0
   itmpveca(:) = -1

   if (idiag_aa > 0) then
   do i = 1, ncol
!     if (lat_ndx(i) /= 50) cycle
!     if ( (lon_ndx(i) /=  30) .and. &
!          (lon_ndx(i) /= 139) ) cycle
      if ( (lchnk == 146 .and. i == 10) .or. &
           (lchnk == 765 .and. i == 14) ) then
         continue
      else
         cycle
      end if

      itmpa = itmpa + 1
      itmpveca(i) = 1
      write(lun,'(a)')
      write(lun,'(a,i9,6i5,2f6.1)') &
         'sqakn lchnk, ncol, i, lat_ndx, lon_ndx, itmpveca, nstep, lat, lon', &
         lchnk, ncol, i, lat_ndx(i), lon_ndx(i), itmpveca(i), nstep, &
         lat_deg(i), lon_deg(i)
   end do ! i
   end if ! (idiag_aa > 0)


! change profiles of first 4 gases
   call cnst_get_ind('H2O2',   ixh2o2)
   call cnst_get_ind('so4_a1', ixso4_a1)
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
   call cnst_get_ind('bc_a5',  ixbc_a5)
#else
   call cnst_get_ind('bc_a1',  ixbc_a1)
#endif

   if (idiag_aa > 0 .and. idiag_gg > 0) then
   if (itmpa .gt. 0) then
      write(lun,'(a,2i4,1p,e12.4,i9,i5)') 'sqako ixh2o2, ixso4_a1, dt, lchnk, nstep ', &
         ixh2o2, ixso4_a1, dt, lchnk, nstep
   end if
   if (ixh2o2   < 6 .or. ixh2o2   > pcnst-4) &
      call endrun( "*** ma_convproc_sh_intr -- bad ixh2o2" )
   if (ixso4_a1 < 6 .or. ixso4_a1 > pcnst) &
      call endrun( "*** ma_convproc_sh_intr -- bad ixso4_a1" )
#if ( defined MODAL_AERO_5MODE_AGEDCARBON )
   if (ixbc_a5  < 6 .or. ixbc_a5  > pcnst) &
      call endrun( "*** ma_convproc_sh_intr -- bad ixbc_a5" )
#else
   if (ixbc_a1  < 6 .or. ixbc_a1  > pcnst) &
      call endrun( "*** ma_convproc_sh_intr -- bad ixbc_a1" )
#endif

   do ll = 1, 4
      l = ixh2o2 + (ll-1)
      kbb = pver - 5*(ll-1) ; kbb = max( kbb, 1 )
      kaa = kbb -4          ; kaa = max( kaa, 1 )
      ! qaa holds "modified q" before convtran
      qaa(1:ncol,:,l) = 0.0_r8
      do k = kaa, kbb
        qaa(1:ncol,k,l) = 1.0e-9_r8*(0.9_r8 - 0.1_r8*(kbb-k))
      end do

      kbb = pver - 1 - 5*(ll-1)
      kaa = kbb - 4          ; kaa = max( kaa, 1 )
      kcc = kbb + 4          ; kcc = min( kcc, pver )
      qaa(1:ncol,:,l) = 0.0_r8
      do k = kaa, kcc
         qaa(1:ncol,k,l) = max( 0.0_r8, 1.0e-9_r8*(0.9_r8 - 0.2_r8*abs(kbb-k)) )
      end do
   end do
   dotend(ixh2o2:ixh2o2+3) = .true.  ! for initial testing
   end if ! (idiag_aa > 0 .and. idiag_gg > 0)


! output profiles of some conv cloud variables
! (when idiag_aa <= 0, itmpa == 0, and this is inactive)
   ido_20000_block_now = -77
   ido_20000_block_now = 1
20000 if (ido_20000_block_now > 0) then

   if (itmpa > 0) then
   do i = 1, ncol
      if (itmpveca(i) <= 0) cycle

      ii = i ; kaa = 11 ; kbb = pver

      tmpveca(:) = 0.0
      if (ii > 0) then
         k = pver
         tmpveca(k) = hund_ovr_g*(eu(ii,k)-du(ii,k))*dpdry(ii,k)
         do k = pver-1, 1, -1
            tmpveca(k) = tmpveca(k+1) + hund_ovr_g*(eu(ii,k)-du(ii,k))*dpdry(ii,k)
         end do
      end if

      tmpvecc(1:pver) = sh_frac(i,1:pver)
      do k = pver, 3, -1
         if (tmpvecc(k) < 1.0e-4_r8) then
            do kk = k-1, 1, -1
               if (tmpvecc(kk) >= 1.0e-4_r8) then
                  tmpvecc(k) = tmpvecc(kk) 
                  exit
               end if
            end do
         end if
         tmpvecc(k) = max( tmpvecc(k), 1.0e-4_r8 ) 
      end do

      write(lun,'(a,i9,4i5,1p,2e12.4)') &
         'sqakq lchnk, i, lat_ndx, lon_ndx, nstep                 ', &
         lchnk, i, lat_ndx(i), lon_ndx(i), nstep
      write(lun,'(2a)') &
         'sqakq zi; mu; sum(eu+du); ', &
         'eu, du; wup, sh_frac; icwmrsh (in-cld, grid-av), rprdsh*dt'

      tmpa = 0.0 ; tmpc = 0.0
      tmpd = 0.0 ; tmpf = 0.0
      do k = 10, pver
         if (ii > 0) then
            tmpa = hund_ovr_g*mu(ii,k)
            tmpc = hund_ovr_g*eu(ii,k) ; tmpd = hund_ovr_g*du(ii,k)
            tmpf = tmpa / ( (state%pmid(i,k)/(state%t(i,k)*rair)) * tmpvecc(k) )
         end if
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i4,1p, 2e11.3, 2x,e11.3, 2x,2e11.3, &
              &2x,2e10.3, 2x,3e10.3)') 'sqakq', k, &
            state%zi(i,k+1), tmpa, &
            tmpveca(k),    tmpc, tmpd, &
            tmpf, sh_frac(i,k),    icwmrsh(i,k), &
            icwmrsh(i,k)*sh_frac(i,k), rprdsh(i,k)*dt
      end do

      write(lun,'(/2a)') &
         'sqakq --- numb and mass for modes 1-3'
      do k = 10, pver
!        exit
         do n = 1, ntot_amode
            tmpveca(n) = qaa(i,k,numptr_amode(n))
            tmpvecb(n) = 0.0_r8
            do ll = 1, nspec_amode(n)
               tmpvecb(n) = tmpvecb(n) + qaa(i,k,lmassptr_amode(ll,n))
            end do
         end do
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i4,1p, 7(2x,2e10.2))' ) 'sqakq ---', k, &
            (tmpveca(l), tmpvecb(l), l=1,min(ntot_amode,7))
      end do

   end do ! i
   end if ! (itmpa > 0)

   end if ! (ido_20000_block_now > 0) then
   if (ido_20000_block_now == 77) goto 40000


!
! do ma_convproc_tend call
!
! question/issue - when computing first-order removal rate for convective cloud water,
!    should dlf be included as is done in wetdepa?
! detrainment does not change the in-cloud (= in updraft) cldwtr mixing ratio
! when you have detrainment, the updraft air mass flux is decreasing with height,
!    and the cldwtr flux may be decreasing also, 
!    but the cldwtr mixing ratio does not change
! this suggests that wetdepa is incorrect, and dlf should not be included
!
! if dlf should be included, then you want to calculate
!    rprddp / (dp_frac*icwmrdp + dt*(rprddp + dlfdp)]
! so need to pass both rprddp and dlfdp to ma_convproc_tend
!

!  tine ma_convproc_tend(                                            &
!                    convtype,                                       &
!                    lchnk,      ncnst,      nstep,      dt,         &
!                    t,          pmid,       pdel,       q,          &   
!                    mu,         md,         du,         eu,         &   
!                    ed,         dp,         dsubcld,    jt,         &   
!                    mx,         ideep,      il1g,       il2g,       &       
!                    concld,     icwmr1,     cmfdqrzh,   fracice,    &
!                    dqdt,       doconvproc, nsrflx,     qsrflx      )

   itmpveca2 = itmpveca
!  itmpveca2 = -1
   call ma_convproc_tend(                                            &
                     'uwsh',                                         &
                     lchnk,      pcnst,      nstep,      dt,         &
                     state%t,    state%pmid, state%pdel, qaa,        &   
                     mu,         md,         du,         eu,         &   
                     ed,         dp,         dpdry,      jt,         &   
                     maxg,       ideep,      1,          lengath,    &       
                     sh_frac,    icwmrsh,    rprdsh,     evapcsh,    &
                     fracice,                                        &
                     dqdt,       dotend,     nsrflx,     qsrflx,     &
                     species_class, mam_prevap_resusp_optaa,         & ! REASTER 08/05/2015
                     xx_mfup_max, xx_wcldbase, xx_kcldbase,          &
                     lun,        itmpveca2                           )



! set qbb = "modified q" with convtran tendency applied,
! and set tendencies to zero
   if (idiag_aa > 0 .and. idiag_gg > 0) then
   do ll = 1, 4
      l = ixh2o2 + (ll-1)
      qbb(1:ncol,:,l) = qaa(1:ncol,:,l) + dqdt(1:ncol,:,l)*dt
      dqdt(1:ncol,:,l) = 0.0  ! set tendency to zero
   end do
   end if ! (idiag_aa > 0 .and. idiag_gg > 0)

! output profiles of the "modified" gases
! (when idiag_aa <= 0, itmpa == 0, and this is inactive)
   if (itmpa > 0) then
   do i = 1, ncol
      if (itmpveca(i) <= 0) cycle

      ii = i ; kaa = 11 ; kbb = pver
      tmpa = 0.0
      if (ii > 0) then
         tmpa = 0.0
         tmpa = hund_ovr_g * sum(mu(ii,kaa:kbb)) / (kbb-kaa+1)
      end if
      write(lun,'(a,i9,4i5,1p,e12.4)') &
         'sqakp lchnk, i, lat_ndx, lon_ndx, nstep, mu-avg ', &
         lchnk, i, lat_ndx(i), lon_ndx(i), nstep, tmpa
      write(lun,'(a,i9,4i5,1p,e12.4)') &
         'sqakp zi and mu; eu and du ; gases 1-4'

      tmpa = 0.0 ; tmpb = 0.0 ; tmpc = 0.0 
      do k = 10, pver
         if (ii > 0) then
            tmpa = hund_ovr_g*mu(ii,k)
            tmpb = hund_ovr_g*eu(ii,k)
            tmpc = hund_ovr_g*du(ii,k)
         end if
         if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
         write(lun,'(a,i3,1p,2e11.3,(2x,2e10.3),4(e11.1,e10.2))') 'sqakp', k, &
            state%zi(i,k+1), tmpa, tmpb, tmpc, &
            (qaa(i,k,l)*1.0e9, (qbb(i,k,l)-qaa(i,k,l))*1.0e9, l=ixh2o2,ixh2o2+3)
      end do

      tmpveca(:) = 0.0 ; tmpa = 0.0
      if (ii > 0) then
         do l = ixh2o2, ixh2o2+3
            tmpa = sum( dpdry(ii,:) )
            tmpveca(l) = sum( dpdry(ii,:) * (qbb(i,:,l)-qaa(i,:,l)) )
            tmpveca(l) = tmpveca(l) / max( tmpa, 1.0e-10_r8 )
         end do
      end if
      write(lun,'(a,i3,1p,e11.3,11x,22x,4(11x,e10.2))') 'sqakp', -1, &
         tmpa, (tmpveca(l), l=ixh2o2,ixh2o2+3)

      k = 26
      write(lun,'(a,i9,4i5,1p,2e12.4)') &
         'sqakt lchnk, i, k, nstep // l, qold, qnew', &
         lchnk, i, k, nstep
      do l = ixh2o2, pcnst
         tmpa = qaa(i,k,l) + dt*dqdt(i,k,l) 
         if (tmpa .le. -1.0e-14_r8) write(lun,'(2a,i4,1p,2e12.4)') &
            'sqakt ', cnst_name(l)(1:8), l, qaa(i,k,l), tmpa
      end do

   end do ! i
   end if ! (itmpa > 0)


   if (ido_20000_block_now == -77) then
      ido_20000_block_now = 77
      goto 20000
   end if
40000 continue


!  if (maxg_minval <= pver) write(lun,'(i3.3,i9,a)') &
!     maxg_minval, lchnk, '  sqak999888 lchnk, maxg_minval'


   call outfld( 'SH_MFUP_MAX', xx_mfup_max, pcols, lchnk )
   call outfld( 'SH_WCLDBASE', xx_wcldbase, pcols, lchnk )
   call outfld( 'SH_KCLDBASE', xx_kcldbase, pcols, lchnk )


end subroutine ma_convproc_sh_intr



!=========================================================================================
subroutine ma_convproc_tend(                                           &
                     convtype,                                       &
                     lchnk,      ncnst,      nstep,      dt,         &
                     t,          pmid,       pdel,       q,          &   
                     mu,         md,         du,         eu,         &   
                     ed,         dp,         dpdry,      jt,         &   
                     mx,         ideep,      il1g,       il2g,       &       
                     cldfrac,    icwmr,      rprd,       evapc,      &
                     fracice,                                        &
                     dqdt,       doconvproc, nsrflx,     qsrflx,     &
                     species_class, mam_prevap_resusp_optaa,         & ! REASTER 08/05/2015
                     xx_mfup_max, xx_wcldbase, xx_kcldbase,          &
                     lun,        idiag_in                            )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Convective transport of trace species.
! The trace species need not be conservative, and source/sink terms for
!    activation, resuspension, aqueous chemistry and gas uptake, and
!    wet removal are all applied.
! Currently this works with the ZM deep convection, but we should be able
!    to adapt it for both Hack and McCaa shallow convection
!
!
! Compare to subr convproc which does conservative trace species.
!
! A distinction between "moist" and "dry" mixing ratios is not currently made.
! (P. Rasch comment:  Note that we are still assuming that the tracers are 
!  in a moist mixing ratio this will change soon)

! 
! Method: 
! Computes tracer mixing ratios in updraft and downdraft "cells" in a
! Lagrangian manner, with source/sinks applied in the updraft other.
! Then computes grid-cell-mean tendencies by considering
!    updraft and downdraft fluxes across layer boundaries
!    environment subsidence/lifting fluxes across layer boundaries
!    sources and sinks in the updraft
!    resuspension of activated species in the grid-cell as a whole
!
! Note1:  A better estimate or calculation of either the updraft velocity
!         or fractional area is needed.
! Note2:  If updraft area is a small fraction of over cloud area, 
!         then aqueous chemistry is underestimated.  These are both
!         research areas.
! 
! Authors: O. Seland and R. Easter, based on convtran by P. Rasch
! 
!-----------------------------------------------------------------------

   use shr_kind_mod, only: r8=>shr_kind_r8
   use ppgrid, only: pcols, pver
   use physconst, only: gravit, rair, rhoh2o
   use constituents, only: pcnst, cnst_name

   use modal_aero_data, only:  cnst_name_cw, &
      lmassptr_amode, lmassptrcw_amode, &
      ntot_amode, ntot_amode, &
      nspec_amode, numptr_amode, numptrcw_amode
!  use units, only: getunit

   implicit none

!-----------------------------------------------------------------------
!
! Input arguments
!
   character(len=*), intent(in) :: convtype  ! identifies the type of
                                             ! convection ("deep", "shcu")
   integer,  intent(in) :: lchnk             ! chunk identifier
   integer,  intent(in) :: ncnst             ! number of tracers to transport
   integer,  intent(in) :: nstep             ! Time step index
   real(r8), intent(in) :: dt                ! Model timestep
   real(r8), intent(in) :: t(pcols,pver)     ! Temperature
   real(r8), intent(in) :: pmid(pcols,pver)  ! Pressure at model levels
   real(r8), intent(in) :: pdel(pcols,pver)  ! Pressure thickness of levels
   real(r8), intent(in) :: q(pcols,pver,ncnst) ! Tracer array including moisture

   real(r8), intent(in) :: mu(pcols,pver)    ! Updraft mass flux (positive)
   real(r8), intent(in) :: md(pcols,pver)    ! Downdraft mass flux (negative)
   real(r8), intent(in) :: du(pcols,pver)    ! Mass detrain rate from updraft
   real(r8), intent(in) :: eu(pcols,pver)    ! Mass entrain rate into updraft
   real(r8), intent(in) :: ed(pcols,pver)    ! Mass entrain rate into downdraft
! *** note1 - mu, md, eu, ed, du, dp, dpdry are GATHERED ARRAYS ***
! *** note2 - mu and md units are (mb/s), which is used in the zm_conv code
!           - eventually these should be changed to (kg/m2/s)
! *** note3 - eu, ed, du are "d(massflux)/dp" (with dp units = mb), and are all >= 0

   real(r8), intent(in) :: dp(pcols,pver)    ! Delta pressure between interfaces (mb)
   real(r8), intent(in) :: dpdry(pcols,pver) ! Delta dry-pressure (mb)
!  real(r8), intent(in) :: dsubcld(pcols)    ! Delta pressure from cloud base to sfc
   integer,  intent(in) :: jt(pcols)         ! Index of cloud top for each column
   integer,  intent(in) :: mx(pcols)         ! Index of cloud top for each column
   integer,  intent(in) :: ideep(pcols)      ! Gathering array indices 
   integer,  intent(in) :: il1g              ! Gathered min lon indices over which to operate
   integer,  intent(in) :: il2g              ! Gathered max lon indices over which to operate
! *** note4 -- for il1g <= i <= il2g,  icol = ideep(i) is the "normal" chunk column index

   real(r8), intent(in) :: cldfrac(pcols,pver)  ! Convective cloud fractional area
   real(r8), intent(in) :: icwmr(pcols,pver)    ! Convective cloud water from zhang
   real(r8), intent(in) :: rprd(pcols,pver)     ! Convective precipitation formation rate
   real(r8), intent(in) :: evapc(pcols,pver)    ! Convective precipitation evaporation rate
   real(r8), intent(in) :: fracice(pcols,pver)  ! Ice fraction of cloud droplets

   real(r8), intent(out):: dqdt(pcols,pver,ncnst)  ! Tracer tendency array
   logical,  intent(in) :: doconvproc(ncnst) ! flag for doing convective transport
   integer,  intent(in) :: nsrflx            ! last dimension of qsrflx
   real(r8), intent(out):: qsrflx(pcols,pcnst,nsrflx)
                              ! process-specific column tracer tendencies
                              !  1 = activation   of interstial to conv-cloudborne
                              !  2 = resuspension of conv-cloudborne to interstital
                              !  3 = aqueous chemistry (not implemented yet, so zero)
                              !  4 = wet removal
                              !  5 = actual precip-evap resuspension (what actually is applied to a species)
                              !  6 = pseudo precip-evap resuspension (for history file) ! REASTER 08/05/2015
   integer,  intent(in) :: species_class(:) ! REASTER 08/05/2015
   integer,  intent(in) :: mam_prevap_resusp_optaa ! REASTER 08/05/2015
   real(r8), intent(out):: xx_mfup_max(pcols)
   real(r8), intent(out):: xx_wcldbase(pcols)
   real(r8), intent(out):: xx_kcldbase(pcols)
   integer,  intent(in) :: lun               ! unit number for diagnostic output
   integer,  intent(in) :: idiag_in(pcols)   ! flag for diagnostic output


!--------------------------Local Variables------------------------------

! cloudborne aerosol, so the arrays are dimensioned with pcnst_extd = pcnst*2
   integer, parameter :: pcnst_extd = pcnst*2

   integer :: i, icol         ! Work index
   integer :: iconvtype       ! 1=deep, 2=uw shallow
   integer :: idiag_act       ! Work index
   integer :: iflux_method    ! 1=as in convtran (deep), 2=simpler
   integer :: ipass_calc_updraft
   integer :: itmpa, itmpb    ! Work variable
   integer :: j, jtsub        ! Work index
   integer :: k               ! Work index
   integer :: kactcnt         ! Counter for no. of levels having activation
   integer :: kactcntb        ! Counter for activation diagnostic output
   integer :: kactfirst       ! Lowest layer with activation (= cloudbase)
   integer :: kbot            ! Cloud-flux bottom layer for current i (=mx(i))
   integer :: kbot_prevap     ! Lowest layer for doing resuspension from evaporating precip ! REASTER 08/05/2015
   integer :: ktop            ! Cloud-flux top    layer for current i (=jt(i))
                              ! Layers between kbot,ktop have mass fluxes
                              !    but not all have cloud water, because the
                              !    updraft starts below the cloud base
   integer :: km1, km1x       ! Work index
   integer :: kp1, kp1x       ! Work index
   integer :: l, ll, la, lc   ! Work index
   integer :: m, n            ! Work index
   integer :: merr            ! number of errors (i.e., failed diagnostics)
                              ! for current column
   integer :: nerr            ! number of errors for entire run
   integer :: nerrmax         ! maximum number of errors to report
   integer :: ncnst_extd
   integer :: npass_calc_updraft
   integer :: ntsub           ! 

   logical  do_act_this_lev             ! flag for doing activation at current level
   logical  doconvproc_extd(pcnst_extd) ! flag for doing convective transport

   real(r8) aqfrac(pcnst_extd)       ! aqueous fraction of constituent in updraft
   real(r8) cldfrac_i(pver)          ! cldfrac at current i (with adjustments)

   real(r8) chat(pcnst_extd,pverp)   ! mix ratio in env at interfaces
   real(r8) cond(pcnst_extd,pverp)   ! mix ratio in downdraft at interfaces
   real(r8) const(pcnst_extd,pver)   ! gathered tracer array
   real(r8) conu(pcnst_extd,pverp)   ! mix ratio in updraft at interfaces

   real(r8) dcondt(pcnst_extd,pver)  ! grid-average TMR tendency for current column
   real(r8) dcondt_prevap(pcnst_extd,pver) ! portion of dcondt from precip evaporation
   real(r8) dcondt_prevap_hist(pcnst_extd,pver) ! similar but used for history output ! REASTER 08/05/2015
   real(r8) dcondt_resusp(pcnst_extd,pver) ! portion of dcondt from resuspension

   real(r8) dcondt_wetdep(pcnst_extd,pver) ! portion of dcondt from wet deposition
   real(r8) dconudt_activa(pcnst_extd,pverp) ! d(conu)/dt by activation
   real(r8) dconudt_aqchem(pcnst_extd,pverp) ! d(conu)/dt by aqueous chem
   real(r8) dconudt_wetdep(pcnst_extd,pverp) ! d(conu)/dt by wet removal

   real(r8) maxflux(pcnst_extd)      ! maximum (over layers) of fluxin and fluxout
   real(r8) maxflux2(pcnst_extd)     ! ditto but computed using method-2 fluxes
   real(r8) maxprevap(pcnst_extd)    ! maximum (over layers) of dcondt_prevap*dp
   real(r8) maxresusp(pcnst_extd)    ! maximum (over layers) of dcondt_resusp*dp
   real(r8) maxsrce(pcnst_extd)      ! maximum (over layers) of netsrce

   real(r8) sumflux(pcnst_extd)      ! sum (over layers) of netflux
   real(r8) sumflux2(pcnst_extd)     ! ditto but computed using method-2 fluxes
   real(r8) sumsrce(pcnst_extd)      ! sum (over layers) of dp*netsrce
   real(r8) sumchng(pcnst_extd)      ! sum (over layers) of dp*dcondt
   real(r8) sumchng3(pcnst_extd)     ! ditto but after call to resusp_conv
   real(r8) sumactiva(pcnst_extd)    ! sum (over layers) of dp*dconudt_activa
   real(r8) sumaqchem(pcnst_extd)    ! sum (over layers) of dp*dconudt_aqchem
   real(r8) sumprevap(pcnst_extd)    ! sum (over layers) of dp*dcondt_prevap
   real(r8) sumprevap_hist(pcnst_extd) ! sum (over layers) of dp*dcondt_prevap_hist ! REASTER 08/05/2015
   real(r8) sumresusp(pcnst_extd)    ! sum (over layers) of dp*dcondt_resusp
   real(r8) sumwetdep(pcnst_extd)    ! sum (over layers) of dp*dconudt_wetdep

   real(r8) cabv                 ! mix ratio of constituent above
   real(r8) cbel                 ! mix ratio of constituent below
   real(r8) cdifr                ! normalized diff between cabv and cbel
   real(r8) cdt(pver)            ! (in-updraft first order wet removal rate) * dt
   real(r8) clw_cut              ! threshold clw value for doing updraft
                                 ! transformation and removal
   real(r8) courantmax           ! maximum courant no.
   real(r8) dddp(pver)           ! dd(i,k)*dp(i,k) at current i
   real(r8) dp_i(pver)           ! dp(i,k) at current i
   real(r8) dt_u(pver)           ! lagrangian transport time in the updraft  
   real(r8) dudp(pver)           ! du(i,k)*dp(i,k) at current i
   real(r8) dqdt_i(pver,pcnst)   ! dqdt(i,k,m) at current i
   real(r8) dtsub                ! dt/ntsub
   real(r8) dz                   ! working layer thickness (m) 
   real(r8) eddp(pver)           ! ed(i,k)*dp(i,k) at current i
   real(r8) eudp(pver)           ! eu(i,k)*dp(i,k) at current i
   real(r8) expcdtm1             ! a work variable
   real(r8) fa_u(pver)           ! fractional area of in the updraft  
   real(r8) fa_u_dp              ! current fa_u(k)*dp_i(k)
   real(r8) f_ent                ! fraction of the "before-detrainment" updraft
                                 ! massflux at k/k-1 interface resulting from
                                 ! entrainment of level k air
   real(r8) fluxin               ! a work variable
   real(r8) fluxout              ! a work variable
   real(r8) maxc                 ! a work variable
   real(r8) mbsth                ! Threshold for mass fluxes
   real(r8) minc                 ! a work variable
   real(r8) md_m_eddp            ! a work variable
   real(r8) md_i(pverp)          ! md(i,k) at current i (note pverp dimension)
   real(r8) md_x(pverp)          ! md(i,k) at current i (note pverp dimension)
   real(r8) mu_i(pverp)          ! mu(i,k) at current i (note pverp dimension)
   real(r8) mu_x(pverp)          ! mu(i,k) at current i (note pverp dimension)
   ! md_i, md_x, mu_i, mu_x are all "dry" mass fluxes
   ! the mu_x/md_x are initially calculated from the incoming mu/md by applying dp/dpdry
   ! the mu_i/md_i are next calculated by applying the mbsth threshold
   real(r8) mu_p_eudp(pver)      ! = mu_i(kp1) + eudp(k)
   real(r8) netflux              ! a work variable
   real(r8) netsrce              ! a work variable
   real(r8) q_i(pver,pcnst)      ! q(i,k,m) at current i
   real(r8) qsrflx_i(pcnst,nsrflx) ! qsrflx(i,m,n) at current i
   real(r8) relerr_cut           ! relative error criterion for diagnostics
   real(r8) rhoair_i(pver)       ! air density at current i
   real(r8) small                ! a small number
   real(r8) tmpa, tmpb, tmpc     ! work variables
   real(r8) tmpf                 ! work variables
   real(r8) tmpveca(pcnst_extd)  ! work variables
   real(r8) tmpmata(pcnst_extd,3) ! work variables
   real(r8) xinv_ntsub           ! 1.0/ntsub
   real(r8) wup(pver)            ! working updraft velocity (m/s)
   real(r8) zmagl(pver)          ! working height above surface (m)
   real(r8) zkm                  ! working height above surface (km)

   character(len=16) :: cnst_name_extd(pcnst_extd)

   !PMA
   real(r8), parameter :: sfadj = 3.5_r8 !2.5_r8
!-----------------------------------------------------------------------
!

!  if (nstep > 1) call endrun()

   if (convtype == 'deep') then
      iconvtype = 1
      iflux_method = 1
   else if (convtype == 'uwsh') then
      iconvtype = 2
      iflux_method = 2
   else
      call endrun( '*** ma_convproc_tend -- convtype is not |deep| or |uwsh|' )
   end if

   nerr = 0
   nerrmax = 99

   ncnst_extd = pcnst_extd


   small = 1.e-36
! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
   mbsth = 1.e-15

   qsrflx(:,:,:) = 0.0_r8
   dqdt(:,:,:) = 0.0_r8
   xx_mfup_max(:) = 0.0_r8
   xx_wcldbase(:) = 0.0_r8
   xx_kcldbase(:) = 0.0_r8

! set doconvproc_extd (extended array) values
! inititialize aqfrac to 1.0 for activated aerosol species, 0.0 otherwise
   doconvproc_extd(:) = .false.
   doconvproc_extd(2:ncnst) = doconvproc(2:ncnst)
   aqfrac(:) = 0.0_r8
   do n = 1, ntot_amode
      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
         end if
         if ( doconvproc(la) ) then
            doconvproc_extd(lc) = .true.
            aqfrac(lc) = 1.0_r8
         end if
      enddo
   enddo ! n

   do l = 1, pcnst_extd
      if (l <= pcnst) then
         cnst_name_extd(l) = cnst_name(l)
      else
         cnst_name_extd(l) = trim(cnst_name(l-pcnst)) // '_cw'
      end if
   end do


! Loop ever each column that has convection
! *** i is index to gathered arrays; ideep(i) is index to "normal" chunk arrays
i_loop_main_aa: &
   do i = il1g, il2g
   icol = ideep(i)


   if ( (jt(i) <= 0) .and. (mx(i) <= 0) .and. (iconvtype /= 1) ) then
! shallow conv case with jt,mx <= 0, which means there is no shallow conv
! in this column -- skip this column
      cycle i_loop_main_aa

   else if ( (jt(i) < 1) .or. (mx(i) > pver) .or. (jt(i) > mx(i)) ) then
! invalid cloudtop and cloudbase indices -- skip this column
      write(lun,9010) 'illegal jt, mx', convtype, lchnk, icol, i,    &
                                      jt(i), mx(i)
9010  format( '*** ma_convproc_tend error -- ', a, 5x, 'convtype = ', a /   &
              '*** lchnk, icol, il, jt, mx = ', 5(1x,i10) )
      cycle i_loop_main_aa

   else if (jt(i) == mx(i)) then
! cloudtop = cloudbase (1 layer cloud) -- skip this column
      write(lun,9010) 'jt == mx', convtype, lchnk, icol, i, jt(i), mx(i)
      cycle i_loop_main_aa

   end if


!
! cloudtop and cloudbase indices are valid so proceed with calculations
!

! Load dp_i and cldfrac_i, and calc rhoair_i
      do k = 1, pver
         dp_i(k) = dpdry(i,k)
         cldfrac_i(k) = cldfrac(icol,k)
         rhoair_i(k) = pmid(icol,k)/(rair*t(icol,k))
      end do

! Calc dry mass fluxes
!    This is approximate because the updraft air is has different temp and qv than
!    the grid mean, but the whole convective parameterization is highly approximate
      mu_x(:) = 0.0
      md_x(:) = 0.0
! (eu-du) = d(mu)/dp -- integrate upwards, multiplying by dpdry
      do k = pver, 1, -1
         mu_x(k) = mu_x(k+1) + (eu(i,k)-du(i,k))*dp_i(k)
         xx_mfup_max(icol) = max( xx_mfup_max(icol), mu_x(k) )
      end do
! (ed) = d(md)/dp -- integrate downwards, multiplying by dpdry
      do k = 2, pver
         md_x(k) = md_x(k-1) - ed(i,k-1)*dp_i(k-1)
      end do

! Load mass fluxes over cloud layers
! (Note - use of arrays dimensioned k=1,pver+1 simplifies later coding)
! Zero out values below threshold
! Zero out values at "top of cloudtop", "base of cloudbase"
      ktop = jt(i)
      kbot = mx(i)
! REASTER 08/05/2015 BEGIN
! usually the updraft ( & downdraft) start ( & end ) at kbot=pver, but sometimes kbot < pver
! transport, activation, resuspension, and wet removal only occur between kbot >= k >= ktop
! resuspension from evaporating precip can occur at k > kbot when kbot < pver
! in the first version of this routine, the precp evap resusp tendencies for k > kbot were ignored, 
!    but that is now fixed
! this was a minor bug with quite minor affects on the aerosol, 
!    because convective precip evap is (or used to be) much less than stratiform precip evap )
      kbot_prevap = kbot
      if ( convproc_prevap_resusp_fixaa ) kbot_prevap = pver
! REASTER 08/05/2015 END
      mu_i(:) = 0.0
      md_i(:) = 0.0
      do k = ktop+1, kbot
         mu_i(k) = mu_x(k)
         if (mu_i(k) <= mbsth) mu_i(k) = 0.0
         md_i(k) = md_x(k)
         if (md_i(k) >= -mbsth) md_i(k) = 0.0
      end do
      mu_i(ktop) = 0.0
      md_i(ktop) = 0.0
      mu_i(kbot+1) = 0.0
      md_i(kbot+1) = 0.0

!  Compute updraft and downdraft "entrainment*dp" from eu and ed
!  Compute "detrainment*dp" from mass conservation
      eudp(:) = 0.0
      dudp(:) = 0.0
      eddp(:) = 0.0
      dddp(:) = 0.0
      courantmax = 0.0
      do k = ktop, kbot
         if ((mu_i(k) > 0) .or. (mu_i(k+1) > 0)) then
            if (du(i,k) <= 0.0) then
               eudp(k) = mu_i(k) - mu_i(k+1) 
            else
               eudp(k) = max( eu(i,k)*dp_i(k), 0.0_r8 )
               dudp(k) = (mu_i(k+1) + eudp(k)) - mu_i(k) 
               if (dudp(k) < 1.0e-12*eudp(k)) then
                  eudp(k) = mu_i(k) - mu_i(k+1) 
                  dudp(k) = 0.0
               end if
            end if
         end if
         if ((md_i(k) < 0) .or. (md_i(k+1) < 0)) then
            eddp(k) = max( ed(i,k)*dp_i(k), 0.0_r8 )
            dddp(k) = (md_i(k+1) + eddp(k)) - md_i(k) 
            if (dddp(k) < 1.0e-12*eddp(k)) then
               eddp(k) = md_i(k) - md_i(k+1) 
               dddp(k) = 0.0
            end if
         end if
!        courantmax = max( courantmax, (eudp(k)+eddp(k))*dt/dp_i(k) )  ! old version - incorrect
         courantmax = max( courantmax, ( mu_i(k+1)+eudp(k)-md_i(k)+eddp(k) )*dt/dp_i(k) )
      end do ! k

! number of time substeps needed to maintain "courant number" <= 1
      ntsub = 1
      if (courantmax > (1.0_r8 + 1.0e-6_r8)) then
         ntsub = 1 + int( courantmax )
      end if
      xinv_ntsub = 1.0_r8/ntsub
      dtsub = dt*xinv_ntsub
      courantmax = courantmax*xinv_ntsub

! zmagl(k) = height above surface for middle of level k
      zmagl(pver) = 0.0
      do k = pver, 1, -1
         if (k < pver) then
            zmagl(k) = zmagl(k+1) + 0.5*dz
         end if
         dz = dp_i(k)*hund_ovr_g/rhoair_i(k)
         zmagl(k) = zmagl(k) + 0.5*dz
      end do

!  load tracer mixing ratio array, which will be updated at the end of each jtsub interation
      q_i(1:pver,1:pcnst) = q(icol,1:pver,1:pcnst)

!
!   when method_reduce_actfrac = 2, need to do the updraft calc twice
!   (1st to get non-adjusted activation amount, 2nd to apply reduction factor)
      npass_calc_updraft = 1
      if ( (method_reduce_actfrac == 2)      .and. &
           (factor_reduce_actfrac >= 0.0_r8) .and. &
           (factor_reduce_actfrac <= 1.0_r8) ) npass_calc_updraft = 2


jtsub_loop_main_aa: &
      do jtsub = 1, ntsub


ipass_calc_updraft_loop: &
      do ipass_calc_updraft = 1, npass_calc_updraft


      if (idiag_in(icol) > 0) &
         write(lun,'(/a,3x,a,1x,i9,5i5)') 'qakr - convtype,lchnk,i,jt,mx,jtsub,ipass=', &
            trim(convtype), lchnk, icol, jt(i), mx(i), jtsub, ipass_calc_updraft

      qsrflx_i(:,:) = 0.0
      dqdt_i(:,:) = 0.0

      const(:,:) = 0.0 ! zero cloud-phase species
      chat(:,:) = 0.0 ! zero cloud-phase species
      conu(:,:) = 0.0
      cond(:,:) = 0.0

      dcondt(:,:) = 0.0
      dcondt_resusp(:,:) = 0.0
      dcondt_wetdep(:,:) = 0.0
      dcondt_prevap(:,:) = 0.0
      dcondt_prevap_hist(:,:) = 0.0 ! REASTER 08/05/2015
      dconudt_aqchem(:,:) = 0.0
      dconudt_wetdep(:,:) = 0.0
! only initialize the activation tendency on ipass=1
      if (ipass_calc_updraft == 1) dconudt_activa(:,:) = 0.0

! initialize mixing ratio arrays (chat, const, conu, cond)
      do m = 2, ncnst
      if ( doconvproc_extd(m) ) then

! Gather up the constituent
         do k = 1,pver
            const(m,k) = q_i(k,m)
         end do
         
! From now on work only with gathered data
! Interpolate environment tracer values to interfaces
         do k = 1,pver
            km1 = max(1,k-1)
            minc = min(const(m,km1),const(m,k))
            maxc = max(const(m,km1),const(m,k))
            if (minc < 0) then
               cdifr = 0.
            else
               cdifr = abs(const(m,k)-const(m,km1))/max(maxc,small)
            endif

! If the two layers differ significantly use a geometric averaging procedure
! But only do that for deep convection.  For shallow, use the simple
! averaging which is used in subr cmfmca
            if (iconvtype /= 1) then
               chat(m,k) = 0.5* (const(m,k)+const(m,km1))
            else if (cdifr > 1.E-6) then
!           if (cdifr > 1.E-6) then
               cabv = max(const(m,km1),maxc*1.e-12_r8)
               cbel = max(const(m,k),maxc*1.e-12_r8)
               chat(m,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel
            else             ! Small diff, so just arithmetic mean
               chat(m,k) = 0.5* (const(m,k)+const(m,km1))
            end if

! Set provisional up and down draft values, and tendencies
            conu(m,k) = chat(m,k)
            cond(m,k) = chat(m,k)
         end do ! k

! Values at surface inferface == values in lowest layer
         chat(m,pver+1) = const(m,pver)
         conu(m,pver+1) = const(m,pver)
         cond(m,pver+1) = const(m,pver)
      end if
      end do ! m




! Compute updraft mixing ratios from cloudbase to cloudtop
! No special treatment is needed at k=pver because arrays 
!    are dimensioned 1:pver+1
! A time-split approach is used.  First, entrainment is applied to produce 
!    an initial conu(m,k) from conu(m,k+1).  Next, chemistry/physics are
!    applied to the initial conu(m,k) to produce a final conu(m,k).
!    Detrainment from the updraft uses this final conu(m,k).
! Note that different time-split approaches would give somewhat different
!    results
      kactcnt = 0 ; kactcntb = 0 ; kactfirst = 1
k_loop_main_bb: &
      do k = kbot, ktop, -1
         kp1 = k+1

! cldfrac = conv cloud fractional area.  This could represent anvil cirrus area, 
!    and may not useful for aqueous chem and wet removal calculations
         cldfrac_i(k) = max( cldfrac_i(k), 0.005_r8 )
! mu_p_eudp(k) = updraft massflux at k, without detrainment between kp1,k
         mu_p_eudp(k) = mu_i(kp1) + eudp(k)

         fa_u(k) = 0.0_r8 !BSINGH(10/15/2014): Initialized so that it has a value if the following "if" check yeilds .false.
         if (mu_p_eudp(k) > mbsth) then
! if (mu_p_eudp(k) <= mbsth) the updraft mass flux is negligible at base and top
!    of current layer, 
! so current layer is a "gap" between two unconnected updrafts,
! so essentially skip all the updraft calculations for this layer

! First apply changes from entrainment
            f_ent = eudp(k)/mu_p_eudp(k)
            f_ent = max( 0.0_r8, min( 1.0_r8, f_ent ) )
            tmpa = 1.0 - f_ent
            do m = 2, ncnst_extd
               if (doconvproc_extd(m)) then
                  conu(m,k) = tmpa*conu(m,kp1) + f_ent*const(m,k)
               end if
            end do

! estimate updraft velocity (wup) 
            if (iconvtype /= 1) then
! shallow - wup = (mup in kg/m2/s) / [rhoair * (updraft area)]
               wup(k) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
                      / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
               wup(k) = max( 0.1_r8, wup(k) )
            else
! deep - the above method overestimates updraft area and underestimate wup
!    the following is based lemone and zipser (j atmos sci, 1980, p. 2455)
!    peak updraft (= 4 m/s) is sort of a "grand median" from their GATE data
!       and Thunderstorm Project data which they also show
!    the vertical profile shape is a crude fit to their median updraft profile
               zkm = zmagl(k)*1.0e-3
               if (zkm .ge. 1.0) then
                  wup(k) = 4.0_r8*((zkm/4.0_r8)**0.21_r8)
               else
                  wup(k) = 2.9897_r8*(zkm**0.5_r8)
               end if
               wup(k) = max( 0.1_r8, min( 4.0_r8, wup(k) ) )
            end if

! compute lagrangian transport time (dt_u) and updraft fractional area (fa_u)
! *** these must obey    dt_u(k)*mu_p_eudp(k) = dp_i(k)*fa_u(k)
            dt_u(k) = dz/wup(k)
            dt_u(k) = min( dt_u(k), dt )
            fa_u(k) = dt_u(k)*(mu_p_eudp(k)/dp_i(k))


! Now apply transformation and removal changes 
!    Skip levels where icwmr(icol,k) <= clw_cut (= 1.0e-6) to eliminate
!    occasional very small icwmr values from the ZM module
            clw_cut = 1.0e-6


            if (convproc_method_activate <= 1) then
! aerosol activation - method 1
!    skip levels that are completely glaciated (fracice(icol,k) == 1.0)
!    when kactcnt=1 (first/lowest layer with cloud water) apply 
!       activatation to the entire updraft
!    when kactcnt>1 apply activatation to the amount entrained at this level
               if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
                  kactcnt = kactcnt + 1

                  idiag_act = idiag_in(icol)
                  if ((kactcnt == 1) .or. (f_ent > 0.0_r8)) then
                     kactcntb = kactcntb + 1
                     if ((kactcntb == 1) .and. (idiag_act > 0)) then
                        write(lun,'(/a,i9,2i4)') &
                           'qaku act_conv lchnk,i,jtsub', lchnk, icol, jtsub
                     end if
                  end if

                  if (kactcnt == 1) then
                     ! diagnostic fields
                     ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
                     xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
                         / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
                     xx_kcldbase(icol) = k

                     kactfirst = k
                     tmpa = 1.0
                     call ma_activate_convproc(                            &
                        conu(:,k),  dconudt_activa(:,k), conu(:,k),        &
                        tmpa,       dt_u(k),             wup(k),           &
                        t(icol,k),  rhoair_i(k),         fracice(icol,k),  &
                        pcnst_extd, lun,                 idiag_act,        &
                        lchnk,      icol,                k,                &
                        ipass_calc_updraft                                 )
                  else if (f_ent > 0.0_r8) then
                     ! current layer is above cloud base (=first layer with activation)
                     !    only allow activation at k = kactfirst thru kactfirst-(method1_activate_nlayers-1)
                     if (k >= kactfirst-(method1_activate_nlayers-1)) then
                        call ma_activate_convproc(                            &
                           conu(:,k),  dconudt_activa(:,k), const(:,k),       &
                           f_ent,      dt_u(k),             wup(k),           &
                           t(icol,k),  rhoair_i(k),         fracice(icol,k),  &
                           pcnst_extd, lun,                 idiag_act,        &
                           lchnk,      icol,                k,                &
                           ipass_calc_updraft                                 )
                     end if
                  end if
! the following was for cam2 shallow convection (hack),
! but is not appropriate for cam5 (uwshcu)
!                 else if ((kactcnt > 0) .and. (iconvtype /= 1)) then
! !    for shallow conv, when you move from activation occuring to 
! !       not occuring, reset kactcnt=0, because the hack scheme can
! !       produce multiple "1.5 layer clouds" separated by clear air
!                    kactcnt = 0
!                 end if
               end if ! ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then

            else ! (convproc_method_activate >= 2)
! aerosol activation - method 2
!    skip levels that are completely glaciated (fracice(icol,k) == 1.0)
!    when kactcnt=1 (first/lowest layer with cloud water) 
!       apply "primary" activatation to the entire updraft
!    when kactcnt>1 
!       apply secondary activatation to the entire updraft
!       do this for all levels above cloud base (even if completely glaciated)
!          (this is something for sensitivity testing)
               do_act_this_lev = .false.
               if (kactcnt <= 0) then
                  if (icwmr(icol,k) > clw_cut) then
                     do_act_this_lev = .true.
                     kactcnt = 1
                     kactfirst = k
                     ! diagnostic fields
                     ! xx_wcldbase = w at first cloudy layer, estimated from mu and cldfrac
                     xx_wcldbase(icol) = (mu_i(kp1) + mu_i(k))*0.5_r8*hund_ovr_g &
                         / (rhoair_i(k) * (cldfrac_i(k)*0.5_r8))
                     xx_kcldbase(icol) = k
                  end if
               else
!                 if ((icwmr(icol,k) > clw_cut) .and. (fracice(icol,k) < 1.0)) then
                     do_act_this_lev = .true.
                     kactcnt = kactcnt + 1
!                 end if
               end if

               idiag_act = idiag_in(icol)
               if ( do_act_this_lev ) then
                  kactcntb = kactcntb + 1
                  if ((kactcntb == 1) .and. (idiag_act > 0)) then
                     write(lun,'(/a,i9,2i4)') &
                        'qaku act_conv lchnk,i,jtsub', lchnk, icol, jtsub
                  end if

                  call ma_activate_convproc_method2(                    &
                     conu(:,k),  dconudt_activa(:,k),                   &
                     f_ent,      dt_u(k),             wup(k),           &
                     t(icol,k),  rhoair_i(k),         fracice(icol,k),  &
                     pcnst_extd, lun,                 idiag_act,        &
                     lchnk,      icol,                k,                &
                     kactfirst,  ipass_calc_updraft                     )
               end if

            end if ! (convproc_method_activate <= 1)

! aqueous chemistry
!    do glaciated levels as aqchem_conv will eventually do acid vapor uptake
!    to ice, and aqchem_conv module checks fracice before doing liquid wtr stuff
            if (icwmr(icol,k) > clw_cut) then
!              call aqchem_conv( conu(1,k), dconudt_aqchem(1,k), aqfrac,  &
!                 t(icol,k), fracice(icol,k), icwmr(icol,k), rhoair_i(k), &
!                 lh2o2(icol,k), lo3(icol,k), dt_u(k)                     )
            end if

! wet removal
!
! mirage2
!    rprd               = precip formation as a grid-cell average (kgW/kgA/s)
!    icwmr              = cloud water MR within updraft area (kgW/kgA)
!    fupdr              = updraft fractional area (--)
!    A = rprd/fupdr     = precip formation rate within updraft area (kgW/kgA/s)
!    B = A/icwmr = rprd/(icwmr*fupdr) 
!                       = first-order removal rate (1/s)
!    C = dp/(mup/fupdr) = updraft air residence time in the layer (s)
!
!    fraction removed = (1.0 - exp(-cdt)) where
!                 cdt = B*C = (dp/mup)*rprd/icwmr
!
!    Note1:  fupdr cancels out in cdt, so need not be specified
!    Note2:  dp & mup units need only be consistent (e.g., mb & mb/s)
!    Note3:  for shallow conv, cdt = 1-beta (beta defined in Hack scheme)
!    Note4:  the "dp" in C above and code below should be the moist dp
!
! cam5
!    clw_preloss = cloud water MR before loss to precip
!                = icwmr + dt*(rprd/fupdr)
!    B = A/clw_preloss  = (rprd/fupdr)/(icwmr + dt*rprd/fupdr) 
!                       = rprd/(fupdr*icwmr + dt*rprd) 
!                       = first-order removal rate (1/s)
!
!    fraction removed = (1.0 - exp(-cdt)) where
!                 cdt = B*C = (fupdr*dp/mup)*[rprd/(fupdr*icwmr + dt*rprd)]
!
!    Note1:  *** cdt is now sensitive to fupdr, which we do not really know,
!                and is not the same as the convective cloud fraction
!    Note2:  dt is appropriate in the above cdt expression, not dtsub
!
!    Apply wet removal at levels where 
!       icwmr(icol,k) > clw_cut  AND  rprd(icol,k) > 0.0
!    as wet removal occurs in both liquid and ice clouds
!
            cdt(k) = 0.0_r8
            if ((icwmr(icol,k) > clw_cut) .and. (rprd(icol,k) > 0.0)) then 
!              if (iconvtype == 1) then
                  tmpf = 0.5_r8*cldfrac_i(k)
                  cdt(k) = sfadj * (tmpf*dp(i,k)/mu_p_eudp(k)) * rprd(icol,k) / &
                        (tmpf*icwmr(icol,k) + dt*rprd(icol,k))
!              else if (k < pver) then
!                 if (eudp(k+1) > 0) cdt(k) =   &
!                       rprd(icol,k)*dp(i,k)/(icwmr(icol,k)*eudp(k+1))
!              end if
            end if
            if (cdt(k) > 0.0_r8) then
               expcdtm1 = exp(-cdt(k)) - 1.0
               do m = 2, ncnst_extd
                  if (doconvproc_extd(m)) then
                     dconudt_wetdep(m,k) = conu(m,k)*aqfrac(m)*expcdtm1
                     conu(m,k) = conu(m,k) + dconudt_wetdep(m,k)
                     dconudt_wetdep(m,k) = dconudt_wetdep(m,k) / dt_u(k)
                  end if
               enddo
            end if

         end if    ! "(mu_p_eudp(k) > mbsth)"
      end do k_loop_main_bb ! "k = kbot, ktop, -1"

! when doing updraft calcs twice, only need to go this far on the first pass
      if ( (ipass_calc_updraft == 1) .and. &
           (npass_calc_updraft == 2) ) cycle ipass_calc_updraft_loop

      if (idiag_in(icol) > 0) then
         ! do wet removal diagnostics here 
         do k = kbot, ktop, -1
            if (mu_p_eudp(k) > mbsth) &
               write(lun,'(a,i9,3i4,1p,6e10.3)') &
                  'qakr - l,i,k,jt; cdt, cldfrac, icwmr, rprd, ...', lchnk, icol, k, jtsub, &
                  cdt(k), cldfrac_i(k), icwmr(icol,k), rprd(icol,k), dp(i,k), mu_p_eudp(k)
         end do
      end if


! Compute downdraft mixing ratios from cloudtop to cloudbase
! No special treatment is needed at k=2
! No transformation or removal is applied in the downdraft
      do k = ktop, kbot
         kp1 = k + 1
! md_m_eddp = downdraft massflux at kp1, without detrainment between k,kp1
         md_m_eddp = md_i(k) - eddp(k)
         if (md_m_eddp < -mbsth) then
            do m = 2, ncnst_extd
               if (doconvproc_extd(m)) then
                  cond(m,kp1) = ( md_i(k)*cond(m,k)		&
                                - eddp(k)*const(m,k) ) / md_m_eddp
               endif
            end do
         end if
      end do ! k


! Now computes fluxes and tendencies
! NOTE:  The approach used in convtran applies to inert tracers and
!        must be modified to include source and sink terms
      sumflux(:) = 0.0
      sumflux2(:) = 0.0
      sumsrce(:) = 0.0
      sumchng(:) = 0.0
      sumchng3(:) = 0.0
      sumactiva(:) = 0.0
      sumaqchem(:) = 0.0
      sumwetdep(:) = 0.0
      sumresusp(:) = 0.0
      sumprevap(:) = 0.0
      sumprevap_hist(:) = 0.0 ! REASTER 08/05/2015

      maxflux(:) = 0.0
      maxflux2(:) = 0.0
      maxresusp(:) = 0.0
      maxsrce(:) = 0.0
      maxprevap(:) = 0.0

k_loop_main_cc: &
      do k = ktop, kbot
         kp1 = k+1
         km1 = k-1
         kp1x = min( kp1, pver )
         km1x = max( km1, 1 )
         fa_u_dp = fa_u(k)*dp_i(k)
         do m = 2, ncnst_extd
            if (doconvproc_extd(m)) then

! First compute fluxes using environment subsidence/lifting and 
! entrainment/detrainment into up/downdrafts, 
! to provide an additional mass balance check
! (this could be deleted after the code is well tested)
               fluxin  = mu_i(k)*min(chat(m,k),const(m,km1x))       &
                       - md_i(kp1)*min(chat(m,kp1),const(m,kp1x))   &
                       + dudp(k)*conu(m,k) + dddp(k)*cond(m,kp1)
               fluxout = mu_i(kp1)*min(chat(m,kp1),const(m,k))      &
                       - md_i(k)*min(chat(m,k),const(m,k))          &
                       + (eudp(k) + eddp(k))*const(m,k)

               netflux = fluxin - fluxout

               sumflux2(m) = sumflux2(m) + netflux
               maxflux2(m) = max( maxflux2(m), abs(fluxin), abs(fluxout) )

! Now compute fluxes as in convtran, and also source/sink terms
! (version 3 limit fluxes outside convection to mass in appropriate layer
! (these limiters are probably only safe for positive definite quantitities
! (it assumes that mu and md already satify a courant number limit of 1)
            if (iflux_method /= 2) then
               fluxin  =     mu_i(kp1)*conu(m,kp1)                     &
                           + mu_i(k  )*min(chat(m,k  ),const(m,km1x))  &
                         - ( md_i(k  )*cond(m,k)                       &
                           + md_i(kp1)*min(chat(m,kp1),const(m,kp1x)) )
               fluxout =     mu_i(k  )*conu(m,k)                       &
                           + mu_i(kp1)*min(chat(m,kp1),const(m,k   ))  &
                         - ( md_i(kp1)*cond(m,kp1)                     &
                           + md_i(k  )*min(chat(m,k  ),const(m,k   )) )
            else
               fluxin  =     mu_i(kp1)*conu(m,kp1)                     &
                         - ( md_i(k  )*cond(m,k) )
               fluxout =     mu_i(k  )*conu(m,k)                       &
                         - ( md_i(kp1)*cond(m,kp1) )
               tmpveca(1) = fluxin ; tmpveca(4) = -fluxout

               ! new method -- simple upstream method for the env subsidence
               ! tmpa = net env mass flux (positive up) at top of layer k
               tmpa = -( mu_i(k  ) + md_i(k  ) )
               if (tmpa <= 0.0_r8) then
                  fluxin  = fluxin  - tmpa*const(m,km1x)
               else
                  fluxout = fluxout + tmpa*const(m,k   )
               end if
               tmpveca(2) = fluxin ; tmpveca(5) = -fluxout
               ! tmpa = net env mass flux (positive up) at base of layer k
               tmpa = -( mu_i(kp1) + md_i(kp1) )
               if (tmpa >= 0.0_r8) then
                  fluxin  = fluxin  + tmpa*const(m,kp1x)
               else
                  fluxout = fluxout - tmpa*const(m,k   )
               end if
               tmpveca(3) = fluxin ; tmpveca(6) = -fluxout
            end if

            netflux = fluxin - fluxout
            netsrce = fa_u_dp*(dconudt_aqchem(m,k) + &
                        dconudt_activa(m,k) + dconudt_wetdep(m,k))
            dcondt(m,k) = (netflux+netsrce)/dp_i(k)

            dcondt_wetdep(m,k) = fa_u_dp*dconudt_wetdep(m,k)/dp_i(k)

            sumflux(m) = sumflux(m) + netflux
            maxflux(m) = max( maxflux(m), abs(fluxin), abs(fluxout) )
            sumsrce(m) = sumsrce(m) + netsrce
            maxsrce(m) = max( maxsrce(m), &
               fa_u_dp*max( abs(dconudt_aqchem(m,k)), &
                            abs(dconudt_activa(m,k)),  abs(dconudt_wetdep(m,k)) ) )
            sumchng(m) = sumchng(m) + dcondt(m,k)*dp_i(k)
            sumactiva(m) = sumactiva(m) + fa_u_dp*dconudt_activa(m,k)
            sumaqchem(m) = sumaqchem(m) + fa_u_dp*dconudt_aqchem(m,k)
            sumwetdep(m) = sumwetdep(m) + fa_u_dp*dconudt_wetdep(m,k)

            if ( idiag_in(icol)>0 .and. k==26 .and. &
                 (m==16 .or. m==23 .or. m==16+pcnst .or. m==23+pcnst) ) then
               if (m==16) &
               write(lun,'(a,i9,4i4,1p,22x,   2x,11x,  2x,6e11.3)') &
                  'qakww0-'//convtype(1:4), lchnk, icol, k, -1, jtsub, &
                  dtsub*mu_i(k+1)/dp_i(k), dtsub*mu_i(k)/dp_i(k), dtsub*eudp(k)/dp_i(k), &
                  dtsub*md_i(k+1)/dp_i(k), dtsub*md_i(k)/dp_i(k), dtsub*eddp(k)/dp_i(k) 

               write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,6e11.3)') &
                  'qakww1-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
                  const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k), &
                  dtsub*fluxin/dp_i(k), -dtsub*fluxout/dp_i(k), &
                  dtsub*fa_u_dp*dconudt_aqchem(m,k)/dp_i(k), &
                  dtsub*fa_u_dp*dconudt_activa(m,k)/dp_i(k), &
                  dtsub*fa_u_dp*dconudt_wetdep(m,k)/dp_i(k)
               write(lun,'(a,i9,4i4,1p,22x,   2x,11x,  2x,6e11.3)') &
                  'qakww1-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
                  dtsub*tmpveca(1:6)/dp_i(k)
            end if

            end if   ! "(doconvproc_extd(m))"
         end do      ! "m = 2,ncnst_extd"
      end do k_loop_main_cc ! "k = ktop, kbot"


! calculate effects of precipitation evaporation
      call ma_precpevap_convproc( dcondt, dcondt_wetdep,  dcondt_prevap,   &
                                  dcondt_prevap_hist,                      & ! REASTER 08/05/2015
                                  rprd,   evapc,          dp_i,            &
                                  icol,   ktop,           pcnst_extd,      &
                                  lun,    idiag_in(icol), lchnk,           &
                                  doconvproc_extd,                         &
                                  species_class, mam_prevap_resusp_optaa   ) ! REASTER 08/05/2015
      if ( idiag_in(icol)>0 ) then
         k = 26
         do m = 16, 23, 7
            write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
               'qakww2-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
               const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
         end do
         do m = 16+pcnst, 23+pcnst, 7
            write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
               'qakww2-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
               const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
         end do
      end if



! make adjustments to dcondt for activated & unactivated aerosol species
!    pairs to account any (or total) resuspension of convective-cloudborne aerosol
      call ma_resuspend_convproc( dcondt, dcondt_resusp,   &
                                  const, dp_i, ktop, kbot_prevap, pcnst_extd ) ! REASTER 08/05/2015
      if ( idiag_in(icol)>0 ) then
         k = 26
         do m = 16, 23, 7
            write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
               'qakww3-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
               const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
         end do
         do m = 16+pcnst, 23+pcnst, 7
            write(lun,'(a,i9,4i4,1p,2e11.3,2x,e11.3,2x,5e11.3)') &
               'qakww3-'//convtype(1:4), lchnk, icol, k, m, jtsub, &
               const(m,k), const(m,k)+dtsub*dcondt(m,k), dtsub*dcondt(m,k)
         end do
      end if


! calculate new column-tendency variables
      do m = 2, ncnst_extd
         if (doconvproc_extd(m)) then
            ! should go to k=pver for dcondt_prevap, and this should be safe for other sums
            do k = ktop, kbot_prevap  ! REASTER 08/05/2015
               sumchng3(m)  = sumchng3(m)  + dcondt(m,k)*dp_i(k)
               sumresusp(m) = sumresusp(m) + dcondt_resusp(m,k)*dp_i(k)
               maxresusp(m) = max( maxresusp(m),   &
                                         abs(dcondt_resusp(m,k)*dp_i(k)) )
               sumprevap(m) = sumprevap(m) + dcondt_prevap(m,k)*dp_i(k)
               maxprevap(m) = max( maxprevap(m),   &
                                         abs(dcondt_prevap(m,k)*dp_i(k)) )
               sumprevap_hist(m) = sumprevap_hist(m) + dcondt_prevap_hist(m,k)*dp_i(k) ! REASTER 08/05/2015
            end do
         end if
      end do ! m


! do checks for mass conservation
! do not expect errors > 1.0e-14, but use a conservative 1.0e-10 here,
! as an error of this size is still not a big concern
      relerr_cut = 1.0e-10
      if (nerr < nerrmax) then
         merr = 0
         if (courantmax > (1.0_r8 + 1.0e-6_r8)) then
            write(lun,9161) '-', trim(convtype), courantmax
            merr = merr + 1
         end if
         do m = 2, ncnst_extd
            if (doconvproc_extd(m)) then
               itmpa = 0
               ! sumflux should be ~=0.0 because fluxout of one layer cancels
               !    fluxin to adjacent layer
               tmpa = sumflux(m)
               tmpb = max( maxflux(m), small )
               if (abs(tmpa) > relerr_cut*tmpb) then
                  write(lun,9151) '1', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
                  itmpa = itmpa + 1
               end if
               ! sumflux2 involve environment fluxes and entrainment/detrainment
               !    to up/downdrafts, and it should be equal to sumchng,
               !    and so (sumflux2 - sumsrce) should be ~=0.0
               tmpa = sumflux2(m) - sumsrce(m)
               tmpb = max( maxflux2(m), maxsrce(m), small )
               if (abs(tmpa) > relerr_cut*tmpb) then
                  write(lun,9151) '2', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
                  itmpa = itmpa + 10
               end if
               ! sunchng = sumflux + sumsrce, so (sumchng - sumsrc) should be ~=0.0
               tmpa = sumchng(m) - sumsrce(m)
               tmpb = max( maxflux(m), maxsrce(m), small )
               if (abs(tmpa) > relerr_cut*tmpb) then
                  write(lun,9151) '3', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
                  itmpa = itmpa + 100
               end if
               ! sumchng3 = sumchng + sumresusp + sumprevap, 
               !    so tmpa (below) should be ~=0.0
               tmpa = sumchng3(m) - (sumsrce(m) + sumresusp(m) + sumprevap(m))
               tmpb = max( maxflux(m), maxsrce(m), maxresusp(m), maxprevap(m), small )
               if (abs(tmpa) > relerr_cut*tmpb) then
                  write(lun,9151) '4', m, cnst_name_extd(m), tmpb, tmpa, (tmpa/tmpb)
                  itmpa = itmpa + 1000
               end if

               if (itmpa > 0) merr = merr + 1
            end if
         end do ! m
         if (merr > 0) write(lun,9181) convtype, lchnk, icol, i, jt(i), mx(i)
         nerr = nerr + merr
         if (nerr >= nerrmax) write(lun,9171) nerr
      end if ! (nerr < nerrmax) then

9151 format( '*** ma_convproc_tend error, massbal', a, 1x, i5,1x,a, &
             ' -- maxflux, sumflux, relerr =', 3(1pe14.6) )
9161 format( '*** ma_convproc_tend error, courantmax', 2a, 3x, 1pe14.6 )
9171 format( '*** ma_convproc_tend error, stopping messages after nerr =', i10 )

9181 format( '*** ma_convproc_tend error -- convtype, lchnk, icol, il, jt, mx =  ', a,2x,5(1x,i10) )


!
! note again the ma_convproc_tend does not apply convective cloud processing
!    to the stratiform-cloudborne aerosol
! within this routine, cloudborne aerosols are convective-cloudborne
!
! before tendencies (dcondt, which is loaded into dqdt) are returned,
!    the convective-cloudborne aerosol tendencies must be combined
!    with the interstitial tendencies
! ma_resuspend_convproc has already done this for the dcondt
!
! the individual process column tendencies (sumwetdep, sumprevap, ...)
!    are just diagnostic fields that can be written to history
! tendencies for interstitial and convective-cloudborne aerosol could
!    both be passed back and output, if desired
! currently, however, the interstitial and convective-cloudborne tendencies
!    are combined (in the next code block) before being passed back (in qsrflx)
!
      do n = 1, ntot_amode
         do ll = 0, nspec_amode(n)
            if (ll == 0) then
               la = numptr_amode(n)
               lc = numptrcw_amode(n) + pcnst
            else
               la = lmassptr_amode(ll,n)
               lc = lmassptrcw_amode(ll,n) + pcnst
            end if
            if (doconvproc(la)) then
               sumactiva(la) = sumactiva(la) + sumactiva(lc)
               sumresusp(la) = sumresusp(la) + sumresusp(lc)
               sumaqchem(la) = sumaqchem(la) + sumaqchem(lc)
               sumwetdep(la) = sumwetdep(la) + sumwetdep(lc)
               sumprevap(la) = sumprevap(la) + sumprevap(lc)
               sumprevap_hist(la) = sumprevap_hist(la) + sumprevap_hist(lc) ! REASTER 08/05/2015
!              if (n==1 .and. ll==1) then
!	           write(lun,*) 'la, sumaqchem(la) =', la, sumaqchem(la)
!              endif
            end if
         enddo ! ll
      enddo ! n

!
! scatter overall tendency back to full array
!
      do m = 2, ncnst
         if (doconvproc(m)) then
            do k = ktop, kbot_prevap  ! should go to k=pver because of prevap ! REASTER 08/05/2015
               dqdt_i(k,m) = dcondt(m,k)
               dqdt(icol,k,m) = dqdt(icol,k,m) + dqdt_i(k,m)*xinv_ntsub
            end do
!           dqdt_i(:,m) = 0.
         end if
      end do ! m

! scatter column burden tendencies for various processes to qsrflx
      do m = 2, ncnst
         if (doconvproc(m)) then
            qsrflx_i(m,1) = sumactiva(m)*hund_ovr_g
            qsrflx_i(m,2) = sumresusp(m)*hund_ovr_g
            qsrflx_i(m,3) = sumaqchem(m)*hund_ovr_g
            qsrflx_i(m,4) = sumwetdep(m)*hund_ovr_g
            qsrflx_i(m,5) = sumprevap(m)*hund_ovr_g
            qsrflx_i(m,6) = sumprevap_hist(m)*hund_ovr_g ! REASTER 08/05/2015
!           qsrflx_i(m,1:4) = 0.
            qsrflx(icol,m,1:6) = qsrflx(icol,m,1:6) + qsrflx_i(m,1:6)*xinv_ntsub ! REASTER 08/05/2015
         end if
      end do ! m


! diagnostic output of profiles before
      if (idiag_in(icol) > 0) then
         write(lun, '(/3a,i9,2i4)' ) 'qakr-', trim(convtype), ' - lchnk,i,jtsub', lchnk, icol, jtsub
         n = 1

         do j = 1, 2
            if (j == 1) then
               write(lun, '(4a,i4)' ) &
                  'qakr-', trim(convtype), ' - k, mu,md; then mode-1 ', &
                  'numb & numbcw for q, const, conu, cond, delq(a/c/ac noresu)', jtsub
            else
               write(lun, '(/4a,i4)' ) &
                  'qakr-', trim(convtype), ' - k, mu,md; then mode-1 ', &
                  'mass & masscw for q, const, conu, cond, delq(a/c/ac noresu)', jtsub
            end if

            do k = 10, pver
               tmpveca(:) = 0.0
               do ll = 1, nspec_amode(n)
                  if (j == 1) then
                     la = numptr_amode(n)
                     lc = numptr_amode(n) + pcnst
                  else
                     la = lmassptr_amode(ll,n)
                     lc = lmassptr_amode(ll,n) + pcnst
                  end if
                  tmpveca(1)  = tmpveca(1) + q_i(k,la)
                  tmpveca(2)  = tmpveca(2) + const(la,k)
                  tmpveca(3)  = tmpveca(3) + const(lc,k)
                  tmpveca(4)  = tmpveca(4) + conu( la,k)
                  tmpveca(5)  = tmpveca(5) + conu( lc,k)
                  tmpveca(6)  = tmpveca(6) + cond( la,k)
                  tmpveca(7)  = tmpveca(7) + cond( lc,k)
                  tmpveca(8)  = tmpveca(8) + (dcondt(la,k)-dcondt_resusp(la,k))*dtsub
                  tmpveca(9)  = tmpveca(9) + (dcondt(lc,k)-dcondt_resusp(lc,k))*dtsub
                  tmpveca(10) = tmpveca(8) + tmpveca(9)
                  if (j == 1) exit
               end do ! ll
               if ((k > 15) .and. (mod(k,5) == 1)) write(lun,'(a)')
               write(lun, '(a,i3,1p,2e10.2, e11.2, 3(2x,2e9.2), 2x,3e10.2 )' ) 'qakr', k, &
                  mu_i(k), md_i(k), tmpveca(1:10)
            end do ! k
         end do ! j

         if (pcnst < 0) then
         write(lun, '(/a,i4)' ) &
            'qakr - name; burden; qsrflx tot, activa,resusp,aqchem,wetdep,resid', jtsub
         do m = 2, ncnst
            if ( .not. doconvproc(m) ) cycle
            tmpveca(1) = sum(    q_i(:,m)*dp_i(:) ) * hund_ovr_g
            tmpveca(2) = sum( dqdt_i(:,m)*dp_i(:) ) * hund_ovr_g
            tmpveca(3:6) = qsrflx_i(m,1:4)
            tmpveca(7) = tmpveca(2) - sum( tmpveca(3:6) )
            write(lun, '(2a,1p,2(2x,e11.3),2x,4e11.3,2x,e11.3)' ) &
               'qakr  ', cnst_name_extd(m)(1:10), tmpveca(1:7)
         end do ! m
         end if ! (pcnst < 0) then

         write(lun, '(/3a,i4)' ) 'qakr-', trim(convtype), &
            ' - name; burden; sumchng3, sumactiva,resusp,aqchem,wetdep, resid,resid*dt/burden', jtsub
!        write(lun, '(/2a)' ) &
!           'qakr - name;  burden;  sumchng3;  ', &
!           'sumactiva,resusp,aqchem,wetdep,prevap;  resid,resid*dtsub/burden'
         tmpb = 0.0_r8
         itmpb = 0
         do m = 2, pcnst
            if ( .not. doconvproc_extd(m) ) cycle

            tmpmata(:,:) = 0.0
            do j = 1, 3
               l = m
               if (j == 3) l = m + pcnst
               if ( .not. doconvproc_extd(l) ) cycle

               if (j == 1) then
                  tmpmata(1,j) = sum(    q_i(:,l)*dp_i(:) ) * hund_ovr_g
                  tmpmata(2,j) = sum( dqdt_i(:,l)*dp_i(:) ) * hund_ovr_g
                  tmpmata(3:7,j) = qsrflx_i(l,1:5)
               else
                  tmpmata(1,j) = sum( const(l,1:pver)*dp_i(1:pver) ) * hund_ovr_g
                  tmpmata(2,j) = sumchng3( l) * hund_ovr_g
                  tmpmata(3,j) = sumactiva(l) * hund_ovr_g
                  tmpmata(4,j) = sumresusp(l) * hund_ovr_g
                  tmpmata(5,j) = sumaqchem(l) * hund_ovr_g
                  tmpmata(6,j) = sumwetdep(l) * hund_ovr_g
                  tmpmata(7,j) = sumprevap(l) * hund_ovr_g
               end if
            end do ! j

            tmpmata(3:7,2) = tmpmata(3:7,2) - tmpmata(3:7,3)  ! because lc values were added onto la
            do j = 1, 3
               tmpmata(8,j) = tmpmata(2,j) - sum( tmpmata(3:7,j) )  ! residual
               tmpa = max( tmpmata(1,min(j,2)), 1.0e-20_r8 )
               tmpmata(9,j) = tmpmata(8,j) * dtsub / tmpa
               if (abs(tmpmata(9,j)) > tmpb) then
                  tmpb = abs(tmpmata(9,j))
                  itmpb = m
               end if
            end do

!           write(lun, '(/2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
!              'qakr1 ', cnst_name_extd(m)(1:10), tmpmata(1:6,1), tmpmata(8:9,1)
            write(lun, '(/2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
               'qakr1 ', cnst_name_extd(m)(1:10), tmpmata(1:9,1)
!           write(lun, '( 2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
!              'qakr2 ', cnst_name_extd(m)(1:10), tmpmata(1:6,2), tmpmata(8:9,2)
            write(lun, '( 2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
               'qakr2 ', cnst_name_extd(m)(1:10), tmpmata(1:9,2)
            if ( .not. doconvproc_extd(l) ) cycle
!           write(lun, '( 2a,1p,2(2x,e11.3),2x,4e11.3,2x,2e11.3)' ) &
!              'qakr3 ', cnst_name_cw(m)(1:10), tmpmata(1:6,3), tmpmata(8:9,3)
            write(lun, '( 2a,1p,2(2x,e11.3),2x,5e11.3,2x,2e11.3)' ) &
               'qakr3 ', cnst_name_cw(m)(1:10), tmpmata(1:9,3)
         end do ! m
         write(lun, '(/3a,2i4,1p,e11.2)' ) 'qakr-', trim(convtype), &
            ' - max(resid*dt/burden)', jtsub, itmpb, tmpb

      end if ! (idiag_in(icol) > 0) then


      if (jtsub < ntsub) then
         ! update the q_i for the next interation of the jtsub loop
         do m = 2, ncnst
            if (doconvproc(m)) then
               do k = ktop, kbot_prevap  ! should go to k=pver because of prevap ! REASTER 08/05/2015
                  q_i(k,m) = max( (q_i(k,m) + dqdt_i(k,m)*dtsub), 0.0_r8 )
               end do
            end if
         end do ! m
      end if

      end do ipass_calc_updraft_loop

      end do jtsub_loop_main_aa  ! of the main "do jtsub = 1, ntsub" loop


   end do i_loop_main_aa  ! of the main "do i = il1g, il2g" loop

   return
end subroutine ma_convproc_tend



!=========================================================================================
   subroutine ma_precpevap_convproc(                           &
              dcondt,  dcondt_wetdep, dcondt_prevap,           &
              dcondt_prevap_hist,                              & ! REASTER 08/05/2015
              rprd,    evapc,         dp_i,                    &
              icol,    ktop,          pcnst_extd,              &
              lun,     idiag_prevap,  lchnk,                   &
              doconvproc_extd,                                 &
              species_class, mam_prevap_resusp_optaa           ) ! REASTER 08/05/2015
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate resuspension of wet-removed aerosol species resulting 
!    precip evaporation
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

   use ppgrid, only: pcols, pver
   use constituents, only: pcnst
   use physconst, only: spec_class_aerosol ! REASTER 08/05/2015

   use modal_aero_data, only:  &
      lmassptrcw_amode, nspec_amode, numptrcw_amode, &
      mmtoo_prevap_resusp ! REASTER 08/05/2015

   implicit none

!-----------------------------------------------------------------------
! arguments
! (note:  TMR = tracer mixing ratio)
   integer,  intent(in)    :: pcnst_extd

   real(r8), intent(inout) :: dcondt(pcnst_extd,pver)
                              ! overall TMR tendency from convection
   real(r8), intent(in)    :: dcondt_wetdep(pcnst_extd,pver)
                              ! portion of TMR tendency due to wet removal
   real(r8), intent(inout) :: dcondt_prevap(pcnst_extd,pver)
                              ! portion of TMR tendency due to precip evaporation
                              ! (actually, due to the adjustments made here)
                              ! (on entry, this is 0.0)
   real(r8), intent(inout) :: dcondt_prevap_hist(pcnst_extd,pver) ! REASTER 08/05/2015
                              ! this determines what goes into the history 
                              !    precip-evap SFSEC variables
                              ! currently, the SFSEC resuspension are attributed
                              !    to the species that got scavenged,
                              !    WHICH IS NOT the species that actually
                              !    receives the resuspension
                              !    when modal_aero_wetdep_resusp_opt > 0
                              ! so when scavenged so4_c1 is resuspended as so4_a1, 
                              !    this resuspension column-tendency shows
                              !    up in so4_c1SFSES
                              ! this is done to allow better tracking of the
                              !    resuspension in the mass-budget post-processing scripts

   real(r8), intent(in)    :: rprd(pcols,pver)  ! conv precip production  rate (gathered)
   real(r8), intent(in)    :: evapc(pcols,pver)  ! conv precip evaporation rate (gathered)
   real(r8), intent(in)    :: dp_i(pver) ! pressure thickness of level (in mb)

   integer,  intent(in)    :: icol  ! normal (ungathered) i index for current column
   integer,  intent(in)    :: ktop  ! index of top cloud level for current column
   integer,  intent(in)    :: lun    ! logical unit for diagnostic output
   integer,  intent(in)    :: idiag_prevap ! flag for diagnostic output
   integer,  intent(in)    :: lchnk  ! chunk index

   logical,  intent(in)    :: doconvproc_extd(pcnst_extd)  ! indicates which species to process
   integer,  intent(in)    :: species_class(:) ! REASTER 08/05/2015
   integer,  intent(in)    :: mam_prevap_resusp_optaa ! REASTER 08/05/2015

!-----------------------------------------------------------------------
! local variables
   integer  :: k, l, ll, m, m2, mmtoo, n ! REASTER 08/05/2015
   real(r8) :: del_pr_flux_prod      ! change to precip flux from production  [(kg/kg/s)*mb]
   real(r8) :: del_pr_flux_evap      ! change to precip flux from evaporation [(kg/kg/s)*mb]
   real(r8) :: del_wd_flux_evap      ! change to wet deposition flux from evaporation [(kg/kg/s)*mb]
   real(r8) :: fdel_pr_flux_evap     ! fractional change to precip flux from evaporation
   real(r8) :: pr_flux               ! precip flux at base of current layer [(kg/kg/s)*mb]
   real(r8) :: pr_flux_old
   real(r8) :: tmpa, tmpb, tmpc, tmpd
   real(r8) :: tmpdp                 ! delta-pressure (mb)
   real(r8) :: wd_flux(pcnst_extd)   ! tracer wet deposition flux at base of current layer [(kg/kg/s)*mb]
   character(len=100) :: msg
   
!-----------------------------------------------------------------------


   if ( mam_prevap_resusp_optaa == 30 ) then
      call ma_precpevap30_convproc(                         &
           dcondt,  dcondt_wetdep, dcondt_prevap,           &
           dcondt_prevap_hist,                              &
           rprd,    evapc,         dp_i,                    &
           icol,    ktop,          pcnst_extd,              &
           lun,     idiag_prevap,  lchnk,                   &
           doconvproc_extd,                                 &
           species_class, mam_prevap_resusp_optaa           )
      return
   else if ( mam_prevap_resusp_optaa /=   0 .and. &
             mam_prevap_resusp_optaa /=  10 .and. &
             mam_prevap_resusp_optaa /=  11 .and. &
             mam_prevap_resusp_optaa /=  20 .and. &
             mam_prevap_resusp_optaa /=  21 ) then
      write(msg,'(a,2(1x,i10))') &
         'ma_precpevap_convproc - bad mam_prevap_resusp_optaa =', &
         mam_prevap_resusp_optaa
      call endrun( msg )
   end if

!
! *** note use of non-standard units
!
! precip
!    tmpdp = dp_i is mb
!    rprd and evapc are kgwtr/kgair/s
!    pr_flux = tmpdp*rprd is mb*kgwtr/kgair/s
! this works ok because the only important thing is fdel_pr_flux_evap which is dimensionless
!
! precip-borne aerosol
!    dcondt_wetdep is kgaero/kgair/s
!    wd_flux = tmpdp*dcondt_wetdep is mb*kgaero/kgair/s
!    dcondt_prevap = del_wd_flux_evap/tmpdp is kgaero/kgair/s
! so this works ok too
!
   pr_flux = 0.0_r8
   wd_flux(:) = 0.0_r8

   if (idiag_prevap > 0) then
      write(lun,'(a,i9,i4,5x,a)') 'qakx - lchnk,i', lchnk, icol, &
         '// k; pr_flux old,new; delprod,devap; mode-1 numb wetdep,prevap; mass ...'
   end if

   do k = ktop, pver
      tmpdp = dp_i(k)

      pr_flux_old = pr_flux
      del_pr_flux_prod = tmpdp*max(0.0_r8, rprd(icol,k))
      pr_flux = pr_flux_old + del_pr_flux_prod

      del_pr_flux_evap = min( pr_flux, tmpdp*max(0.0_r8, evapc(icol,k)) )
      fdel_pr_flux_evap = del_pr_flux_evap / max(pr_flux, 1.0e-35_r8)

      if ( mam_prevap_resusp_optaa <= 0 ) then
         fdel_pr_flux_evap = 0.0_r8  ! REASTER 08/05/2015 - turn off resuspension from precip evap
      end if

      do m = 2, pcnst_extd
         if ( doconvproc_extd(m) ) then

            ! dcondt_wetdep(m,k) is negative (or zero), so use -dcondt_wetdep(m,k) here
            wd_flux(m) = wd_flux(m) + tmpdp*max(0.0_r8, -dcondt_wetdep(m,k)) 
            del_wd_flux_evap = wd_flux(m)*fdel_pr_flux_evap

! REASTER 08/05/2015 BEGIN
!           wd_flux(m) = max( 0.0_r8, wd_flux(m)-del_wd_flux_evap )

!           dcondt_prevap(m,k) = del_wd_flux_evap/tmpdp
!           dcondt(m,k) = dcondt(m,k) + dcondt_prevap(m,k)

            ! do this for mam_prevap_resusp_optaa = 0,10,11,20,21
            ! (also for trace gases when mam_prevap_resusp_opt = ???,
            !  but currently modal_aero_convproc does not do trace gases)
            dcondt_prevap_hist(m,k) = del_wd_flux_evap/tmpdp
            dcondt_prevap(m,k) = del_wd_flux_evap/tmpdp
            dcondt(m,k) = dcondt(m,k) + dcondt_prevap(m,k)

            wd_flux(m) = max( 0.0_r8, wd_flux(m)-del_wd_flux_evap )

         end if ! ( doconvproc_extd(m) ) then
      end do ! m = 2, pcnst_extd
! REASTER 08/05/2015 END

      pr_flux = max( 0.0_r8, pr_flux-del_pr_flux_evap )

      if (idiag_prevap > 0) then
         n = 1
         l = numptrcw_amode(n) + pcnst
         tmpa = dcondt_wetdep(l,k)
         tmpb = dcondt_prevap(l,k)
         tmpc = 0.0_r8
         tmpd = 0.0_r8
         do ll = 1, nspec_amode(n)
            l = lmassptrcw_amode(ll,n) + pcnst
            tmpc = tmpc + dcondt_wetdep(l,k)
            tmpd = tmpd + dcondt_prevap(l,k)
         end do
         write(lun,'(a,i4,1p,4(2x,2e10.2))') 'qakx', k, &
            pr_flux_old, pr_flux, del_pr_flux_prod, -del_pr_flux_evap, &
            -tmpa, tmpb, -tmpc, tmpd
      end if
   end do ! k

   return
   end subroutine ma_precpevap_convproc



!=========================================================================================
   subroutine ma_precpevap30_convproc(                         &
              dcondt,  dcondt_wetdep, dcondt_prevap,           &
              dcondt_prevap_hist,                              & ! REASTER 08/05/2015
              rprd,    evapc,         dp_i,                    &
              icol,    ktop,          pcnst_extd,              &
              lun,     idiag_prevap,  lchnk,                   &
              doconvproc_extd,                                 &
              species_class, mam_prevap_resusp_optaa           ) ! REASTER 08/05/2015
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate resuspension of wet-removed aerosol species resulting precip evaporation
! for mam_prevap_resusp_optaa = 30
!
!     for aerosol mass   species, do non-linear resuspension to coarse mode
!     for aerosol number species, all the resuspension is done in wetdepa_v2, so do nothing here
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

   use ppgrid, only: pcols, pver
   use constituents, only: pcnst
   use physconst, only: pi, spec_class_aerosol ! REASTER 08/05/2015

   use modal_aero_data, only:  &
      lmassptr_amode, lmassptrcw_amode, lspectype_amode, &
      nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, &
      mmtoo_prevap_resusp, ntoo_prevap_resusp, & ! REASTER 08/05/2015
      specdens_amode
   use wetdep, only:  faer_resusp_vs_fprec_evap_mpln

   implicit none

!-----------------------------------------------------------------------
! arguments
! (note:  TMR = tracer mixing ratio)
   integer,  intent(in)    :: pcnst_extd

   real(r8), intent(inout) :: dcondt(pcnst_extd,pver)
                              ! overall TMR tendency from convection
   real(r8), intent(in)    :: dcondt_wetdep(pcnst_extd,pver)
                              ! portion of TMR tendency due to wet removal
   real(r8), intent(inout) :: dcondt_prevap(pcnst_extd,pver)
                              ! portion of TMR tendency due to precip evaporation
                              ! (actually, due to the adjustments made here)
                              ! (on entry, this is 0.0)
   real(r8), intent(inout) :: dcondt_prevap_hist(pcnst_extd,pver) ! REASTER 08/05/2015
                              ! this determines what goes into the history 
                              !    precip-evap SFSEC variables
                              ! currently, the SFSEC resuspension are attributed
                              !    to the species that got scavenged,
                              !    WHICH IS NOT the species that actually
                              !    receives the resuspension
                              !    when modal_aero_wetdep_resusp_opt > 0
                              ! so when scavenged so4_c1 is resuspended as so4_a1, 
                              !    this resuspension column-tendency shows
                              !    up in so4_c1SFSES
                              ! this is done to allow better tracking of the
                              !    resuspension in the mass-budget post-processing scripts

   real(r8), intent(in)    :: rprd(pcols,pver)  ! conv precip production  rate (gathered)
   real(r8), intent(in)    :: evapc(pcols,pver)  ! conv precip evaporation rate (gathered)
   real(r8), intent(in)    :: dp_i(pver) ! pressure thickness of level (in mb)

   integer,  intent(in)    :: icol  ! normal (ungathered) i index for current column
   integer,  intent(in)    :: ktop  ! index of top cloud level for current column
   integer,  intent(in)    :: lun    ! logical unit for diagnostic output
   integer,  intent(in)    :: idiag_prevap ! flag for diagnostic output
   integer,  intent(in)    :: lchnk  ! chunk index

   logical,  intent(in)    :: doconvproc_extd(pcnst_extd)  ! indicates which species to process
   integer,  intent(in)    :: species_class(:) ! REASTER 08/05/2015
   integer,  intent(in)    :: mam_prevap_resusp_optaa ! REASTER 08/05/2015

!-----------------------------------------------------------------------
! local variables
   integer  :: k, l, ll, lspec, lspectype
   integer  :: m, m2, mmtoo, n, numtoo
   real(r8) :: d1p_prevap_resusp     ! dry-diameter of one resuspended aerosol particle (m) 
   real(r8) :: del_pr_flux_prod      ! change to precip flux from production  [(kg/kg/s)*mb]
   real(r8) :: del_pr_flux_evap      ! change to precip flux from evaporation [(kg/kg/s)*mb]
   real(r8) :: del_wd_flux_evap      ! change to wet deposition flux from evaporation [(kg/kg/s)*mb]
   real(r8) :: fdel_pr_flux_evap     ! fractional change to precip flux from evaporation
   real(r8) :: pr_flux               ! precip flux at base of current layer [(kg/kg/s)*mb]
   real(r8) :: pr_flux_old
   real(r8) :: pr_flux_tmp
   real(r8) :: pr_flux_base          ! precip flux at an effective cloud base for calculations in a particular layer
   real(r8) :: specdens_prevap_resusp(pcnst)
   real(r8) :: tmpa, tmpb, tmpc, tmpd
   real(r8) :: tmpdp                 ! delta-pressure (mb)
   real(r8) :: u_old, u_tmp
   real(r8) :: v1p_prevap_resusp     ! dry-volume of one resuspended aerosol particle (m^3)
   real(r8) :: wd_flux(pcnst_extd)   ! tracer wet deposition flux at base of current layer [(kg/kg/s)*mb]
   real(r8) :: wd_flux_tmp(pcnst_extd)
   real(r8) :: x_old, x_tmp, x_ratio
!-----------------------------------------------------------------------

!
! *** note use of non-standard units
!
! precip
!    tmpdp = dp_i is mb
!    rprd and evapc are kgwtr/kgair/s
!    pr_flux = tmpdp*rprd is mb*kgwtr/kgair/s
! this works ok because the only important thing is fdel_pr_flux_evap which is dimensionless
!
! precip-borne aerosol
!    dcondt_wetdep is kgaero/kgair/s
!    wd_flux = tmpdp*dcondt_wetdep is mb*kgaero/kgair/s
!    dcondt_prevap = del_wd_flux_evap/tmpdp is kgaero/kgair/s
! so this works ok too
!
! *** dilip switched from tmpdg to tmpdpg = tmpdp/gravit
! that is incorrect, but probably does not matter
!    for precip, the u_old and u_tmp are dimensionless
!    for aerosol, wd_flux units do not matter
!        only important thing is that tmpdp (or tmpdpg) is used
!        consistently when going from dcondt to wd_flux then to dcondt
! 

   pr_flux = 0.0_r8
   pr_flux_base = 0.0_r8
   wd_flux(:) = 0.0_r8

   if (idiag_prevap > 0) then
      write(lun,'(a,i9,i4,5x,a)') 'qakx - lchnk,i', lchnk, icol, &
         '// k; pr_flux old,new; delprod,devap; mode-1 numb wetdep,prevap; mass ...'
   end if

   do k = ktop, pver
      tmpdp = dp_i(k)

! note - setting pr_flux_tmp, wd_flux_tmp, u_old, u_tmp, x_old, x_tmp
!    to zero at this point is not necessary since they are all calculated
!    "fresh" in each iteration of the "do k" loop
! no big deal except it clutters up the code

!     pr_flux_old = pr_flux
!     del_pr_flux_prod = tmpdp*max(0.0_r8, rprd(icol,k))
!     pr_flux = pr_flux_old + del_pr_flux_prod

!     del_pr_flux_evap = min( pr_flux, tmpdp*max(0.0_r8, evapc(icol,k)) )
!     fdel_pr_flux_evap = del_pr_flux_evap / max(pr_flux, 1.0e-35_r8)

! step 1 - precip evaporation and aerosol resuspension
      tmpa = max( 0.0_r8, evapc(icol,k)*tmpdp )
      pr_flux_tmp = max( 0.0_r8, pr_flux - tmpa )
      pr_flux_tmp = min( pr_flux_base, pr_flux_tmp )

      if (pr_flux_base < 1.0e-30_r8) then
         ! when pr_flux_base=0, set u=0 to force 100% resuspension
         u_old = 1.0_r8 ; x_old = 1.0_r8
         u_tmp = 0.0_r8 ; x_tmp = 0.0_r8
         x_ratio = 0.0_r8
         pr_flux_base = 0.0_r8    ! this will start things fresh at the next layer 
         pr_flux_tmp  = 0.0_r8    ! (the next layer will then have u_old = 1)
      else
         u_old = pr_flux/pr_flux_base
         u_old = max( 0.0_r8, min( 1.0_r8, u_old ) )
         x_old = 1.0_r8 - faer_resusp_vs_fprec_evap_mpln( 1.0_r8-u_old, 2)
         x_old = max( 0.0_r8, min( 1.0_r8, x_old ) )

         u_tmp = pr_flux_tmp/pr_flux_base
         u_tmp = max( 0.0_r8, min( 1.0_r8, u_tmp ) )
         u_tmp = min( u_tmp, u_old )
         x_tmp = 1.0_r8 - faer_resusp_vs_fprec_evap_mpln( 1.0_r8-u_tmp, 2)
         x_tmp = max( 0.0_r8, min( 1.0_r8, x_tmp ) )
         x_tmp = min( x_tmp, x_old )
   
         if (x_tmp < 1.0e-30_r8) then  ! or check on x?  note that should have x_tmp >= x
            x_ratio = 0.0_r8
            pr_flux_base = 0.0_r8    ! this will start things fresh at the next layer 
            pr_flux_tmp  = 0.0_r8    ! (the next layer will then have u_old = 1)
         else
            x_ratio = x_tmp/x_old
         end if
      end if

! step 2 - precip production and aerosol scavenging
      tmpa = max( 0.0_r8, rprd(icol,k)*tmpdp )
      pr_flux_base = max( 0.0_r8, pr_flux_base + tmpa )
      pr_flux      = max( 0.0_r8, pr_flux_tmp  + tmpa )
      pr_flux      = min( pr_flux_base, pr_flux )


      do m = 2, pcnst_extd
         if ( .not. doconvproc_extd(m) ) cycle

! step 1 again, but only the aerosol resuspension
! wd_flux_tmp (updated) = (wd_flux coming into the layer) - (resuspension ! decrement)
         wd_flux_tmp(m) = max( 0.0_r8, wd_flux(m) * x_ratio )
         del_wd_flux_evap = max( 0.0_r8, wd_flux(m) - wd_flux_tmp(m) )

! step 2 again, but only the aerosol scavenging
! wd_flux (updated) = (wd_flux after resuspension) - (scavenging increment)
         tmpa = max( 0.0_r8, -dcondt_wetdep(m,k)*tmpdp )
         wd_flux(m) = max( 0.0_r8, wd_flux_tmp(m) + tmpa )

         tmpc = del_wd_flux_evap/tmpdp

         m2 = mod( m-1, pcnst ) + 1 ! for interstitial m2=m;  for activated m2=m-pcnst
         mmtoo = mmtoo_prevap_resusp(m2)

         if ( species_class(m2) == spec_class_aerosol ) then
            if (mmtoo > 0) then
               ! current species is an aerosol mass species 
               !    because mmtoo_resusp <= 0 for aerosol number species

               ! add the precip-evap (resuspension) to the history-tendency of the current species
               dcondt_prevap_hist(m,k) = dcondt_prevap_hist(m,k) + tmpc
               ! add the precip-evap (resuspension) to the actual tendencies
               ! of appropriate 
               ! coarse-mode species
               dcondt_prevap(mmtoo,k) = dcondt_prevap(mmtoo,k) + tmpc
               dcondt(mmtoo,k) = dcondt(mmtoo,k) + tmpc
            end if

         else ! ( species_class(m2) /= spec_class_aerosol ) 
            ! do this for trace gases (although currently modal_aero_convproc
            ! does not treat trace gases)
            dcondt_prevap_hist(m,k) = dcondt_prevap_hist(m,k) + tmpc
            dcondt_prevap(m,k) = dcondt_prevap(m,k) + tmpc
            dcondt(m,k) = dcondt(m,k) + tmpc
         end if

      end do ! m = 2, pcnst_extd

      if (idiag_prevap > 0) then
         n = 1
         l = numptrcw_amode(n) + pcnst
         tmpa = dcondt_wetdep(l,k)
         tmpb = dcondt_prevap(l,k)
         tmpc = 0.0_r8
         tmpd = 0.0_r8
         do ll = 1, nspec_amode(n)
            l = lmassptrcw_amode(ll,n) + pcnst
            tmpc = tmpc + dcondt_wetdep(l,k)
            tmpd = tmpd + dcondt_prevap(l,k)
         end do
         write(lun,'(a,i4,1p,4(2x,2e10.2))') 'qakx', k, &
            pr_flux_old, pr_flux, del_pr_flux_prod, -del_pr_flux_evap, &
            -tmpa, tmpb, -tmpc, tmpd
      end if

   end do ! k

   return
   end subroutine ma_precpevap30_convproc



!=========================================================================================
   subroutine ma_activate_convproc(             &
              conu,       dconudt,   conent,    &
              f_ent,      dt_u,      wup,       &
              tair,       rhoair,    fracice,   &
              pcnst_extd, lun,       idiag_act, &
              lchnk,      i,         k,         &
              ipass_calc_updraft                )
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate activation of aerosol species in convective updraft
! for a single column and level
!
! Method:
! conu(l)    = Updraft TMR (tracer mixing ratio) at k/k-1 interface
! conent(l)  = TMR of air that is entrained into the updraft from level k
! f_ent      = Fraction of the "before-detrainment" updraft massflux at 
!              k/k-1 interface" resulting from entrainment of level k air
!              (where k is the current level in subr ma_convproc_tend)
!
! On entry to this routine, the conu(l) represents the updraft TMR
! after entrainment, but before chemistry/physics and detrainment, 
! and is equal to
!    conu(l) = f_ent*conent(l) + (1.0-f_ent)*conu_below(l)
! where 
!    conu_below(l) = updraft TMR at the k+1/k interface, and
!    f_ent   = (eudp/mu_p_eudp) is the fraction of the updraft massflux 
!              from level k entrainment
!
! This routine applies aerosol activation to the entrained tracer,
! then adjusts the conu so that on exit,
!   conu(la) = conu_incoming(la) - f_ent*conent(la)*f_act(la)
!   conu(lc) = conu_incoming(lc) + f_ent*conent(la)*f_act(la)
! where 
!   la, lc   = indices for an unactivated/activated aerosol component pair
!   f_act    = fraction of conent(la) that is activated.  The f_act are
!              calculated with the Razzak-Ghan activation parameterization.
!              The f_act differ for each mode, and for number/surface/mass.
!
! Note:  At the lowest layer with cloud water, subr convproc calls this 
! routine with conent==conu and f_ent==1.0, with the result that 
! activation is applied to the entire updraft tracer flux
!
! *** The updraft velocity used for activation calculations is rather 
!     uncertain and needs more work.  However, an updraft of 1-3 m/s 
!     will activate essentially all of accumulation and coarse mode particles.
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

   use ppgrid, only: pver
   use constituents, only: pcnst, cnst_name
   use ndrop, only: activate_modal

   use modal_aero_data, only:  lmassptr_amode, lmassptrcw_amode, &
      lspectype_amode, ntot_amode, &
      nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, &
      sigmag_amode, specdens_amode, spechygro, &
      voltonumblo_amode, voltonumbhi_amode

   implicit none

!-----------------------------------------------------------------------
! arguments  (note:  TMR = tracer mixing ratio)
   integer, intent(in)     :: pcnst_extd
   ! conu = tracer mixing ratios in updraft at top of this (current) level
   !        The conu are changed by activation
   real(r8), intent(inout) :: conu(pcnst_extd)    
   ! conent = TMRs in the entrained air at this level
   real(r8), intent(in)    :: conent(pcnst_extd)
   real(r8), intent(inout) :: dconudt(pcnst_extd) ! TMR tendencies due to activation

   real(r8), intent(in)    :: f_ent  ! fraction of updraft massflux that was
                                     ! entrained across this layer == eudp/mu_p_eudp
   real(r8), intent(in)    :: dt_u   ! lagrangian transport time (s) in the 
                                     ! updraft at current level
   real(r8), intent(in)    :: wup    ! mean updraft vertical velocity (m/s)
                                     ! at current level updraft

   real(r8), intent(in)    :: tair   ! Temperature in Kelvin
   real(r8), intent(in)    :: rhoair ! air density (kg/m3)

   real(r8), intent(in)    :: fracice ! Fraction of ice within the cloud
                                     ! used as in-cloud wet removal rate
   integer,  intent(in)    :: lun    ! logical unit for diagnostic output
   integer,  intent(in)    :: idiag_act ! flag for diagnostic output
   integer,  intent(in)    :: lchnk  ! chunk index
   integer,  intent(in)    :: i      ! column index
   integer,  intent(in)    :: k      ! level index
   integer,  intent(in)    :: ipass_calc_updraft

!-----------------------------------------------------------------------
! local variables
   integer  :: l, ll, la, lc, n

   real(r8) :: delact                ! working variable
   real(r8) :: dt_u_inv              ! 1.0/dt_u
   real(r8) :: fluxm(ntot_amode)      ! to understand this, see subr activate_modal
   real(r8) :: fluxn(ntot_amode)      ! to understand this, see subr activate_modal
   real(r8) :: flux_fullact           ! to understand this, see subr activate_modal
   real(r8) :: fm(ntot_amode)         ! mass fraction of aerosols activated
   real(r8) :: fn(ntot_amode)         ! number fraction of aerosols activated
   real(r8) :: hygro(ntot_amode)      ! current hygroscopicity for int+act
   real(r8) :: naerosol(ntot_amode)   ! interstitial+activated number conc (#/m3)
   real(r8) :: sigw                  ! standard deviation of updraft velocity (cm/s)
   real(r8) :: tmpa, tmpb, tmpc      ! working variable
   real(r8) :: tmp_fact              ! working variable
   real(r8) :: vaerosol(ntot_amode)   ! int+act volume (m3/m3)
   real(r8) :: wbar                  ! mean updraft velocity (cm/s)
   real(r8) :: wdiab                 ! diabatic vertical velocity (cm/s)
   real(r8) :: wminf, wmaxf          ! limits for integration over updraft spectrum (cm/s)
   real(r8) :: lnsmax

!-----------------------------------------------------------------------


! when ipass_calc_updraft == 2, apply the activation tendencies
!    from pass 1, but multiplied by factor_reduce_actfrac
! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
   if (ipass_calc_updraft == 2) then

   dt_u_inv = 1.0_r8/dt_u
   do n = 1, ntot_amode
      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
         end if

         delact = dconudt(lc)*dt_u * factor_reduce_actfrac
         delact = min( delact, conu(la) )
         delact = max( delact, 0.0_r8 )
         conu(la) = conu(la) - delact
         conu(lc) = conu(lc) + delact
         dconudt(la) = -delact*dt_u_inv
         dconudt(lc) =  delact*dt_u_inv
      end do
   end do   ! "n = 1, ntot_amode"
   return

   end if ! (ipass_calc_updraft == 2)


! check f_ent > 0
   if (f_ent <= 0.0) return


   do n = 1, ntot_amode
! compute a (or a+cw) volume and hygroscopicity
      tmpa = 0.0_r8
      tmpb = 0.0_r8
      do ll = 1, nspec_amode(n)
         tmpc = max( conent(lmassptr_amode(ll,n)), 0.0_r8 )
         if ( use_cwaer_for_activate_maxsat ) &
         tmpc = tmpc + max( conent(lmassptrcw_amode(ll,n)+pcnst), 0.0_r8 )
         tmpc = tmpc / specdens_amode(lspectype_amode(ll,n))
         tmpa = tmpa + tmpc
         tmpb = tmpb + tmpc * spechygro(lspectype_amode(ll,n))
      end do
      vaerosol(n) = tmpa * rhoair
      if (tmpa < 1.0e-35_r8) then
         hygro(n) = 0.2_r8
      else
         hygro(n) = tmpb/tmpa
      end if

! load a (or a+cw) number and bound it
      tmpa = max( conent(numptr_amode(n)), 0.0_r8 )
      if ( use_cwaer_for_activate_maxsat ) &
      tmpa = tmpa + max( conent(numptrcw_amode(n)+pcnst), 0.0_r8 )
      naerosol(n) = tmpa * rhoair
      naerosol(n) = max( naerosol(n),   &
                         vaerosol(n)*voltonumbhi_amode(n) )
      naerosol(n) = min( naerosol(n),   &
                         vaerosol(n)*voltonumblo_amode(n) )

! diagnostic output for testing/development
!      if (lun > 0) then
!         if (n == 1) then
!            write(lun,9500)
!            write(lun,9510) (cnst_name(l), conu(l), l=1,pcnst_extd)
!            write(lun,9520) tair, rhoaircgs, airconcgs
!         end if
!         write(lun,9530) n, ntype(n), vaerosol
!         write(lun,9540) naerosol(n), tmp*airconcgs, &
!                         voltonumbhi_amode(n), voltonumblo_amode(n)
!         write(lun,9550) (maerosol(l,n), l=1,ntype(n))
!9500     format( / 'activate_conv output -- conu values' )
!9510     format( 3( a, 1pe11.3, 4x ) )
!9520     format( 'ta, rhoa, acon     ', 3(1pe11.3) )
!9530     format( 'n, ntype, sg, vol  ', i6, i5, 2(1pe11.3) )
!9540     format( 'num, num0, v2nhi&lo', 4(1pe11.3) )
!9550     format( 'masses             ', 6(1pe11.3) )
!      end if

   end do


! call Razzak-Ghan activation routine with single updraft
   wbar = max( wup, 0.5_r8 )  ! force wbar >= 0.5 m/s for now
   sigw = 0.0
   wdiab = 0.0
   wminf = wbar
   wmaxf = wbar

!  -ubroutine activate_modal(                            &
!        wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
!        na, pmode, nmode, volume, sigman, hygro,        &
!        fn, fm, fluxn, fluxm, flux_fullact              )
!     real(r8) wbar          ! grid cell mean vertical velocity (m/s)
!     real(r8) sigw          ! subgrid standard deviation of vertical vel (m/s)
!     real(r8) wdiab         ! diabatic vertical velocity (0 if adiabatic)
!     real(r8) wminf         ! minimum updraft velocity for integration (m/s)
!     real(r8) wmaxf         ! maximum updraft velocity for integration (m/s)
!     real(r8) tair          ! air temperature (K)
!     real(r8) rhoair        ! air density (kg/m3)
!     real(r8) na(pmode)     ! aerosol number concentration (/m3)
!     integer pmode          ! dimension of modes
!     integer nmode          ! number of aerosol modes
!     real(r8) volume(pmode) ! aerosol volume concentration (m3/m3)
!     real(r8) sigman(pmode) ! geometric standard deviation of aerosol size distribution
!     real(r8) hygro(pmode)  ! hygroscopicity of aerosol mode

   call activate_modal(                                                    &
         wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
         naerosol, ntot_amode, vaerosol, hygro, lnsmax,                    &
         fn, fm, fluxn, fluxm, flux_fullact                                )
   


! diagnostic output for testing/development
   if (idiag_act > 0) then
      n = min( ntot_amode, 3 )
      write(lun, '(a,i3,2f6.3, 1p,2(2x,3e10.2), 0p,3(2x,3f6.3) )' ) &
         'qaku k,w,qn,qm,hy,fn,fm', k, wup, wbar, &
         naerosol(1:n)/rhoair, vaerosol(1:n)*1.8e3/rhoair, &
         hygro(1:n), fn(1:n), fm(1:n)
         ! convert naer, vaer to number and (approx) mass TMRs
   end if
!   if (lun > 0) then
!      write(lun,9560) (fn(n), n=1,ntot_amode)
!      write(lun,9570) (fm(n), n=1,ntot_amode)
!9560  format( 'fnact values       ', 6(1pe11.3) )
!9570  format( 'fmact values       ', 6(1pe11.3) )
!   end if

      
! apply the activation fractions to the updraft aerosol mixing ratios
   dt_u_inv = 1.0_r8/dt_u

   do n = 1, ntot_amode
      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
            tmp_fact = fn(n)
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
            tmp_fact = fm(n)
         end if

         if ( (method_reduce_actfrac == 1)      .and. &
              (factor_reduce_actfrac >= 0.0_r8) .and. &
              (factor_reduce_actfrac <  1.0_r8) )     &
              tmp_fact = tmp_fact * factor_reduce_actfrac

         delact = min( conent(la)*tmp_fact*f_ent, conu(la) )
         delact = max( delact, 0.0_r8 )
         conu(la) = conu(la) - delact
         conu(lc) = conu(lc) + delact
         dconudt(la) = -delact*dt_u_inv
         dconudt(lc) =  delact*dt_u_inv
      end do
   end do   ! "n = 1, ntot_amode"

   return
   end subroutine ma_activate_convproc



!=========================================================================================
   subroutine ma_activate_convproc_method2(     &
              conu,       dconudt,              &
              f_ent,      dt_u,      wup,       &
              tair,       rhoair,    fracice,   &
              pcnst_extd, lun,       idiag_act, &
              lchnk,      i,         k,         &
              kactfirst,  ipass_calc_updraft    )
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate activation of aerosol species in convective updraft
! for a single column and level
!
! Method:
! conu(l)    = Updraft TMR (tracer mixing ratio) at k/k-1 interface
! f_ent      = Fraction of the "before-detrainment" updraft massflux at 
!              k/k-1 interface" resulting from entrainment of level k air
!              (where k is the current level in subr ma_convproc_tend)
!
! On entry to this routine, the conu(l) represents the updraft TMR
! after entrainment, but before chemistry/physics and detrainment. 
!
! This routine applies aerosol activation to the conu tracer mixing ratios,
! then adjusts the conu so that on exit,
!   conu(la) = conu_incoming(la) - conu(la)*f_act(la)
!   conu(lc) = conu_incoming(lc) + conu(la)*f_act(la)
! where 
!   la, lc   = indices for an unactivated/activated aerosol component pair
!   f_act    = fraction of conu(la) that is activated.  The f_act are
!              calculated with the Razzak-Ghan activation parameterization.
!              The f_act differ for each mode, and for number/surface/mass.
!
! At cloud base (k==kactfirst), primary activation is done using the
! "standard" code in subr activate do diagnose maximum supersaturation.
! Above cloud base, secondary activation is done using a
! prescribed supersaturation.
!
! *** The updraft velocity used for activation calculations is rather 
!     uncertain and needs more work.  However, an updraft of 1-3 m/s 
!     will activate essentially all of accumulation and coarse mode particles.
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

   use ppgrid, only: pver
   use constituents, only: pcnst, cnst_name
   use ndrop, only: activate_modal

   use modal_aero_data, only:  lmassptr_amode, lmassptrcw_amode, &
      lspectype_amode, ntot_amode, &
      nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, &
      sigmag_amode, specdens_amode, spechygro, &
      voltonumblo_amode, voltonumbhi_amode

   implicit none

!-----------------------------------------------------------------------
! arguments  (note:  TMR = tracer mixing ratio)
   integer, intent(in)     :: pcnst_extd
   ! conu = tracer mixing ratios in updraft at top of this (current) level
   !        The conu are changed by activation
   real(r8), intent(inout) :: conu(pcnst_extd)    
   real(r8), intent(inout) :: dconudt(pcnst_extd) ! TMR tendencies due to activation

   real(r8), intent(in)    :: f_ent  ! fraction of updraft massflux that was
                                     ! entrained across this layer == eudp/mu_p_eudp
   real(r8), intent(in)    :: dt_u   ! lagrangian transport time (s) in the 
                                     ! updraft at current level
   real(r8), intent(in)    :: wup    ! mean updraft vertical velocity (m/s)
                                     ! at current level updraft

   real(r8), intent(in)    :: tair   ! Temperature in Kelvin
   real(r8), intent(in)    :: rhoair ! air density (kg/m3)

   real(r8), intent(in)    :: fracice ! Fraction of ice within the cloud
                                     ! used as in-cloud wet removal rate
   integer,  intent(in)    :: lun    ! logical unit for diagnostic output
   integer,  intent(in)    :: idiag_act ! flag for diagnostic output
   integer,  intent(in)    :: lchnk  ! chunk index
   integer,  intent(in)    :: i      ! column index
   integer,  intent(in)    :: k      ! level index
   integer,  intent(in)    :: kactfirst ! k at cloud base
   integer,  intent(in)    :: ipass_calc_updraft

!-----------------------------------------------------------------------
! local variables
   integer  :: l, ll, la, lc, n

   real(r8) :: delact                ! working variable
   real(r8) :: dt_u_inv              ! 1.0/dt_u
   real(r8) :: fluxm(ntot_amode)      ! to understand this, see subr activate_modal
   real(r8) :: fluxn(ntot_amode)      ! to understand this, see subr activate_modal
   real(r8) :: flux_fullact           ! to understand this, see subr activate_modal
   real(r8) :: fm(ntot_amode)         ! mass fraction of aerosols activated
   real(r8) :: fn(ntot_amode)         ! number fraction of aerosols activated
   real(r8) :: hygro(ntot_amode)      ! current hygroscopicity for int+act
   real(r8) :: naerosol(ntot_amode)   ! interstitial+activated number conc (#/m3)
   real(r8) :: sigw                  ! standard deviation of updraft velocity (cm/s)
   real(r8) :: smax_prescribed       ! prescribed supersaturation for secondary activation (0-1 fraction)
   real(r8) :: tmpa, tmpb, tmpc      ! working variable
   real(r8) :: tmp_fact              ! working variable
   real(r8) :: vaerosol(ntot_amode)   ! int+act volume (m3/m3)
   real(r8) :: wbar                  ! mean updraft velocity (cm/s)
   real(r8) :: wdiab                 ! diabatic vertical velocity (cm/s)
   real(r8) :: wminf, wmaxf          ! limits for integration over updraft spectrum (cm/s)
   real(r8) :: lnsmax

!-----------------------------------------------------------------------


! when ipass_calc_updraft == 2, apply the activation tendencies
!    from pass 1, but multiplied by factor_reduce_actfrac
! (can only have ipass_calc_updraft == 2 when method_reduce_actfrac = 2)
   if (ipass_calc_updraft == 2) then

   dt_u_inv = 1.0_r8/dt_u
   do n = 1, ntot_amode
      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
         end if

         delact = dconudt(lc)*dt_u * factor_reduce_actfrac
         delact = min( delact, conu(la) )
         delact = max( delact, 0.0_r8 )
         conu(la) = conu(la) - delact
         conu(lc) = conu(lc) + delact
         dconudt(la) = -delact*dt_u_inv
         dconudt(lc) =  delact*dt_u_inv
      end do
   end do   ! "n = 1, ntot_amode"
   return

   end if ! (ipass_calc_updraft == 2)


! check f_ent > 0
   if (f_ent <= 0.0) return


   do n = 1, ntot_amode
! compute a (or a+cw) volume and hygroscopicity
      tmpa = 0.0_r8
      tmpb = 0.0_r8
      do ll = 1, nspec_amode(n)
         tmpc = max( conu(lmassptr_amode(ll,n)), 0.0_r8 )
         if ( use_cwaer_for_activate_maxsat ) &
         tmpc = tmpc + max( conu(lmassptrcw_amode(ll,n)+pcnst), 0.0_r8 )
         tmpc = tmpc / specdens_amode(lspectype_amode(ll,n))
         tmpa = tmpa + tmpc
         tmpb = tmpb + tmpc * spechygro(lspectype_amode(ll,n))
      end do
      vaerosol(n) = tmpa * rhoair
      if (tmpa < 1.0e-35_r8) then
         hygro(n) = 0.2_r8
      else
         hygro(n) = tmpb/tmpa
      end if

! load a (or a+cw) number and bound it
      tmpa = max( conu(numptr_amode(n)), 0.0_r8 )
      if ( use_cwaer_for_activate_maxsat ) &
      tmpa = tmpa + max( conu(numptrcw_amode(n)+pcnst), 0.0_r8 )
      naerosol(n) = tmpa * rhoair
      naerosol(n) = max( naerosol(n),   &
                         vaerosol(n)*voltonumbhi_amode(n) )
      naerosol(n) = min( naerosol(n),   &
                         vaerosol(n)*voltonumblo_amode(n) )

! diagnostic output for testing/development
!      if (lun > 0) then
!         if (n == 1) then
!            write(lun,9500)
!            write(lun,9510) (cnst_name(l), conu(l), l=1,pcnst_extd)
!            write(lun,9520) tair, rhoaircgs, airconcgs
!         end if
!         write(lun,9530) n, ntype(n), vaerosol
!         write(lun,9540) naerosol(n), tmp*airconcgs, &
!                         voltonumbhi_amode(n), voltonumblo_amode(n)
!         write(lun,9550) (maerosol(l,n), l=1,ntype(n))
!9500     format( / 'activate_conv output -- conu values' )
!9510     format( 3( a, 1pe11.3, 4x ) )
!9520     format( 'ta, rhoa, acon     ', 3(1pe11.3) )
!9530     format( 'n, ntype, sg, vol  ', i6, i5, 2(1pe11.3) )
!9540     format( 'num, num0, v2nhi&lo', 4(1pe11.3) )
!9550     format( 'masses             ', 6(1pe11.3) )
!      end if

   end do


! call Razzak-Ghan activation routine with single updraft
   wbar = max( wup, 0.5_r8 )  ! force wbar >= 0.5 m/s for now
   sigw = 0.0
   wdiab = 0.0
   wminf = wbar
   wmaxf = wbar

!  -ubroutine activate_modal(                                &
!        wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,      &
!        na, pmode, nmode, volume, sigman, hygro,            &
!        fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed )
!     real(r8) wbar          ! grid cell mean vertical velocity (m/s)
!     real(r8) sigw          ! subgrid standard deviation of vertical vel (m/s)
!     real(r8) wdiab         ! diabatic vertical velocity (0 if adiabatic)
!     real(r8) wminf         ! minimum updraft velocity for integration (m/s)
!     real(r8) wmaxf         ! maximum updraft velocity for integration (m/s)
!     real(r8) tair          ! air temperature (K)
!     real(r8) rhoair        ! air density (kg/m3)
!     real(r8) na(pmode)     ! aerosol number concentration (/m3)
!     integer pmode          ! dimension of modes
!     integer nmode          ! number of aerosol modes
!     real(r8) volume(pmode) ! aerosol volume concentration (m3/m3)
!     real(r8) sigman(pmode) ! geometric standard deviation of aerosol size distribution
!     real(r8) hygro(pmode)  ! hygroscopicity of aerosol mode
!     real(r8), optional :: smax_prescribed  ! prescribed max. supersaturation for secondary activation
   if (k == kactfirst) then
! at cloud base - do primary activation
      call activate_modal(                                                 &
         wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
         naerosol, ntot_amode, vaerosol, hygro, lnsmax,                    &
         fn, fm, fluxn, fluxm, flux_fullact                                )


   else
! above cloud base - do secondary activation with prescribed supersat 
! that is constant with height
      smax_prescribed = method2_activate_smaxmax
      call activate_modal(                                                 &
         wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,                    &
         naerosol, ntot_amode, vaerosol, hygro, lnsmax,                    &
         fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed               )
   end if


! diagnostic output for testing/development
   if (idiag_act > 0) then
      n = min( ntot_amode, 3 )
      write(lun, '(a,i3,2f6.3, 1p,2(2x,3e10.2), 0p,3(2x,3f6.3) )' ) &
         'qaku k,w,qn,qm,hy,fn,fm', k, wup, wbar, &
         naerosol(1:n)/rhoair, vaerosol(1:n)*1.8e3/rhoair, &
         hygro(1:n), fn(1:n), fm(1:n)
         ! convert naer, vaer to number and (approx) mass TMRs
   end if
!   if (lun > 0) then
!      write(lun,9560) (fn(n), n=1,ntot_amode)
!      write(lun,9570) (fm(n), n=1,ntot_amode)
!9560  format( 'fnact values       ', 6(1pe11.3) )
!9570  format( 'fmact values       ', 6(1pe11.3) )
!   end if

      
! apply the activation fractions to the updraft aerosol mixing ratios
   dt_u_inv = 1.0_r8/dt_u

   do n = 1, ntot_amode
      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
            tmp_fact = fn(n)
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
            tmp_fact = fm(n)
         end if

         if ( (method_reduce_actfrac == 1)      .and. &
              (factor_reduce_actfrac >= 0.0_r8) .and. &
              (factor_reduce_actfrac <  1.0_r8) )     &
              tmp_fact = tmp_fact * factor_reduce_actfrac

         delact = min( conu(la)*tmp_fact, conu(la) )
         delact = max( delact, 0.0_r8 )
         conu(la) = conu(la) - delact
         conu(lc) = conu(lc) + delact
         dconudt(la) = -delact*dt_u_inv
         dconudt(lc) =  delact*dt_u_inv
      end do
   end do   ! "n = 1, ntot_amode"

   return
   end subroutine ma_activate_convproc_method2



!=========================================================================================
   subroutine ma_resuspend_convproc(                           &
              dcondt,  dcondt_resusp,                          &
              const,   dp_i,          ktop,  kbot_prevap,  pcnst_extd ) ! REASTER 08/05/2015
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate resuspension of activated aerosol species resulting from both
!    detrainment from updraft and downdraft into environment
!    subsidence and lifting of environment, which may move air from
!       levels with large-scale cloud to levels with no large-scale cloud
!
! Method:
! Three possible approaches were considered:
!
! 1. Ad-hoc #1 approach.  At each level, adjust dcondt for the activated 
!    and unactivated portions of a particular aerosol species so that the 
!    ratio of dcondt (activated/unactivate) is equal to the ratio of the 
!    mixing ratios before convection.  
!    THIS WAS IMPLEMENTED IN MIRAGE2
!
! 2. Ad-hoc #2 approach.  At each level, adjust dcondt for the activated 
!    and unactivated portions of a particular aerosol species so that the 
!    change to the activated portion is minimized (zero if possible).  The
!    would minimize effects of convection on the large-scale cloud.
!    THIS IS CURRENTLY IMPLEMENTED IN CAM5 where we assume that convective
!    clouds have no impact on the stratiform-cloudborne aerosol
!
! 3. Mechanistic approach that treats the details of interactions between
!    the large-scale and convective clouds.  (Something for the future.)
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

   use ppgrid, only: pver
   use constituents, only: pcnst

   use modal_aero_data, only:  lmassptr_amode, lmassptrcw_amode, &
      nspec_amode, ntot_amode, numptr_amode, numptrcw_amode

   implicit none

!-----------------------------------------------------------------------
! arguments
! (note:  TMR = tracer mixing ratio)
   integer,  intent(in)    :: pcnst_extd
   real(r8), intent(inout) :: dcondt(pcnst_extd,pver)
                              ! overall TMR tendency from convection
   real(r8), intent(inout) :: dcondt_resusp(pcnst_extd,pver)
                              ! portion of TMR tendency due to resuspension
                              ! (actually, due to the adjustments made here)
   real(r8), intent(in)    :: const(pcnst_extd,pver)  ! TMRs before convection

   real(r8), intent(in)    :: dp_i(pver) ! pressure thickness of level (in mb)
   integer,  intent(in)    :: ktop, kbot_prevap ! indices of top and bottom cloud levels ! REASTER 08/05/2015

!-----------------------------------------------------------------------
! local variables
   integer  :: k, ll, la, lc, n
   real(r8) :: qa, qc, qac           ! working variables (mixing ratios)
   real(r8) :: qdota, qdotc, qdotac  ! working variables (MR tendencies)
!-----------------------------------------------------------------------


   do n = 1, ntot_amode

      do ll = 0, nspec_amode(n)
         if (ll == 0) then
            la = numptr_amode(n)
            lc = numptrcw_amode(n) + pcnst
         else
            la = lmassptr_amode(ll,n)
            lc = lmassptrcw_amode(ll,n) + pcnst
         end if

! apply adjustments to dcondt for pairs of unactivated (la) and 
! activated (lc) aerosol species
         if ( (la <= 0) .or. (la > pcnst_extd) ) cycle
         if ( (lc <= 0) .or. (lc > pcnst_extd) ) cycle

         do k = ktop, kbot_prevap ! REASTER 08/05/2015
            qdota = dcondt(la,k)
            qdotc = dcondt(lc,k)
            qdotac = qdota + qdotc

! mirage2 approach
!           qa = max( const(la,k), 0.0_r8 )
!           qc = max( const(lc,k), 0.0_r8 )
!           qac = qa + qc
!           if (qac <= 0.0) then
!              dcondt(la,k) = qdotac
!              dcondt(lc,k) = 0.0
!           else
!              dcondt(la,k) = qdotac*(qa/qac)
!              dcondt(lc,k) = qdotac*(qc/qac)
!           end if

! cam5 approach
            dcondt(la,k) = qdotac
            dcondt(lc,k) = 0.0

            dcondt_resusp(la,k) = (dcondt(la,k) - qdota)
            dcondt_resusp(lc,k) = (dcondt(lc,k) - qdotc)
         end do

      end do   ! "ll = -1, nspec_amode(n)"
   end do      ! "n = 1, ntot_amode"

   return
   end subroutine ma_resuspend_convproc



!=========================================================================================



end module modal_aero_convproc
