/wrfv2_fire/phys/module_mixactivate.F

Large files files are truncated, but you can click here to view the full file

!************************************************************************
! This computer software was prepared by Battelle Memorial Institute, 
! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
! LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see chem/module_mosaic_driver.F for references and terms
! of use
!************************************************************************

MODULE module_mixactivate
PRIVATE
PUBLIC prescribe_aerosol_mixactivate, mixactivate
CONTAINS


!----------------------------------------------------------------------
!----------------------------------------------------------------------
! 06-nov-2005 rce - grid_id & ktau added to arg list
! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
      subroutine prescribe_aerosol_mixactivate (                      &
		grid_id, ktau, dtstep, naer,                          &
		rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old,       &
		z, dz8w, p_at_w, t_at_w, exch_h,                      &
        qv, qc, qi, qndrop3d,                                 &
        nsource,                                              &
		ids,ide, jds,jde, kds,kde,                            &
		ims,ime, jms,jme, kms,kme,                            &
		its,ite, jts,jte, kts,kte,                            &
		f_qc, f_qi                                            )

!        USE module_configure

! wrapper to call mixactivate for mosaic description of aerosol

	implicit none

!   subr arguments
	integer, intent(in) ::                  &
		grid_id, ktau,                  &
		ids, ide, jds, jde, kds, kde,   &
		ims, ime, jms, jme, kms, kme,   &
		its, ite, jts, jte, kts, kte

	real, intent(in) :: dtstep
	real, intent(inout) :: naer ! aerosol number (/kg)

	real, intent(in),   &
		dimension( ims:ime, kms:kme, jms:jme ) :: &
		rho_phy, th_phy, pi_phy, w,  &
		z, dz8w, p_at_w, t_at_w, exch_h

	real, intent(inout),   &
		dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old

	real, intent(in),   &
		dimension( ims:ime, kms:kme, jms:jme ) :: &
		qv, qc, qi

	real, intent(inout),   &
		dimension( ims:ime, kms:kme, jms:jme ) :: &
		qndrop3d

	real, intent(out),   &
		dimension( ims:ime, kms:kme, jms:jme) :: nsource

    LOGICAL, OPTIONAL :: f_qc, f_qi

! local vars
	integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem
	parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10)
	real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity
	real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol
	real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array
	integer i,j,k,l,m,n,p
	real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk
	integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer
      	integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
      	  waterptr_aer( maxd_asize, maxd_atype ),   &
      	  numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
	  ai_phase, cw_phase
        real dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
             dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
	     sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
	     dgnum_aer(maxd_asize, maxd_atype),    & ! median diameter (cm) of number distrib of mode
	     dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
	     mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
	     dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
! terminology:  (pi/6) * (mean-volume diameter)**3 ==
!       (volume mixing ratio of section/mode)/(number mixing ratio)
	real, dimension(ims:ime,kms:kme,jms:jme) :: &
	     ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat
	integer idrydep_onoff
	real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy
	integer msectional


	  integer ptr
	  real maer

      if(naer.lt.1.)then
	     naer=1000.e6 ! #/kg default value
      endif
	  ai_phase=1
	  cw_phase=2
	  idrydep_onoff = 0
	  msectional = 0

	  t_phy(its:ite,kts:kte,jts:jte)=th_phy(its:ite,kts:kte,jts:jte)*pi_phy(its:ite,kts:kte,jts:jte)

      ntype_aer=maxd_atype
      do n=1,ntype_aer
         nsize_aer(n)=maxd_asize
	 ncomp_aer(n)=maxd_acomp
      end do
      nphase_aer=maxd_aphase

! set properties for each type and size
       do n=1,ntype_aer
       do m=1,nsize_aer(n)
          dlo_sect( m,n )=0.01e-4    ! minimum size of section (cm)
          dhi_sect( m,n )=0.5e-4    ! maximum size of section (cm)
	  sigmag_aer(m,n)=2.      ! geometric standard deviation of aerosol size dist
	  dgnum_aer(m,n)=0.1e-4       ! median diameter (cm) of number distrib of mode
	  dpvolmean_aer(m,n) = dgnum_aer(m,n) * exp( 1.5 * (log(sigmag_aer(m,n)))**2 )
	  end do
	  do l=1,ncomp_aer(n)
	     dens_aer( l, n)=1.0   ! density (g/cm3) of material
	     mw_aer( l, n)=132. ! molecular weight (g/mole)
	  end do
      end do
       ptr=0
       do p=1,nphase_aer
       do n=1,ntype_aer
       do m=1,nsize_aer(n)
          ptr=ptr+1
          numptr_aer( m, n, p )=ptr
	  if(p.eq.ai_phase)then
	     chem(its:ite,kts:kte,jts:jte,ptr)=naer
	  else
	     chem(its:ite,kts:kte,jts:jte,ptr)=0.
	  endif
	end do ! size
	end do ! type
	end do ! phase
       do p=1,maxd_aphase
       do n=1,ntype_aer
       do m=1,nsize_aer(n)
	  do l=1,ncomp_aer(n)
          ptr=ptr+1
	     if(ptr.gt.max_chem)then
	        write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
	        call wrf_error_fatal("1")
	     endif
	     massptr_aer(l, m, n, p)=ptr
! maer is ug/kg-air;  naer is #/kg-air;  dgnum is cm;  dens_aer is g/cm3
! 1.e6 factor converts g to ug
	     maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) *   &
                 (dgnum_aer(m,n)**3) * exp( 4.5*((log(sigmag_aer(m,n)))**2) ) )
	     if(p.eq.ai_phase)then
	        chem(its:ite,kts:kte,jts:jte,ptr)=maer
	     else
	        chem(its:ite,kts:kte,jts:jte,ptr)=0.
	     endif
	  end do
	end do ! size
	end do ! type
	end do ! phase
       do n=1,ntype_aer
       do m=1,nsize_aer(n)
          ptr=ptr+1
	  if(ptr.gt.max_chem)then
	     write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
	     call wrf_error_fatal("1")
	  endif
!wig	  waterptr_aer(m, n)=ptr
	  waterptr_aer(m, n)=-1
	end do ! size
	end do ! type
	ddvel(its:ite,jts:jte,:)=0.
    hygro(its:ite,kts:kte,jts:jte,:,:) = 0.5

! 06-nov-2005 rce - grid_id & ktau added to arg list
      call mixactivate(  msectional,     &
            chem,max_chem,qv,qc,qi,qndrop3d,        &
            t_phy, w, ddvel, idrydep_onoff,  &
            maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
            ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
            numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
            dens_aer, mw_aer,           &
            waterptr_aer, hygro,  ai_phase, cw_phase,                &
            ids,ide, jds,jde, kds,kde,                            &
            ims,ime, jms,jme, kms,kme,                            &
            its,ite, jts,jte, kts,kte,                            &
            rho_phy, z, dz8w, p_at_w, t_at_w, exch_h,      &
            cldfra, cldfra_old, qsrflx,         &
            ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
            grid_id, ktau, dtstep, &
            F_QC=f_qc, F_QI=f_qi                              )


      end subroutine prescribe_aerosol_mixactivate

!----------------------------------------------------------------------
!----------------------------------------------------------------------
!   nov-04 sg ! replaced amode with aer and expanded aerosol dimension to include type and phase

! 06-nov-2005 rce - grid_id & ktau added to arg list
! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
subroutine mixactivate(  msectional,            &
           chem, num_chem, qv, qc, qi, qndrop3d,         &
           temp, w, ddvel, idrydep_onoff,  &
           maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
           ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
           numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
           dens_aer, mw_aer,               &
           waterptr_aer, hygro, ai_phase, cw_phase,              &
           ids,ide, jds,jde, kds,kde,                            &
           ims,ime, jms,jme, kms,kme,                            &
           its,ite, jts,jte, kts,kte,                            &
           rho, zm, dz8w, p_at_w, t_at_w, kvh,      &
           cldfra, cldfra_old, qsrflx,          &
           ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
           grid_id, ktau, dtstep, &
           f_qc, f_qi                       )


!     vertical diffusion and nucleation of cloud droplets
!     assume cloud presence controlled by cloud fraction
!     doesn't distinguish between warm, cold clouds

  USE module_model_constants, only: g, rhowater, xlv, cp, rvovrd, r_d, r_v, mwdry, ep_2
  USE module_radiation_driver, only: cal_cldfra

  implicit none

!     input

  INTEGER, intent(in) ::         grid_id, ktau
  INTEGER, intent(in) ::         num_chem
  integer, intent(in) ::         ids,ide, jds,jde, kds,kde,    &
                                 ims,ime, jms,jme, kms,kme,    &
                                 its,ite, jts,jte, kts,kte

  integer maxd_aphase, nphase_aer, maxd_atype, ntype_aer
  integer maxd_asize, maxd_acomp, nsize_aer(maxd_atype)
  integer, intent(in) ::   &
       ncomp_aer( maxd_atype  ),   &
       massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
       waterptr_aer( maxd_asize, maxd_atype ),   &
       numptr_aer( maxd_asize, maxd_atype, maxd_aphase), &
       ai_phase, cw_phase
  integer, intent(in) :: msectional ! 1 for sectional, 0 for modal
  integer, intent(in) :: idrydep_onoff
  real, intent(in)  ::                       &
       dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
       dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
       sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
       dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
       mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
       dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
! terminology:  (pi/6) * (mean-volume diameter)**3 ==
!       (volume mixing ratio of section/mode)/(number mixing ratio)


  REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) :: &
       chem ! aerosol molar mixing ratio (ug/kg or #/kg)

  REAL, intent(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
       qv, qc, qi ! water species (vapor, cloud drops, cloud ice) mixing ratio (g/g)

  LOGICAL, OPTIONAL :: f_qc, f_qi

  REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
       qndrop3d    ! water species mixing ratio (g/g)

  real, intent(in) :: dtstep             ! time step for microphysics (s)
  real, intent(in) :: temp(ims:ime, kms:kme, jms:jme)    ! temperature (K)
  real, intent(in) :: w(ims:ime, kms:kme, jms:jme)   ! vertical velocity (m/s)
  real, intent(in) :: rho(ims:ime, kms:kme, jms:jme)    ! density at mid-level  (kg/m3)
  REAL, intent(in) :: ddvel( its:ite, jts:jte, num_chem ) ! deposition velocity  (m/s)
  real, intent(in) :: zm(ims:ime, kms:kme, jms:jme)     ! geopotential height of level (m)
  real, intent(in) :: dz8w(ims:ime, kms:kme, jms:jme) ! layer thickness (m)
  real, intent(in) :: p_at_w(ims:ime, kms:kme, jms:jme) ! pressure at layer interface (Pa)
  real, intent(in) :: t_at_w(ims:ime, kms:kme, jms:jme) ! temperature at layer interface (K)
  real, intent(in) :: kvh(ims:ime, kms:kme, jms:jme)    ! vertical diffusivity (m2/s)
  real, intent(inout) :: cldfra_old(ims:ime, kms:kme, jms:jme)! cloud fraction on previous time step
  real, intent(inout) :: cldfra(ims:ime, kms:kme, jms:jme)    ! cloud fraction
  real, intent(in) :: hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk hygroscopicity   &

  REAL, intent(out), DIMENSION( ims:ime, jms:jme, num_chem ) ::   qsrflx ! dry deposition rate for aerosol
  real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, &  ! droplet number source (#/kg/s)
       ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat


!--------------------Local storage-------------------------------------
!
  real :: dgnum_aer(maxd_asize, maxd_atype) ! median diameter (cm) of number distrib of mode
  real :: qndrop(kms:kme)      ! cloud droplet number mixing ratio (#/kg)
  real :: lcldfra(kms:kme)     ! liquid cloud fraction
  real :: lcldfra_old(kms:kme) ! liquid cloud fraction for previous timestep
  real :: wtke(kms:kme)        ! turbulent vertical velocity at base of layer k (m2/s)
  real zn(kms:kme)             ! g/pdel (m2/g) for layer
  real zs(kms:kme)             ! inverse of distance between levels (m)
  real, parameter :: zkmin = 0.01
  real, parameter :: zkmax = 100.
  real cs(kms:kme)             ! air density (kg/m3) at layer center
  real csbot(kms:kme)          ! air density (kg/m3) at layer bottom
  real csbot_cscen(kms:kme)    ! csbot(k)/cs(k)
  real dz(kms:kme)             ! geometric thickness of layers (m)

  real wdiab                   ! diabatic vertical velocity
!      real, parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
  real, parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
!      real, parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
  real :: qndrop_new(kms:kme)  ! droplet number nucleated on cloud boundaries
  real :: ekd(kms:kme)         ! diffusivity for droplets (m2/s)
  real :: ekk(kms:kme)         ! density*diffusivity for droplets (kg/m3 m2/s)
  real :: srcn(kms:kme)        ! droplet source rate (/s)
  real, parameter :: sq2pi = 2.5066282746
  real dtinv

  integer km1,kp1
  real wbar,wmix,wmin,wmax
  real dum
  real tmpa, tmpb, tmpc, tmpc1, tmpc2, tmpd, tmpe, tmpf
  real tmpcourno
  real dact
  real fluxntot         ! (#/cm2/s)
  real fac_srflx
  real depvel_drop, depvel_tmp
  real, parameter :: depvel_uplimit = 1.0 ! upper limit for dep vels (m/s)
  real :: surfrate(num_chem) ! surface exchange rate (/s)
  real surfratemax      ! max surfrate for all species treated here
  real surfrate_drop    ! surfade exchange rate for droplelts
  real dtmin,tinv,dtt
  integer nsubmix,nsubmix_bnd
  integer i,j,k,m,n,nsub
  real dtmix
  real alogarg
  real qcld
  real pi
  integer nnew,nsav,ntemp
  real :: overlapp(kms:kme),overlapm(kms:kme) ! cloud overlap
  real ::  ekkp(kms:kme),ekkm(kms:kme) ! zn*zs*density*diffusivity
!  integer, save :: count_submix(100)=0 ! wig: Note that this is a no-no for tile threads with OMP

  integer lnum,lnumcw,l,lmass,lmasscw,lsfc,lsfccw,ltype,lsig,lwater
  integer :: ntype(maxd_asize)

  real ::  naerosol(maxd_asize, maxd_atype)    ! interstitial aerosol number conc (/m3)
  real ::  naerosolcw(maxd_asize, maxd_atype)  ! activated number conc (/m3)
  real ::   maerosol(maxd_acomp,maxd_asize, maxd_atype)   ! interstit mass conc (kg/m3)
  real ::   maerosolcw(maxd_acomp,maxd_asize, maxd_atype) ! activated mass conc (kg/m3)
  real ::   maerosol_tot(maxd_asize, maxd_atype)     ! species-total interstit mass conc (kg/m3)
  real ::   maerosol_totcw(maxd_asize, maxd_atype)   ! species-total activated mass conc (kg/m3)
  real ::   vaerosol(maxd_asize, maxd_atype) ! interstit+activated aerosol volume conc (m3/m3)
  real ::   vaerosolcw(maxd_asize, maxd_atype) ! activated aerosol volume conc (m3/m3)
  real ::   raercol(kms:kme,num_chem,2) ! aerosol mass, number mixing ratios
  real ::   source(kms:kme) !

  real ::   fn(maxd_asize, maxd_atype)         ! activation fraction for aerosol number
  real ::   fs(maxd_asize, maxd_atype)         ! activation fraction for aerosol sfcarea
  real ::   fm(maxd_asize, maxd_atype)         ! activation fraction for aerosol mass
  integer ::   ncomp(maxd_atype)

  real ::   fluxn(maxd_asize, maxd_atype)      ! number  activation fraction flux (m/s)
  real ::   fluxs(maxd_asize, maxd_atype)      ! sfcarea activation fraction flux (m/s)
  real ::   fluxm(maxd_asize, maxd_atype)      ! mass    activation fraction flux (m/s)
  real ::   flux_fullact(kms:kme)              ! 100%    activation fraction flux (m/s)
!     note:  activation fraction fluxes are defined as
!     fluxn = [flux of activated aero. number into cloud (#/m2/s)]
!           / [aero. number conc. in updraft, just below cloudbase (#/m3)]

  real :: nact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. number  activation rate (/s)
  real :: mact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. mass    activation rate (/s)
  real :: npv(maxd_asize, maxd_atype) ! number per volume concentration (/m3)
  real scale

  real :: hygro_aer(maxd_asize, maxd_atype)  ! hygroscopicity of aerosol mode
  real :: exp45logsig     ! exp(4.5*alogsig**2)
  real :: alogsig(maxd_asize, maxd_atype) ! natl log of geometric standard dev of aerosol
  integer, parameter :: psat=6  ! number of supersaturations to calc ccn concentration
  real ccn(kts:kte,psat)        ! number conc of aerosols activated at supersat
  real, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
       (/0.02,0.05,0.1,0.2,0.5,1.0/)
  real super(psat) ! supersaturation
  real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m)
  real :: ccnfact(psat,maxd_asize, maxd_atype)
  real :: amcube(maxd_asize, maxd_atype) ! cube of dry mode radius (m)
  real :: argfactor(maxd_asize, maxd_atype)
  real aten ! surface tension parameter
  real t0 ! reference temperature
  real sm ! critical supersaturation
  real arg

  integer,parameter :: icheck_colmass = 0
           ! icheck_colmass > 0 turns on mass/number conservation checking
           ! values of 1, 10, 100 produce less to more diagnostics
  integer :: colmass_worst_ij( 2, 0:maxd_acomp, maxd_asize, maxd_atype )
  integer :: colmass_maxworst_i(3)
  real :: colmass_bgn( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
  real :: colmass_end( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
  real :: colmass_sfc( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
  real :: colmass_worst( 0:maxd_acomp, maxd_asize, maxd_atype )
  real :: colmass_maxworst_r
  real :: rhodz( kts:kte ), rhodzsum

!!$#if (defined AIX)
!!$#define ERF erf
!!$#define ERFC erfc
!!$#else
!!$#define ERF erf
!!$    real erf
!!$#define ERFC erfc
!!$    real erfc
!!$#endif

  character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)


  colmass_worst(:,:,:) = 0.0
  colmass_worst_ij(:,:,:,:) = -1


  arg = 1.0
  if (abs(0.8427-ERF_ALT(arg))/0.8427>0.001) then
     write (6,*) 'erf_alt(1.0) = ',ERF_ALT(arg)
     call wrf_error_fatal('dropmixnuc: Error function error')
  endif
  arg = 0.0
  if (ERF_ALT(arg) /= 0.0) then
     write (6,*) 'erf_alt(0.0) = ',ERF_ALT(arg)
     call wrf_error_fatal('dropmixnuc: Error function error')
  endif

  pi = 4.*atan(1.0)
  dtinv=1./dtstep

  depvel_drop =  0.1 ! prescribed here rather than getting it from dry_dep_driver
  if (idrydep_onoff .le. 0) depvel_drop =  0.0
  depvel_drop =  min(depvel_drop,depvel_uplimit)

  do n=1,ntype_aer
     do m=1,nsize_aer(n)
        ncomp(n)=ncomp_aer(n)
        alogsig(m,n)=alog(sigmag_aer(m,n))
        dgnum_aer(m,n) = dpvolmean_aer(m,n) * exp( -1.5*alogsig(m,n)*alogsig(m,n) )
!    print *,'sigmag_aer,dgnum_aer=',sigmag_aer(m,n),dgnum_aer(m,n)
        ! npv is used only if number is diagnosed from volume
        npv(m,n)=6./(pi*(0.01*dgnum_aer(m,n))**3*exp(4.5*alogsig(m,n)*alogsig(m,n)))
     end do
  end do
  t0=273.15   !wig, 1-Mar-2009: Added .15
  aten=2.*surften/(r_v*t0*rhowater)
  super(:)=0.01*supersat(:)
  do n=1,ntype_aer
     do m=1,nsize_aer(n)
        exp45logsig=exp(4.5*alogsig(m,n)*alogsig(m,n))
        argfactor(m,n)=2./(3.*sqrt(2.)*alogsig(m,n))
        amcube(m,n)=3./(4.*pi*exp45logsig*npv(m,n))
     enddo
  enddo

  IF( PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
     CALL cal_cldfra(CLDFRA,qc,qi,f_qc,f_qi,      &
          ids,ide, jds,jde, kds,kde,              &
          ims,ime, jms,jme, kms,kme,              &
          its,ite, jts,jte, kts,kte               )
  END IF

  qsrflx(its:ite,jts:jte,:) = 0.

!     start loop over columns

OVERALL_MAIN_J_LOOP: do j=jts,jte
OVERALL_MAIN_I_LOOP: do i=its,ite

!        load number nucleated into qndrop on cloud boundaries

! initialization for current i .........................................

     do k=kts+1,kte
	    zs(k)=1./(zm(i,k,j)-zm(i,k-1,j))
	 enddo
	 zs(kts)=zs(kts+1)
     zs(kte+1)=0.

     do k=kts,kte
!!$	    if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
!!$!	       call wrf_error_fatal("1")
!!$	    endif
        if(f_qi)then
           qcld=qc(i,k,j)+qi(i,k,j)
        else
           qcld=qc(i,k,j)
        endif
        if(qcld.lt.-1..or.qcld.gt.1.)then
           write(6,'(a,g12.2,a,3i5)')'qcld=',qcld,' for i,k,j=',i,k,j
           call wrf_error_fatal("1")
        endif
        if(qcld.gt.1.e-20)then
           lcldfra(k)=cldfra(i,k,j)*qc(i,k,j)/qcld
           lcldfra_old(k)=cldfra_old(i,k,j)*qc(i,k,j)/qcld
        else
           lcldfra(k)=0.
           lcldfra_old(k)=0.
        endif
        qndrop(k)=qndrop3d(i,k,j)
!	    qndrop(k)=1.e5
        cs(k)=rho(i,k,j) ! air density (kg/m3)
        dz(k)=dz8w(i,k,j)
        do n=1,ntype_aer
           do m=1,nsize_aer(n)
              nact(k,m,n)=0.
              mact(k,m,n)=0.
           enddo
        enddo
        zn(k)=1./(cs(k)*dz(k))
        if(k>kts)then
           ekd(k)=kvh(i,k,j)
           ekd(k)=max(ekd(k),zkmin)
           ekd(k)=min(ekd(k),zkmax)
        else
           ekd(k)=0
        endif
!           diagnose subgrid vertical velocity from diffusivity
        if(k.eq.kts)then
           wtke(k)=sq2pi*depvel_drop
!               wtke(k)=sq2pi*kvh(i,k,j)
!               wtke(k)=max(wtke(k),wmixmin)
        else
           wtke(k)=sq2pi*ekd(k)/dz(k)
        endif
        wtke(k)=max(wtke(k),wmixmin)
        nsource(i,k,j)=0.
     enddo
     nsource(i,kte+1,j) = 0.
     qndrop(kte+1)      = 0.
     zn(kte+1)          = 0.

     do k = kts+1, kte
        tmpa = dz(k-1) ; tmpb = dz(k)
        tmpc = tmpa/(tmpa + tmpb)
        csbot(k) = cs(k-1)*(1.0-tmpc) + cs(k)*tmpc
        csbot_cscen(k) = csbot(k)/cs(k)
     end do
     csbot(kts) = cs(kts)
     csbot_cscen(kts) = 1.0
     csbot(kte+1) = cs(kte)
     csbot_cscen(kte+1) = 1.0

    !  calculate surface rate and mass mixing ratio for aerosol

     surfratemax = 0.0
     nsav=1
     nnew=2
     surfrate_drop=depvel_drop/dz(kts)
     surfratemax = max( surfratemax, surfrate_drop )
     do n=1,ntype_aer
        do m=1,nsize_aer(n)
           lnum=numptr_aer(m,n,ai_phase)
           lnumcw=numptr_aer(m,n,cw_phase)
           if(lnum>0)then
              depvel_tmp = max( 0.0, min( ddvel(i,j,lnum), depvel_uplimit ) )
              surfrate(lnum)=depvel_tmp/dz(kts)
              surfrate(lnumcw)=surfrate_drop
              surfratemax = max( surfratemax, surfrate(lnum) )
!             scale = 1000./mwdry ! moles/kg
              scale = 1.
              raercol(kts:kte,lnumcw,nsav)=chem(i,kts:kte,j,lnumcw)*scale ! #/kg
              raercol(kts:kte,lnum,nsav)=chem(i,kts:kte,j,lnum)*scale
           endif
           do l=1,ncomp(n)
              lmass=massptr_aer(l,m,n,ai_phase)
              lmasscw=massptr_aer(l,m,n,cw_phase)
!             scale = mw_aer(l,n)/mwdry
              scale = 1.e-9 ! kg/ug
              depvel_tmp = max( 0.0, min( ddvel(i,j,lmass), depvel_uplimit ) )
              surfrate(lmass)=depvel_tmp/dz(kts)
              surfrate(lmasscw)=surfrate_drop
              surfratemax = max( surfratemax, surfrate(lmass) )
              raercol(kts:kte,lmasscw,nsav)=chem(i,kts:kte,j,lmasscw)*scale ! kg/kg
              raercol(kts:kte,lmass,nsav)=chem(i,kts:kte,j,lmass)*scale ! kg/kg
           enddo
           lwater=waterptr_aer(m,n)
           if(lwater>0)then
              depvel_tmp = max( 0.0, min( ddvel(i,j,lwater), depvel_uplimit ) )
              surfrate(lwater)=depvel_tmp/dz(kts)
              surfratemax = max( surfratemax, surfrate(lwater) )
              raercol(kts:kte,lwater,nsav)=chem(i,kts:kte,j,lwater) ! don't bother to convert units,
             ! because it doesn't contribute to aerosol mass
           endif
        enddo ! size
     enddo ! type


! mass conservation checking
     if (icheck_colmass > 0) then

! calc initial column burdens
        colmass_bgn(:,:,:,:) = 0.0
        colmass_end(:,:,:,:) = 0.0
        colmass_sfc(:,:,:,:) = 0.0
        rhodz(kts:kte) = 1.0/zn(kts:kte)
        rhodzsum = sum( rhodz(kts:kte) )
        do n=1,ntype_aer
           do m=1,nsize_aer(n)
              lnum=numptr_aer(m,n,ai_phase)
              lnumcw=numptr_aer(m,n,cw_phase)
              if(lnum>0)then
                 colmass_bgn(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
                 colmass_bgn(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
              endif
              do l=1,ncomp(n)
                 lmass=massptr_aer(l,m,n,ai_phase)
                 lmasscw=massptr_aer(l,m,n,cw_phase)
                 colmass_bgn(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
                 colmass_bgn(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
              enddo
           enddo ! size
        enddo ! type
     endif ! (icheck_colmass > 0)


!        droplet nucleation/aerosol activation

! k-loop for growing/shrinking cloud calcs .............................
GROW_SHRINK_MAIN_K_LOOP: do k=kts,kte
        km1=max0(k-1,1)
        kp1=min0(k+1,kde-1)


!       if(lcldfra(k)-lcldfra_old(k).gt.0.01)then   ! this line is the "old" criterion
!       go to 10

!       growing cloud PLUS
!       upwards vertical advection when lcldfra(k-1) < lcldfra(k)
!
! tmpc1 = cloud fraction increase from previous time step
        tmpc1 = max( (lcldfra(k)-lcldfra_old(k)), 0.0 )
        if (k > kts) then
! tmpc2 = fraction of layer for which vertical advection from below
!             (over dtstep) displaces cloudy air with clear air
!       = (courant number using upwards w at layer bottom)*(difference in cloud fraction)
           tmpcourno = dtstep*max(w(i,k,j),0.0)/dz(k)
           tmpc2 = max( (lcldfra(k)-lcldfra(km1)), 0.0 ) * tmpcourno
           tmpc2 = min( tmpc2, 1.0 )
!          tmpc2 = 0.0   ! this turns off the vertical advect part
        else
           tmpc2 = 0.0
        endif

        if ((tmpc1 > 0.001) .or. (tmpc2 > 0.001)) then

!                wmix=wtke(k)
           wbar=w(i,k,j)+wtke(k)
           wmix=0.
           wmin=0.
! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
           wmax=50.
           wdiab=0

!                load aerosol properties, assuming external mixtures
           do n=1,ntype_aer
              do m=1,nsize_aer(n)
                 call loadaer(raercol(1,1,nsav),k,kms,kme,num_chem,    &
                      cs(k), npv(m,n), dlo_sect(m,n),dhi_sect(m,n),             &
                      maxd_acomp, ncomp(n), &
                      grid_id, ktau, i, j, m, n,   &
                      numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
                      dens_aer(1,n),    &
                      massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
                      maerosol(1,m,n), maerosolcw(1,m,n),          &
                      maerosol_tot(m,n), maerosol_totcw(m,n),      &
                      naerosol(m,n), naerosolcw(m,n),                  &
                      vaerosol(m,n), vaerosolcw(m,n) )

                 hygro_aer(m,n)=hygro(i,k,j,m,n)
              enddo
           enddo

! 06-nov-2005 rce - grid_id & ktau added to arg list
           call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
                msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
                naerosol, vaerosol,  &
                dlo_sect,dhi_sect,sigmag_aer,hygro_aer,              &
                fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )

           do n = 1,ntype_aer
              do m = 1,nsize_aer(n)
                 lnum   = numptr_aer(m,n,ai_phase)
                 lnumcw = numptr_aer(m,n,cw_phase)
                 if (tmpc1 > 0.0) then
                    dact = tmpc1*fn(m,n)*raercol(k,lnum,nsav) ! interstitial only
                 else
                    dact = 0.0
                 endif
                 if (tmpc2 > 0.0) then
                    dact = dact + tmpc2*fn(m,n)*raercol(km1,lnum,nsav) ! interstitial only
                 endif
                 dact = min( dact, 0.99*raercol(k,lnum,nsav) )
                 raercol(k,lnumcw,nsav) = raercol(k,lnumcw,nsav)+dact
                 raercol(k,lnum,  nsav) = raercol(k,lnum,  nsav)-dact
                 qndrop(k) = qndrop(k)+dact
                 nsource(i,k,j) = nsource(i,k,j)+dact*dtinv
                 do l = 1,ncomp(n)
                    lmass   = massptr_aer(l,m,n,ai_phase)
                    lmasscw = massptr_aer(l,m,n,cw_phase)
                    if (tmpc1 > 0.0) then
                       dact = tmpc1*fm(m,n)*raercol(k,lmass,nsav) ! interstitial only
                    else
                       dact = 0.0
                    endif
                    if (tmpc2 > 0.0) then
                       dact = dact + tmpc2*fm(m,n)*raercol(km1,lmass,nsav) ! interstitial only
                    endif
                    dact = min( dact, 0.99*raercol(k,lmass,nsav) )
                    raercol(k,lmasscw,nsav)  =  raercol(k,lmasscw,nsav)+dact
                    raercol(k,lmass,  nsav)  =  raercol(k,lmass,  nsav)-dact
                 enddo
              enddo
           enddo
!   10 continue
        endif   ! ((tmpc1 > 0.001) .or. (tmpc2 > 0.001))


        if(lcldfra(k) < lcldfra_old(k) .and. lcldfra_old(k) > 1.e-20)then   ! this line is the "old" criterion
!         go to 20

!       shrinking cloud ......................................................

!                droplet loss in decaying cloud
           nsource(i,k,j)=nsource(i,k,j)+qndrop(k)*(lcldfra(k)-lcldfra_old(k))*dtinv
           qndrop(k)=qndrop(k)*(1.+lcldfra(k)-lcldfra_old(k))
!                 convert activated aerosol to interstitial in decaying cloud

           tmpc = (lcldfra(k)-lcldfra_old(k))/lcldfra_old(k)
           do n=1,ntype_aer
              do m=1,nsize_aer(n)
                 lnum=numptr_aer(m,n,ai_phase)
                 lnumcw=numptr_aer(m,n,cw_phase)
                 if(lnum.gt.0)then
                    dact=raercol(k,lnumcw,nsav)*tmpc
                    raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
                    raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
                 endif
                 do l=1,ncomp(n)
                    lmass=massptr_aer(l,m,n,ai_phase)
                    lmasscw=massptr_aer(l,m,n,cw_phase)
                    dact=raercol(k,lmasscw,nsav)*tmpc
                    raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
                    raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
                 enddo
              enddo
           enddo
!             20 continue
        endif

     enddo GROW_SHRINK_MAIN_K_LOOP
! end of k-loop for growing/shrinking cloud calcs ......................



! ......................................................................
! start of main k-loop for calc of old cloud activation tendencies ..........
! this loop does "set up" for the nsubmix loop
!
! rce-comment
!    changed this part of code to use current cloud fraction (lcldfra) exclusively

OLD_CLOUD_MAIN_K_LOOP: do k=kts,kte
        km1=max0(k-1,kts)
        kp1=min0(k+1,kde-1)
        flux_fullact(k) = 0.0
        if(lcldfra(k).gt.0.01)then
!          go to 30

!               old cloud
              if(lcldfra(k)-lcldfra(km1).gt.0.01.or.k.eq.kts)then

!                   interior cloud
!                   cloud base

                 wdiab=0
                 wmix=wtke(k) ! spectrum of updrafts
                 wbar=w(i,k,j) ! spectrum of updrafts
!                    wmix=0. ! single updraft
!               wbar=wtke(k) ! single updraft
! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
                 wmax=50.
                 ekd(k)=wtke(k)*dz(k)/sq2pi
                 alogarg=max(1.e-20,1/lcldfra(k)-1.)
                 wmin=wbar+wmix*0.25*sq2pi*alog(alogarg)

                 do n=1,ntype_aer
                    do m=1,nsize_aer(n)
                       call loadaer(raercol(1,1,nsav),km1,kms,kme,num_chem,    &
                            cs(k), npv(m,n),dlo_sect(m,n),dhi_sect(m,n),               &
                            maxd_acomp, ncomp(n), &
                            grid_id, ktau, i, j, m, n,   &
                            numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
                            dens_aer(1,n),   &
                            massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
                            maerosol(1,m,n), maerosolcw(1,m,n),          &
                            maerosol_tot(m,n), maerosol_totcw(m,n),      &
                            naerosol(m,n), naerosolcw(m,n),                  &
                            vaerosol(m,n), vaerosolcw(m,n) )

                       hygro_aer(m,n)=hygro(i,k,j,m,n)

                    enddo
                 enddo
!          print *,'old cloud wbar,wmix=',wbar,wmix

                 call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
                      msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
                      naerosol, vaerosol,  &
                      dlo_sect,dhi_sect, sigmag_aer,hygro_aer,                    &
                      fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )
                 
! rce-comment
!    the activation-fraction fluxes (fluxn, fluxm) from subr activate assume that
!       wbar << wmix, which is valid for global-model scale but not mesoscale
!    for wrf-chem application, divide these by flux_fullact to get a 
!       "flux-weighted-average" activation fraction, then multiply by (ekd(k)*zs(k)) 
!       which is the local "turbulent vertical-mixing velocity"
                 if (k > kts) then
                    if (flux_fullact(k) > 1.0e-20) then
                       tmpa = ekd(k)*zs(k)
                       tmpf = flux_fullact(k)
                       do n=1,ntype_aer
                       do m=1,nsize_aer(n)
                          tmpb = max( fluxn(m,n), 0.0 ) / max( fluxn(m,n), tmpf )
                          fluxn(m,n) = tmpa*tmpb
                          tmpb = max( fluxm(m,n), 0.0 ) / max( fluxm(m,n), tmpf )
                          fluxm(m,n) = tmpa*tmpb
                       enddo
                       enddo
                    else
                       fluxn(:,:) = 0.0
                       fluxm(:,:) = 0.0
                    endif
                 endif

                 if(k.gt.kts)then
                    tmpc = lcldfra(k)-lcldfra(km1)
                 else
                    tmpc=lcldfra(k)
                 endif
! rce-comment
!    flux of activated mass into layer k (in kg/m2/s)
!       = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
!    source of activated mass (in kg/kg/s) = flux divergence
!       = actmassflux/(cs(i,k)*dz(i,k))
!    so need factor of csbot_cscen = csbot(k)/cs(i,k)
!                tmpe=1./(dz(k))
                 tmpe = csbot_cscen(k)/(dz(k))
                 fluxntot=0.
                 do n=1,ntype_aer
                 do m=1,nsize_aer(n)
                    fluxn(m,n)=fluxn(m,n)*tmpc
!                   fluxs(m,n)=fluxs(m,n)*tmpc
                    fluxm(m,n)=fluxm(m,n)*tmpc
                    lnum=numptr_aer(m,n,ai_phase)
                    fluxntot=fluxntot+fluxn(m,n)*raercol(km1,lnum,nsav)
!             print *,'fn=',fn(m,n),' for m,n=',m,n
!             print *,'old cloud tmpc=',tmpc,' fn=',fn(m,n),' for m,n=',m,n
                    nact(k,m,n)=nact(k,m,n)+fluxn(m,n)*tmpe
                    mact(k,m,n)=mact(k,m,n)+fluxm(m,n)*tmpe
                 enddo
                 enddo
                 flux_fullact(k) = flux_fullact(k)*tmpe
                 nsource(i,k,j)=nsource(i,k,j)+fluxntot*zs(k)
                 fluxntot=fluxntot*cs(k)
              endif
!       30 continue

        else
!       go to 40
!              no cloud
           if(qndrop(k).gt.10000.e6)then
              print *,'i,k,j,lcldfra,qndrop=',i,k,j,lcldfra(k),qndrop(k)
              print *,'cldfra,ql,qi',cldfra(i,k,j),qc(i,k,j),qi(i,k,j)
           endif
           nsource(i,k,j)=nsource(i,k,j)-qndrop(k)*dtinv
           qndrop(k)=0.
!              convert activated aerosol to interstitial in decaying cloud
           do n=1,ntype_aer
              do m=1,nsize_aer(n)
                 lnum=numptr_aer(m,n,ai_phase)
                 lnumcw=numptr_aer(m,n,cw_phase)
                 if(lnum.gt.0)then
                    raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol(k,lnumcw,nsav)
                    raercol(k,lnumcw,nsav)=0.
                 endif
                 do l=1,ncomp(n)
                    lmass=massptr_aer(l,m,n,ai_phase)
                    lmasscw=massptr_aer(l,m,n,cw_phase)
                    raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol(k,lmasscw,nsav)
                    raercol(k,lmasscw,nsav)=0.
                 enddo
              enddo
           enddo
!      40 continue
        endif

     enddo OLD_CLOUD_MAIN_K_LOOP

!    cycle OVERALL_MAIN_I_LOOP


!    switch nsav, nnew so that nnew is the updated aerosol

     ntemp=nsav
     nsav=nnew
     nnew=ntemp

!    load new droplets in layers above, below clouds

     dtmin=dtstep
     ekk(kts)=0.0
! rce-comment -- ekd(k) is eddy-diffusivity at k/k-1 interface
!   want ekk(k) = ekd(k) * (density at k/k-1 interface)
     do k=kts+1,kte
        ekk(k)=ekd(k)*csbot(k)
     enddo
     ekk(kte+1)=0.0
     do k=kts,kte
        ekkp(k)=zn(k)*ekk(k+1)*zs(k+1)
        ekkm(k)=zn(k)*ekk(k)*zs(k)
        tinv=ekkp(k)+ekkm(k)
        if(k.eq.kts)tinv=tinv+surfratemax
        if(tinv.gt.1.e-6)then
           dtt=1./tinv
           dtmin=min(dtmin,dtt)
        endif
     enddo
     dtmix=0.9*dtmin
     nsubmix=dtstep/dtmix+1
     if(nsubmix>100)then
        nsubmix_bnd=100
     else
        nsubmix_bnd=nsubmix
     endif
!     count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
     dtmix=dtstep/nsubmix
     fac_srflx = -1.0/(zn(1)*nsubmix)
     
     do k=kts,kte
        kp1=min(k+1,kde-1)
        km1=max(k-1,1)
        if(lcldfra(kp1).gt.0)then
           overlapp(k)=min(lcldfra(k)/lcldfra(kp1),1.)
        else
           overlapp(k)=1.
        endif
        if(lcldfra(km1).gt.0)then
           overlapm(k)=min(lcldfra(k)/lcldfra(km1),1.)
        else
           overlapm(k)=1.
        endif
     enddo



! ......................................................................
! start of nsubmix-loop for calc of old cloud activation tendencies ....
OLD_CLOUD_NSUBMIX_LOOP: do nsub=1,nsubmix
        qndrop_new(kts:kte)=qndrop(kts:kte)
!           switch nsav, nnew so that nsav is the updated aerosol
        ntemp=nsav
        nsav=nnew
        nnew=ntemp
        srcn(:)=0.0
        do n=1,ntype_aer
           do m=1,nsize_aer(n)
              lnum=numptr_aer(m,n,ai_phase)
!              update droplet source
! rce-comment - activation source in layer k involves particles from k-1
!             srcn(kts  :kte)=srcn(kts  :kte)+nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
              srcn(kts+1:kte)=srcn(kts+1:kte)+nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
! rce-comment - new formulation for k=kts should be implemented
              srcn(kts      )=srcn(kts      )+nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
           enddo
        enddo
        call explmix(qndrop,srcn,ekkp,ekkm,overlapp,overlapm,   &
             qndrop_new,surfrate_drop,kms,kme,kts,kte,dtmix,.false.)
        do n=1,ntype_aer
           do m=1,nsize_aer(n)
              lnum=numptr_aer(m,n,ai_phase)
              lnumcw=numptr_aer(m,n,cw_phase)
              if(lnum>0)then
! rce-comment - activation source in layer k involves particles from k-1
!                source(kts  :kte)= nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
                 source(kts+1:kte)= nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
! rce-comment - new formulation for k=kts should be implemented
                 source(kts      )= nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
                 call explmix(raercol(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
                      raercol(1,lnumcw,nsav),surfrate(lnumcw),kms,kme,kts,kte,dtmix,&
                      .false.)
                 call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
                      raercol(1,lnum,nsav),surfrate(lnum),kms,kme,kts,kte,dtmix, &
                      .true.,raercol(1,lnumcw,nsav))
                 qsrflx(i,j,lnum) = qsrflx(i,j,lnum) + fac_srflx*            &
                      raercol(kts,lnum,nsav)*surfrate(lnum)
                 qsrflx(i,j,lnumcw) = qsrflx(i,j,lnumcw) + fac_srflx*        &
                      raercol(kts,lnumcw,nsav)*surfrate(lnumcw)
                 if (icheck_colmass > 0) then
                    tmpf = dtmix*rhodz(kts)
                    colmass_sfc(0,m,n,1) = colmass_sfc(0,m,n,1) &
                          + raercol(kts,lnum  ,nsav)*surfrate(lnum  )*tmpf
                    colmass_sfc(0,m,n,2) = colmass_sfc(0,m,n,2) &
                          + raercol(kts,lnumcw,nsav)*surfrate(lnumcw)*tmpf
                 endif
              endif
              do l=1,ncomp(n)
                 lmass=massptr_aer(l,m,n,ai_phase)
                 lmasscw=massptr_aer(l,m,n,cw_phase)
! rce-comment - activation source in layer k involves particles from k-1
!                source(kts  :kte)= mact(kts  :kte,m,n)*(raercol(kts:kte  ,lmass,nsav))
                 source(kts+1:kte)= mact(kts+1:kte,m,n)*(raercol(kts:kte-1,lmass,nsav))
! rce-comment - new formulation for k=kts should be implemented
                 source(kts      )= mact(kts      ,m,n)*(raercol(kts      ,lmass,nsav))
                 call explmix(raercol(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
                      raercol(1,lmasscw,nsav),surfrate(lmasscw),kms,kme,kts,kte,dtmix,  &
                      .false.)
                 call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
                      raercol(1,lmass,nsav),surfrate(lmass),kms,kme,kts,kte,dtmix,  &
                      .true.,raercol(1,lmasscw,nsav))
                 qsrflx(i,j,lmass) = qsrflx(i,j,lmass) + fac_srflx*          &
                      raercol(kts,lmass,nsav)*surfrate(lmass)
                 qsrflx(i,j,lmasscw) = qsrflx(i,j,lmasscw) + fac_srflx*      &
                      raercol(kts,lmasscw,nsav)*surfrate(lmasscw)
                 if (icheck_colmass > 0) then
                    ! colmass_sfc calculation
                    !    colmass_bgn/end = bgn/end column burden = sum.over.k.of{ rho(k)*dz(k)*chem(k,l) }
                    !    colmass_sfc = surface loss over dtstep
                    !       = sum.over.nsubmix.substeps{ depvel(l)*rho(kts)*chem(kts,l)*dtmix }
                    !    surfrate(l) = depvel(l)/dz(kts) so need to multiply by dz(kts)
                    !    for mass, raercol(k,l) = chem(k,l)*1.0e-9, so need to multiply by 1.0e9
                    tmpf = dtmix*rhodz(kts)*1.0e9
                    colmass_sfc(l,m,n,1) = colmass_sfc(l,m,n,1) &
                          + raercol(kts,lmass  ,nsav)*surfrate(lmass  )*tmpf
                    colmass_sfc(l,m,n,2) = colmass_sfc(l,m,n,2) &
                          + raercol(kts,lmasscw,nsav)*surfrate(lmasscw)*tmpf
                 endif
              enddo
              lwater=waterptr_aer(m,n)  ! aerosol water
              if(lwater>0)then
                 source(:)=0.
                 call explmix(   raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm,   &
                      raercol(1,lwater,nsav),surfrate(lwater),kms,kme,kts,kte,dtmix,  &
                      .true.,source)
              endif
           enddo ! size
        enddo ! type

     enddo OLD_CLOUD_NSUBMIX_LOOP

!    cycle OVERALL_MAIN_I_LOOP

!        evaporate particles again if no cloud

     do k=kts,kte
        if(lcldfra(k).eq.0.)then

!              no cloud

           qndrop(k)=0.
!              convert activated aerosol to interstitial in decaying cloud
           do n=1,ntype_aer
              do m=1,nsize_aer(n)
                 lnum=numptr_aer(m,n,ai_phase)
                 lnumcw=numptr_aer(m,n,cw_phase)
                 if(lnum.gt.0)then
                    raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew)
                    raercol(k,lnumcw,nnew)=0.
                 endif
                 do l=1,ncomp(n)
                    lmass=massptr_aer(l,m,n,ai_phase)
                    lmasscw=massptr_aer(l,m,n,cw_phase)
                    raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
                    raercol(k,lmasscw,nnew)=0.
                 enddo
              enddo
           enddo
        endif
     enddo

!         cycle OVERALL_MAIN_I_LOOP

!        droplet number

     do k=kts,kte
!       if(lcldfra(k).gt.0.1)then
!           write(6,'(a,3i5,f12.1)')'i,j,k,qndrop=',i,j,k,qndrop(k)
!       endif
        if(qndrop(k).lt.-10.e6.or.qndrop(k).gt.1.e12)then
           write(6,'(a,g12.2,a,3i5)')'after qndrop=',qndrop(k),' for i,k,j=',i,k,j
        endif

        qndrop3d(i,k,j) = max(qndrop(k),1.e-6)

        if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
           write(6,'(a,g12.2,a,3i5)')'after qndrop3d=',qndrop3d(i,k,j),' for i,k,j=',i,k,j
        endif
        if(qc(i,k,j).lt.-1..or.qc(i,k,j).gt.1.)then
           write(6,'(a,g12.2,a,3i5)')'qc=',qc(i,k,j),' for i,k,j=',i,k,j
           call wrf_error_fatal("1")
        endif
        if(qi(i,k,j).lt.-1..or.qi(i,k,j).gt.1.)then
           write(6,'(a,g12.2,a,3i5)')'qi=',qi(i,k,j),' for i,k,j=',i,k,j
           call wrf_error_fatal("1")
        endif
        if(qv(i,k,j).lt.-1..or.qv(i,k,j).gt.1.)then
           write(6,'(a,g12.2,a,3i5)')'qv=',qv(i,k,j),' for i,k,j=',i,k,j
           call wrf_error_fatal("1")
        endif
        cldfra_old(i,k,j) = cldfra(i,k,j)
!       if(k.gt.6.and.k.lt.11)cldfra_old(i,k,j)=1.
     enddo

!    cycle OVERALL_MAIN_I_LOOP

!        update chem and convert back to mole/mole

     ccn(:,:) = 0.
     do n=1,ntype_aer
        do m=1,nsize_aer(n)
           lnum=numptr_aer(m,n,ai_phase)
           lnumcw=numptr_aer(m,n,cw_phase)
           if(lnum.gt.0)then
              !          scale=mwdry*0.001
              scale = 1.
              chem(i,kts:kte,j,lnumcw)= raercol(kts:kte,lnumcw,nnew)*scale
              chem(i,kts:kte,j,lnum)= raercol(kts:kte,lnum,nnew)*scale
           endif
           do l=1,ncomp(n)
              lmass=massptr_aer(l,m,n,ai_phase)
              lmasscw=massptr_aer(l,m,n,cw_phase)
!          scale = mwdry/mw_aer(l,n)
              scale = 1.e9
              chem(i,kts:kte,j,lmasscw)=raercol(kts:kte,lmasscw,nnew)*scale ! ug/kg
              chem(i,kts:kte,j,lmass)=raercol(kts:kte,lmass,nnew)*scale ! ug/kg
           enddo
           lwater=waterptr_aer(m,n)
           if(lwater>0)chem(i,kts:kte,j,lwater)=raercol(kts:kte,lwater,nnew) ! don't convert units
           do k=kts,kte
              sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n)))
              do l=1,psat
                 arg=argfactor(m,n)*log(sm/super(l))
                 if(arg<2)then
                    if(arg<-2)then
                       ccnfact(l,m,n)=1.e-6 ! convert from #/m3 to #/cm3
                    else
                       ccnfact(l,m,n)=1.e-6*0.5*ERFC_NUM_RECIPES(arg)
                    endif
                 else
                    ccnfact(l,m,n) = 0.
                 endif
!                 ccn concentration as diagnostic
!                 assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles
                 ccn(k,l)=ccn(k,l)+(raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew))*cs(k)*ccnfact(l,m,n)
              enddo
           enddo
        enddo
     enddo
     do l=1,psat
        !wig, 22-Nov-2006: added vertical bounds to prevent out-of-bounds at top
        if(l.eq.1)ccn1(i,kts:kte,j)=ccn(:,l)
        if(l.eq.2)ccn2(i,kts:kte,j)=ccn(:,l)
        if(l.eq.3)ccn3(i,kts:kte,j)=ccn(:,l)
        if(l.eq.4)ccn4(i,kts:kte,j)=ccn(:,l)
        if(l.eq.5)ccn5(i,kts:kte,j)=ccn(:,l)
        if(l.eq.6)ccn6(i,kts:kte,j)=ccn(:,l)
     end do

! mass conservation checking
     if (icheck_colmass > 0) then
! calc final column burdens
        do n=1,ntype_aer
        do m=1,nsize_aer(n)
           lnum=numptr_aer(m,n,ai_phase)
           lnumcw=numptr_aer(m,n,cw_phase)
           if(lnum>0)then
              colmass_end(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
              colmass_end(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
           endif
           do l=1,ncomp(n)
              lmass=massptr_aer(l,m,n,ai_phase)
              lmasscw=massptr_aer(l,m,n,cw_phase)
              colmass_end(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
              colmass_end(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
           enddo
        enddo ! size
        enddo ! type
! calc burden change errors for each interstitial/activated pair
        do n=1,ntype_aer
        do m=1,nsize_aer(n)
           do l=0,ncomp(n)
     …