phys/module_cam_mp_ndrop.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   4       module ndrop
   5 #ifdef MODAL_AERO
   6       use shr_kind_mod,  only: r8 => shr_kind_r8
   7 #ifndef WRF_PORT
   8       use abortutils,    only: endrun
   9       use modal_aero_data, only: ntot_amode, qqcw_get_field
  10       use cam_logfile,     only: iulog
  11 #else
  12       use module_cam_support, only: endrun, iulog
  13       use modal_aero_data, only: ntot_amode
  14
  15 #endif
  16
  17       implicit none
  18
  19       private
  20       save
  21
  22       public activate_init, dropmixnuc
  23
  24       real(r8) :: npv(ntot_amode) ! number per volume concentration
  25       real(r8) :: alogsig(ntot_amode) ! natl log of geometric standard dev of aerosol
  26       real(r8) :: exp45logsig(ntot_amode)
  27       real(r8) :: argfactor(ntot_amode)
  28       real(r8) :: f1(ntot_amode),f2(ntot_amode)  ! abdul-razzak functions of width
  29       real(r8) :: t0 ! reference temperature
  30       real(r8) :: aten
  31       real(r8) :: surften       ! surface tension of water w/respect to air (N/m)
  32       real(r8) :: alogten,alog2,alog3,alogaten
  33       real(r8) :: third, twothird, sixth, zero
  34       real(r8) :: sq2, sqpi, pi
  35
  36       type qqcw_type
  37          real(r8), pointer :: fldcw(:,:)
  38       end type qqcw_type
  39
  40 contains
  41
  42       subroutine dropmixnuc(lchnk, ncol, ncldwtr,tendnd, temp,omega,  &
  43                     pmid,pint,pdel,rpdel,zm,kvh,wsub,cldn,cldo,     &
  44                     raer, cflx, raertend, dtmicro                   &
  45 #ifdef WRF_PORT
  46                     , qqcw_inout                                    &
  47 #endif
  48                                                                     )
  49
  50 !     vertical diffusion and nucleation of cloud droplets
  51 !     assume cloud presence controlled by cloud fraction
  52 !     doesn't distinguish between warm, cold clouds
  53 #ifndef WRF_PORT
  54       use ppgrid,        only: pcols, pver, pverp
  55 #else
  56       use module_cam_support, only: pcols, pver, pverp
  57 #endif
  58       use physconst,     only: gravit, rhoh2o, latvap, cpair, epsilo, rair, mwh2o, r_universal
  59 #ifndef WRF_PORT
  60       use time_manager,  only: get_nstep
  61       use constituents,  only: pcnst, qmin, cnst_get_ind, cnst_name
  62 #else
  63       use module_cam_support, only: pcnst => pcnst_mp
  64       use constituents,  only: cnst_get_ind, cnst_name
  65 #endif
  66       use error_function, only: erf
  67
  68 !      use aerosol_radiation_interface, only: collect_sw_aer_masses, aer_mass
  69
  70 #ifndef WRF_PORT
  71       use modal_aero_data
  72       use cam_history, only: fieldname_len
  73       use phys_grid,   only: get_lat_all_p, get_lon_all_p
  74       use physics_types, only: physics_state
  75       use cam_history,   only: outfld
  76 #else
  77       use modal_aero_data,    only: numptrcw_amode => numptrcw_amode_mp, &
  78            nspec_amode  => nspec_amode_mp, lmassptrcw_amode => lmassptrcw_amode_mp, &
  79            numptr_amode => numptr_amode_mp  , lmassptr_amode   => lmassptr_amode_mp, &
  80            cnst_name_cw => cnst_name_cw_mp
  81       use module_cam_support, only: fieldname_len, outfld
  82 #endif
  83
  84       implicit none
  85
  86 !     input
  87
  88       integer, intent(in) :: lchnk                ! chunk identifier
  89       integer, intent(in) :: ncol                 ! number of columns
  90 !      type(physics_state), intent(in) :: state      ! Physics state variables
  91       real(r8), intent(in) :: dtmicro             ! time step for microphysics (s)
  92       real(r8), intent(in) :: temp(pcols,pver)    ! temperature (K)
  93       real(r8), intent(inout) :: ncldwtr(pcols,pver) ! droplet number concentration (#/kg)
  94       real(r8), intent(inout) :: tendnd(pcols,pver) ! change in droplet number concentration (#/kg/s)
  95       real(r8), intent(in) :: omega(pcols,pver)   ! vertical velocity (Pa/s)
  96       real(r8), intent(in) :: pmid(pcols,pver)    ! mid-level pressure (Pa)
  97       real(r8), intent(in) :: pint(pcols,pverp)   ! pressure at layer interfaces (Pa)
  98       real(r8), intent(in) :: pdel(pcols,pver)    ! pressure thickess of layer (Pa)
  99       real(r8), intent(in) :: rpdel(pcols,pver)   ! inverse of pressure thickess of layer (/Pa)
 100       real(r8), intent(in) :: zm(pcols,pver)      ! geopotential height of level (m)
 101       real(r8), intent(in) :: kvh(pcols,pverp)    ! vertical diffusivity (m2/s)
 102       real(r8), intent(in) :: wsub(pcols,pver)    ! subgrid vertical velocity
 103       real(r8), intent(in) :: cldo(pcols,pver)    ! cloud fraction on previous time step
 104       real(r8), intent(in) :: cldn(pcols,pver)    ! cloud fraction
 105
 106
 107 !     output
 108
 109       real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
 110       real(r8), intent(in) :: cflx(pcols,pcnst) ! constituent flux from surface
 111       real(r8), intent(out) :: raertend(pcols,pver,pcnst) ! tendency of aerosol mass, number mixing ratios
 112 #ifdef WRF_PORT
 113       real(r8), target, intent(inout) :: qqcw_inout(pcols,pver,pcnst) ! cloud-borne aerosol
 114 #endif
 115
 116 !--------------------Local storage-------------------------------------
 117 !
 118       type(qqcw_type) :: QQCW(pcnst)
 119
 120       real(r8) depvel(pcols,pcnst)! deposition velocity for droplets (m/s)
 121       real(r8) qcld(pver) ! cloud droplet number mixing ratio (#/kg)
 122       real(r8) wtke(pcols,pver)     ! turbulent vertical velocity at base of layer k (m/s)
 123       real(r8) wtke_cen(pcols,pver) ! turbulent vertical velocity at center of layer k (m/s)
 124       real(r8) zn(pver) ! g/pdel (m2/g) for layer
 125       real(r8) zs(pver) ! inverse of distance between levels (m)
 126       real(r8), parameter :: zkmin=0.01_r8,zkmax=100._r8
 127       real(r8) w(pcols,pver)       ! vertical velocity (m/s)
 128       real(r8) cs(pcols,pver)      ! air density (kg/m3)
 129       real(r8) dz(pcols,pver)      ! geometric thickness of layers (m)
 130       real(r8) zero,flxconv ! convergence of flux into lowest layer
 131
 132       real(r8) wdiab           ! diabatic vertical velocity
 133       real(r8), parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
 134 !       real(r8), parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
 135 !      real(r8), parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
 136       real(r8) qncld(pver)     ! droplet number nucleated on cloud boundaries
 137       real(r8) ekd(pver)       ! diffusivity for droplets (m2/s)
 138       real(r8) ekk(0:pver)       ! density*diffusivity for droplets (kg/m3 m2/s)
 139       real(r8) srcn(pver) ! droplet source rate (/s)
 140       real(r8), parameter :: sq2pi=2.5066283_r8
 141       real(r8) dtinv
 142
 143       logical top        ! true if cloud top, false if cloud base or new cloud
 144       integer km1,kp1
 145       real(r8) wbar,wmix,wmin,wmax
 146       real(r8) dum,dumc
 147       real(r8) dact
 148       real(r8) fluxntot         ! (#/cm2/s)
 149       real(r8) fac_srflx
 150       real(r8) surfrate(pcnst) ! surface exchange rate (/s)
 151       real(r8) surfratemax      ! max surfrate for all species treated here
 152       real(r8) dtmin,tinv,dtt
 153       integer nsubmix,nsubmix_bnd
 154       integer i,k,m,n
 155       real(r8) tempr4 ! temperature
 156       real(r8) dtmix
 157       real(r8) alogarg
 158       integer nstep
 159       real(r8) pi
 160       integer nnew,nsav,ntemp
 161       real(r8) overlapp(pver),overlapm(pver) ! cloud overlap
 162       real(r8) ekkp(pver),ekkm(pver) ! zn*zs*density*diffusivity
 163       integer count_submix(100)
 164       save count_submix
 165       real(r8) kvhmax
 166       save kvhmax
 167       real(r8) nsource(pcols,pver)            ! droplet number source (#/kg/s)
 168       real(r8) ndropmix(pcols,pver)           ! droplet number mixing (#/kg/s)
 169       real(r8) ndropcol(pcols)               ! column droplet number (#/m2)
 170
 171       integer lnum, lnumcw, lmass, lmasscw, lphase, &
 172               lsfc, lsfccw, lsig, lspec, ltype, lwater
 173       integer ntype(ntot_amode)
 174
 175       real(r8) na(pcols),va(pcols),hy(pcols)
 176       real(r8) naermod(ntot_amode) ! (/m3)
 177       real(r8) sigmag(ntot_amode)  ! geometric standard deviation of aerosol size dist
 178       real(r8) hygro(ntot_amode)  ! hygroscopicity of aerosol mode
 179       real(r8) vaerosol(ntot_amode) ! interstit+activated aerosol volume conc (cm3/cm3)
 180       real(r8) raercol(pver,pcnst,2) ! aerosol mass, number mixing ratios
 181       real(r8) raercol_cw(pver,pcnst,2) ! cloud-borne aerosol mass, number mixing ratios
 182       real(r8) surfrate_cw(pcnst) ! surface exchange rate for cloud-borne aerosol (/s)
 183       real(r8) source(pver) !
 184
 185       real(r8) fn(ntot_amode)         ! activation fraction for aerosol number
 186       real(r8) fm(ntot_amode)         ! activation fraction for aerosol mass
 187
 188       real(r8) fluxn(ntot_amode)      ! number  activation fraction flux (cm/s)
 189       real(r8) fluxm(ntot_amode)      ! mass    activation fraction flux (cm/s)
 190       real(r8) sum,sumcw,flux,fluxcw
 191 !     note:  activation fraction fluxes are defined as
 192 !     fluxn = [flux of activated aero. number into cloud (#/cm2/s)]
 193 !           / [aero. number conc. in updraft, just below cloudbase (#/cm3)]
 194
 195       real(r8) nact(pver,ntot_amode)  ! fractional aero. number  activation rate (/s)
 196       real(r8) mact(pver,ntot_amode)  ! fractional aero. mass    activation rate (/s)
 197       real(r8) sigmag_amode_cur(ntot_amode) ! sigmag for current grid cell
 198       save sigmag_amode_cur
 199       real(r8) :: qqcwtend(pcols,pver,pcnst) ! tendency of cloudborne aerosol mass, number mixing ratios
 200       real(r8) :: coltend(pcols)    ! column tendency
 201       real(r8) :: tmpa
 202       integer :: lc
 203       character(len=fieldname_len)   :: tmpname
 204       character(len=fieldname_len+3) :: fieldname
 205       real(r8) :: csbot(pver)       ! air density at bottom (interface) of layer (kg/m3)
 206       real(r8) :: csbot_cscen(pver) ! csbot(i)/cs(i,k)
 207       real(r8) :: flux_fullact(pver)     ! 100%    activation fraction flux (cm/s)
 208       real(r8) :: taumix_internal_pver_inv ! 1/(internal mixing time scale for k=pver) (1/s)
 209       real(r8) dactn,taunuc,damp
 210       real(r8) :: cldo_tmp, cldn_tmp
 211       real(r8) :: tau_cld_regenerate
 212
 213       logical cldinterior
 214
 215       integer ixndrop, l
 216       real(r8) naer_tot,naer(pcols)
 217       integer, parameter :: psat=6 ! number of supersaturations to calc ccn concentration
 218       real(r8)  :: supersat(psat)= & ! supersaturation (%) to determine ccn concentration
 219                (/0.02,0.05,0.1,0.2,0.5,1.0/)
 220       real(r8) ccn(pcols,pver,psat)        ! number conc of aerosols activated at supersat
 221       character(len=8), dimension(psat) :: ccn_name(psat)= &
 222                (/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
 223       real(r8) arg
 224       integer phase ! phase of aerosol
 225
 226
 227
 228     arg = 1.0_r8
 229     if (abs(0.8427_r8-erf(arg))/0.8427_r8>0.001_r8) then
 230        write (iulog,*) 'erf(1.0) = ',ERF(arg)
 231 #ifdef WRF_PORT
 232        call wrf_message(iulog)
 233 #endif
 234        call endrun('dropmixnuc: Error function error')
 235     endif
 236     arg = 0.0
 237     if (erf(arg) /= 0.0) then
 238        write (iulog,*) 'erf(0.0) = ',erf(arg)
 239 #ifdef WRF_PORT
 240        call wrf_message(iulog)
 241 #endif
 242        write (iulog,*) 'dropmixnuc: Error function error'
 243 #ifdef WRF_PORT
 244        call wrf_message(iulog)
 245 #endif
 246        call endrun('dropmixnuc: Error function error')
 247     endif
 248      zero=0._r8
 249
 250
 251        pi = 4._r8*atan(1.0_r8)
 252        dtinv=1./dtmicro
 253 #ifndef WRF_PORT
 254        nstep = get_nstep()
 255 #endif
 256
 257 !       call t_startf('load_aerosol')
 258
 259
 260        call cnst_get_ind('NUMLIQ', ixndrop)
 261        depvel(:,:) = 0.0_r8        ! droplet number is done in pkg_cld_sediment, aerosols in mz_aerosols_intr
 262 !      depvel(:,ixndrop) =  0.1_r8 ! prescribed here rather than getting it from MIRAGE depvel_part
 263
 264 !    call t_stopf('load_aerosol')
 265
 266        do m=1,ntot_amode
 267 #ifndef WRF_PORT
 268           QQCW(numptrcw_amode(m))%fldcw => qqcw_get_field(numptrcw_amode(m),lchnk)
 269 #else
 270           QQCW(numptrcw_amode(m))%fldcw => qqcw_inout(:,:,numptrcw_amode(m))
 271 #endif
 272           do l=1,nspec_amode(m)
 273 #ifndef WRF_PORT
 274              QQCW(lmassptrcw_amode(l,m))%fldcw => qqcw_get_field(lmassptrcw_amode(l,m),lchnk)
 275 #else
 276              QQCW(lmassptrcw_amode(l,m))%fldcw => qqcw_inout(:,:,lmassptrcw_amode(l,m))
 277 #endif
 278           end do
 279        end do
 280
 281
 282 !     start loop over columns
 283 overall_main_i_loop: &
 284       do i=1,ncol
 285
 286 !        load number nucleated into qcld on cloud boundaries
 287
 288 ! initialization for current i .........................................
 289
 290 !    call t_startf('calc_wtke')
 291          do k=1,pver-1
 292             zs(k)=1._r8/(zm(i,k)-zm(i,k+1))
 293          enddo
 294          zs(pver)=zs(pver-1)
 295
 296          do k=1,pver
 297             qcld(k)=ncldwtr(i,k)
 298             qncld(k)=0._r8
 299             srcn(k)=0._r8
 300             cs(i,k)=pmid(i,k)/(rair*temp(i,k)) ! air density (kg/m3)
 301             dz(i,k)=1._r8/(cs(i,k)*gravit*rpdel(i,k)) ! layer thickness in m
 302             w(i,k)=-1._r8*omega(i,k)/(cs(i,k)*gravit) ! large-scale vertical velocity m/s
 303             do m=1,ntot_amode
 304                nact(k,m)=0._r8
 305                mact(k,m)=0._r8
 306             enddo
 307             zn(k)=gravit*rpdel(i,k)
 308             kvhmax=max(kvh(i,k),kvhmax)
 309             if(k<pver)then
 310                ekd(k)=kvh(i,k+1)
 311                ekd(k)=max(ekd(k),zkmin)
 312                ekd(k)=min(ekd(k),zkmax)
 313                csbot(k)=2.0_r8*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
 314                csbot_cscen(k) = csbot(k)/cs(i,k)
 315             else
 316                ekd(k)=0._r8
 317                csbot(k)=cs(i,k)
 318                csbot_cscen(k) = 1.0_r8
 319             endif
 320 !           diagnose subgrid vertical velocity from diffusivity
 321             if(k.eq.pver)then
 322                wtke(i,k)=sq2pi*depvel(i,ixndrop)
 323 !               wtke(i,k)=sq2pi*kvh(i,k+1)
 324                wtke(i,k)=max(wtke(i,k),wmixmin)
 325             else
 326                wtke(i,k)=sq2pi*ekd(k)/dz(i,k)
 327             endif
 328 ! get sub-grid vertical velocity from diff coef.
 329 ! following morrison et al. 2005, JAS
 330 ! assume mixing length of 30 m
 331 ! rce-comment - define wtke at layer centers for new-cloud activation
 332 !    and at layer boundaries for old-cloud activation
 333
 334 !++ag
 335 !            wtke_cen(i,k)=0.5_r8*(kvh(i,k)+kvh(i,k+1))/30._r8
 336 !            wtke(i,k)=kvh(i,k+1)/30._r8
 337             wtke_cen(i,k)=wsub(i,k)
 338             wtke(i,k)=wsub(i,k)
 339 !--ag
 340             wtke_cen(i,k)=max(wtke_cen(i,k),wmixmin)
 341             wtke(i,k)=max(wtke(i,k),wmixmin)
 342             nsource(i,k)=0._r8
 343          enddo
 344
 345             surfratemax = 0.0_r8
 346             nsav=1
 347             nnew=2
 348             surfrate(ixndrop)=depvel(i,ixndrop)/dz(i,pver)
 349             surfratemax = max( surfratemax, surfrate(ixndrop) )
 350             do m=1,ntot_amode
 351                lnum=numptr_amode(m)
 352                lnumcw=numptrcw_amode(m)
 353                if(lnum>0)then
 354                   surfrate(lnum)=depvel(i,lnum)/dz(i,pver)
 355 !                 surfrate(lnumcw)=surfrate(ixndrop)
 356                   surfrate_cw(lnumcw)=surfrate(ixndrop)
 357                   surfratemax = max( surfratemax, surfrate(lnum) )
 358 !                 raercol(:,lnumcw,nsav)=raer(i,:,lnumcw)
 359                   raercol_cw(:,lnumcw,nsav)=qqcw(lnumcw)%fldcw(i,:)   ! use cloud-borne aerosol array
 360                   raercol(:,lnum,nsav)=raer(i,:,lnum)
 361                endif
 362                do l=1,nspec_amode(m)
 363                   lmass=lmassptr_amode(l,m)
 364                   lmasscw=lmassptrcw_amode(l,m)
 365                   surfrate(lmass)=depvel(i,lmass)/dz(i,pver)
 366 !                 surfrate(lmasscw)=surfrate(ixndrop)
 367                   surfrate_cw(lmasscw)=surfrate(ixndrop)
 368                   surfratemax = max( surfratemax, surfrate(lmass) )
 369 !                 raercol(:,lmasscw,nsav)=raer(i,:,lmasscw)
 370                   raercol_cw(:,lmasscw,nsav)=qqcw(lmasscw)%fldcw(i,:)  ! use cloud-borne aerosol array
 371                   raercol(:,lmass,nsav)=raer(i,:,lmass)
 372                enddo
 373 !              lwater=lwaterptr_amode(m)
 374 !              surfrate(lwater)=depvel(i,lwater)/dz(i,pver)
 375 !              surfratemax = max( surfratemax, surfrate(lwater) )
 376 !              raercol(:,lwater,nsav)=raer(i,:,lwater)
 377             enddo
 378 !    call t_stopf('calc_wtke')
 379
 380 !        droplet nucleation/aerosol activation
 381
 382 ! tau_cld_regenerate = time scale for regeneration of cloudy air
 383 !    by (horizontal) exchange with clear air
 384             tau_cld_regenerate = 3600.0_r8 * 3.0_r8
 385
 386 ! k-loop for growing/shrinking cloud calcs .............................
 387 grow_shrink_main_k_loop: &
 388          do k=1,pver
 389             km1=max0(k-1,1)
 390             kp1=min0(k+1,pver)
 391
 392 ! shrinking cloud ......................................................
 393 !    treat the reduction of cloud fraction from when cldn(i,k) < cldo(i,k)
 394 !    and also dissipate the portion of the cloud that will be regenerated
 395               cldo_tmp = cldo(i,k)
 396               cldn_tmp = cldn(i,k) * exp( -dtmicro/tau_cld_regenerate )
 397 !    alternate formulation
 398 !             cldn_tmp = cldn(i,k) * max( 0.0_r8, (1.0_r8-dtmicro/tau_cld_regenerate) )
 399
 400               if(cldn_tmp.lt.cldo_tmp)then
 401 !                droplet loss in decaying cloud
 402 !++ sungsup
 403                  nsource(i,k)=nsource(i,k)+qcld(k)*(cldn_tmp-cldo_tmp)/cldo_tmp*dtinv
 404                  qcld(k)=qcld(k)*(1.+(cldn_tmp-cldo_tmp)/cldo_tmp)
 405 !-- sungsup
 406
 407 !                 convert activated aerosol to interstitial in decaying cloud
 408
 409                   dumc=(cldn_tmp-cldo_tmp)/cldo_tmp
 410                   do m=1,ntot_amode
 411                      lnum=numptr_amode(m)
 412                      lnumcw=numptrcw_amode(m)
 413                      if(lnum.gt.0)then
 414                         dact=raercol_cw(k,lnumcw,nsav)*dumc
 415 !                       raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
 416                         raercol_cw(k,lnumcw,nsav)=raercol_cw(k,lnumcw,nsav)+dact   ! cloud-borne aerosol
 417                         raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
 418                      endif
 419                      do l=1,nspec_amode(m)
 420                         lmass=lmassptr_amode(l,m)
 421                         lmasscw=lmassptrcw_amode(l,m)
 422                         dact=raercol_cw(k,lmasscw,nsav)*dumc
 423 !                       raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
 424                         raercol_cw(k,lmasscw,nsav)=raercol_cw(k,lmasscw,nsav)+dact  ! cloud-borne aerosol
 425                         raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
 426                      enddo
 427                   enddo
 428                endif
 429
 430 ! growing cloud ......................................................
 431 !    treat the increase of cloud fraction from when cldn(i,k) > cldo(i,k)
 432 !    and also regenerate part of the cloud
 433               cldo_tmp = cldn_tmp
 434               cldn_tmp = cldn(i,k)
 435
 436               if(cldn_tmp-cldo_tmp.gt.0.01)then
 437 !                wmix=wtke(i,k)
 438 !                wbar=wtke(i,k)
 439 ! rce-comment - use wtke at layer centers for new-cloud activation
 440                  wbar=wtke_cen(i,k)
 441                  wmix=0._r8
 442                  wmin=0._r8
 443                  wmax=10._r8
 444                  wdiab=0
 445 !                load aerosol properties, assuming external mixtures
 446
 447                 phase=1 ! interstitial
 448                  do m=1,ntot_amode
 449                     call loadaer(raer,qqcw,i,i,k,m,cs,npv(m),phase, &
 450                          na, va,  hy )
 451                          naermod(m)=na(i)
 452                          vaerosol(m)=va(i)
 453                          hygro(m)=hy(i)
 454                  end do
 455 !                m=1
 456 !                naermod(m)=1000.e6
 457 !                vaerosol(m)=naermod(m)/(density*num_to_mass_aer(m))
 458
 459 !                call t_startf ('activate')
 460
 461                  call activate_modal(wbar,wmix,wdiab,wmin,wmax,temp(i,k),cs(i,k), &
 462                       naermod, ntot_amode,ntot_amode, vaerosol, sigmag,hygro,       &
 463                       fn,fm,fluxn,fluxm,flux_fullact(k))
 464
 465 !                call t_stopf ('activate')
 466 !                write(iulog,'(a,15f8.2)')'wtke,naer,fn=',wtke(i,k),naer_tot*1.e-6,(fn(m),m=1,ntot_amode)
 467
 468
 469                  dumc=(cldn_tmp-cldo_tmp)
 470                  do m=1,ntot_amode
 471                     lnum=numptr_amode(m)
 472                     lnumcw=numptrcw_amode(m)
 473                     dact=dumc*fn(m)*raer(i,k,lnum) ! interstitial only
 474                     qcld(k)=qcld(k)+dact
 475                     nsource(i,k)=nsource(i,k)+dact*dtinv
 476                     if(lnum.gt.0)then
 477                        raercol_cw(k,lnumcw,nsav)=raercol_cw(k,lnumcw,nsav)+dact  ! cloud-borne aerosol
 478                        raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
 479                     endif
 480                     dum=dumc*fm(m)
 481                     do l=1,nspec_amode(m)
 482                         lmass=lmassptr_amode(l,m)
 483                         lmasscw=lmassptrcw_amode(l,m)
 484                         dact=dum*(raer(i,k,lmass)) ! interstitial only
 485                         raercol_cw(k,lmasscw,nsav)=raercol_cw(k,lmasscw,nsav)+dact  ! cloud-borne aerosol
 486                         raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
 487                     enddo
 488                  enddo
 489                endif
 490
 491          enddo grow_shrink_main_k_loop
 492 ! end of k-loop for growing/shrinking cloud calcs ......................
 493
 494
 495 ! ......................................................................
 496 ! start of k-loop for calc of old cloud activation tendencies ..........
 497 !
 498 ! rce-comment
 499 !    changed this part of code to use current cloud fraction (cldn) exclusively
 500 !    consider case of cldo(:)=0, cldn(k)=1, cldn(k+1)=0
 501 !    previous code (which used cldo below here) would have no cloud-base activation
 502 !       into layer k.  however, activated particles in k mix out to k+1,
 503 !       so they are incorrectly depleted with no replacement
 504 !
 505 old_cloud_main_k_loop: &
 506          do k=1,pver
 507             km1=max0(k-1,1)
 508             kp1=min0(k+1,pver)
 509             taumix_internal_pver_inv = 0.0
 510
 511             if(cldn(i,k).gt.0.01)then
 512 ! rce-comment
 513 !             if(cldo(i,k).gt.0.01)then   ! with cldo changed to cldn this is redundant
 514                wdiab=0
 515 !               wmix=wtke(i,k) ! spectrum of updrafts
 516 !              wbar=w(i,k) ! spectrum of updrafts
 517                wmix=0._r8 ! single updraft
 518                wbar=wtke(i,k) ! single updraft
 519 ! rce-comment - different treatment for k=pver
 520                if (k == pver) wbar=wtke_cen(i,k) ! single updraft
 521                wmax=10._r8
 522                wmin=0._r8
 523 !              old cloud
 524                   if(cldn(i,k)-cldn(i,kp1).gt.0.01.or.k.eq.pver)then
 525
 526 !                   cloud base
 527
 528                     ekd(k)=wtke(i,k)*dz(i,k)/sq2pi
 529 ! rce-comments
 530 !   first, should probably have 1/zs(k) here rather than dz(i,k) because
 531 !      the turbulent flux is proportional to ekd(k)*zs(k),
 532 !      while the dz(i,k) is used to get flux divergences
 533 !      and mixing ratio tendency/change
 534 !   second and more importantly, using a single updraft velocity here
 535 !      means having monodisperse turbulent updraft and downdrafts.
 536 !      The sq2pi factor assumes a normal draft spectrum.
 537 !      The fluxn/fluxm from activate must be consistent with the
 538 !      fluxes calculated in explmix.
 539                     ekd(k)=wbar/zs(k)
 540
 541                     alogarg=max(1.e-20_r8,1/cldn(i,k)-1._r8)
 542                     wmin=wbar+wmix*0.25*sq2pi*log(alogarg)
 543                     phase=1 ! interstitial
 544                     do m=1,ntot_amode
 545 ! rce-comment - use kp1 here as old-cloud activation involves
 546 !   aerosol from layer below
 547 !                      call loadaer(raer,qqcw,i,i,k,m,cs, npv(m),phase, &
 548                        call loadaer(raer,qqcw,i,i,kp1,m,cs, npv(m),phase, &
 549                          na, va,  hy )
 550                          naermod(m)=na(i)
 551                          vaerosol(m)=va(i)
 552                          hygro(m)=hy(i)
 553                     end do
 554 !                   m=1
 555 !                   naermod(m)=1000.e6
 556 !                   vaerosol(1,m)=naermod(m)/(density*num_to_mass_aer(m))
 557 !                   call t_startf ('activate')
 558
 559                     call activate_modal(wbar,wmix,wdiab,wmin,wmax,temp(i,k),cs(i,k), &
 560                       naermod, ntot_amode,ntot_amode, vaerosol, sigmag, hygro,    &
 561                       fn,fm,fluxn,fluxm,flux_fullact(k))
 562
 563 !                   call t_stopf ('activate')
 564 !                   write(iulog,'(a,15f8.2)')'wtke,naer,fn=',wtke(i,k),naer_tot*1.e-6,(fn(m),m=1,ntot_amode)
 565
 566                       if(k.lt.pver)then
 567                           dumc = cldn(i,k)-cldn(i,kp1)
 568                       else
 569                           dumc=cldn(i,k)
 570                       endif
 571                     fluxntot=0
 572 ! rce-comment 1
 573 !    flux of activated mass into layer k (in kg/m2/s)
 574 !       = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
 575 !    source of activated mass (in kg/kg/s) = flux divergence
 576 !       = actmassflux/(cs(i,k)*dz(i,k))
 577 !    so need factor of csbot_cscen = csbot(k)/cs(i,k)
 578 !                   dum=1./(dz(i,k))
 579                     dum=csbot_cscen(k)/(dz(i,k))
 580 ! rce-comment 2
 581 !    code for k=pver was changed to use the following conceptual model
 582 !    in k=pver, there can be no cloud-base activation unless one considers
 583 !       a scenario such as the layer being partially cloudy,
 584 !       with clear air at bottom and cloudy air at top
 585 !    assume this scenario, and that the clear/cloudy portions mix with
 586 !       a timescale taumix_internal = dz(i,pver)/wtke_cen(i,pver)
 587 !    in the absence of other sources/sinks, qact (the activated particle
 588 !       mixratio) attains a steady state value given by
 589 !          qact_ss = fcloud*fact*qtot
 590 !       where fcloud is cloud fraction, fact is activation fraction,
 591 !       qtot=qact+qint, qint is interstitial particle mixratio
 592 !    the activation rate (from mixing within the layer) can now be
 593 !       written as
 594 !          d(qact)/dt = (qact_ss - qact)/taumix_internal
 595 !                     = qtot*(fcloud*fact*wtke/dz) - qact*(wtke/dz)
 596 !    note that (fcloud*fact*wtke/dz) is equal to the nact/mact
 597 !    also, d(qact)/dt can be negative.  in the code below
 598 !       it is forced to be >= 0
 599 !
 600 ! steve --
 601 !    you will likely want to change this.  i did not really understand
 602 !       what was previously being done in k=pver
 603 !    in the cam3_5_3 code, wtke(i,pver) appears to be equal to the
 604 !       droplet deposition velocity which is quite small
 605 !    in the cam3_5_37 version, wtke is done differently and is much
 606 !       larger in k=pver, so the activation is stronger there
 607 !
 608                     if (k == pver) then
 609                        taumix_internal_pver_inv = flux_fullact(k)/dz(i,k)
 610                     end if
 611                     do m=1,ntot_amode
 612                        fluxn(m)=fluxn(m)*dumc
 613                        fluxm(m)=fluxm(m)*dumc
 614                        lnum=numptr_amode(m)
 615                        nact(k,m)=nact(k,m)+fluxn(m)*dum
 616                        mact(k,m)=mact(k,m)+fluxm(m)*dum
 617                        if (k > pver) then
 618 ! note that kp1 is used here
 619                           fluxntot = fluxntot &
 620                                    +fluxn(m)*raercol(kp1,lnum,nsav)*cs(i,k)
 621                        else
 622                           lnumcw=numptrcw_amode(m)
 623                           tmpa = raercol(kp1,lnum,nsav)*fluxn(m) &
 624                                + raercol(kp1,lnumcw,nsav)*(fluxn(m) &
 625                                     -taumix_internal_pver_inv)
 626                           fluxntot = fluxntot + max(0.0_r8,tmpa)*cs(i,k)
 627                        end if
 628                     enddo
 629                       srcn(k)=srcn(k)+fluxntot/(cs(i,k)*dz(i,k))
 630                       nsource(i,k)=nsource(i,k)+fluxntot/(cs(i,k)*dz(i,k))
 631
 632                  endif  ! (cldn(i,k)-cldn(i,kp1).gt.0.01 .or. k.eq.pver)
 633
 634 ! rce-comment
 635 !             endif ! (cldo(i,k).gt.0.01)   ! with cldo changed to cldn this is redundant
 636
 637             else
 638 !              no cloud
 639                nsource(i,k)=nsource(i,k)-qcld(k)*dtinv
 640                qcld(k)=0
 641 !              convert activated aerosol to interstitial in decaying cloud
 642                do m=1,ntot_amode
 643                   lnum=numptr_amode(m)
 644                   lnumcw=numptrcw_amode(m)
 645                   if(lnum.gt.0)then
 646                      raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol_cw(k,lnumcw,nsav)  ! cloud-borne aerosol
 647                      raercol_cw(k,lnumcw,nsav)=0.
 648                   endif
 649                   do l=1,nspec_amode(m)
 650                      lmass=lmassptr_amode(l,m)
 651                      lmasscw=lmassptrcw_amode(l,m)
 652                      raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol_cw(k,lmasscw,nsav) ! cloud-borne aerosol
 653                      raercol_cw(k,lmasscw,nsav)=0.
 654                   enddo
 655                 enddo
 656             endif
 657          enddo old_cloud_main_k_loop
 658 !           switch nsav, nnew so that nnew is the updated aerosol
 659             ntemp=nsav
 660             nsav=nnew
 661             nnew=ntemp
 662 !          call t_startf ('nsubmix')
 663
 664 !        load new droplets in layers above, below clouds
 665
 666          dtmin=dtmicro
 667          ekk(0)=0.0
 668          ekk(pver)=0.0
 669          do k=1,pver-1
 670 ! rce-comment -- ekd(k) is eddy-diffusivity at k/k+1 interface
 671 !   want ekk(k) = ekd(k) * (density at k/k+1 interface)
 672 !   so use pint(i,k+1) as pint is 1:pverp
 673 !           ekk(k)=ekd(k)*2.*pint(i,k)/(rair*(temp(i,k)+temp(i,k+1)))
 674 !           ekk(k)=ekd(k)*2.*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
 675             ekk(k)=ekd(k)*csbot(k)
 676          enddo
 677          do k=1,pver
 678             km1=max0(k-1,1)
 679             ekkp(k)=zn(k)*ekk(k)*zs(k)
 680             ekkm(k)=zn(k)*ekk(k-1)*zs(km1)
 681             tinv=ekkp(k)+ekkm(k)
 682
 683             if(k.eq.pver)tinv=tinv+surfratemax
 684 ! rce-comment -- tinv is the sum of all first-order-loss-rates
 685 !    for the layer.  for most layers, the activation loss rate
 686 !    (for interstitial particles) is accounted for by the loss by
 687 !    turb-transfer to the layer above.
 688 !    k=pver is special, and the loss rate for activation within
 689 !    the layer must be added to tinv.  if not, the time step
 690 !    can be too big, and explmix can produce negative values.
 691 !    the negative values are reset to zero, resulting in an
 692 !    artificial source.
 693             if(k.eq.pver)tinv=tinv+taumix_internal_pver_inv
 694
 695             if(tinv.gt.1.e-6)then
 696                dtt=1./tinv
 697                dtmin=min(dtmin,dtt)
 698             endif
 699          enddo
 700          dtmix=0.9*dtmin
 701          nsubmix=dtmicro/dtmix+1
 702          if(nsubmix>100)then
 703             nsubmix_bnd=100
 704          else
 705             nsubmix_bnd=nsubmix
 706          endif
 707          count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
 708          dtmix=dtmicro/nsubmix
 709          fac_srflx = -1.0/(zn(pver)*nsubmix)
 710
 711          do k=1,pver
 712             kp1=min(k+1,pver)
 713             km1=max(k-1,1)
 714             ! maximum overlap assumption
 715             if(cldn(i,kp1).gt.1.e-10_r8)then
 716                overlapp(k)=min(cldn(i,k)/cldn(i,kp1),1._r8)
 717             else
 718                overlapp(k)=1.
 719             endif
 720             if(cldn(i,km1).gt.1.e-10_r8)then
 721                overlapm(k)=min(cldn(i,k)/cldn(i,km1),1._r8)
 722             else
 723                overlapm(k)=1.
 724             endif
 725          enddo
 726
 727
 728 ! rce-comment
 729 !    the activation source(k) = mact(k,m)*raercol(kp1,lmass)
 730 !       should not exceed the rate of transfer of unactivated particles
 731 !       from kp1 to k which = ekkp(k)*raercol(kp1,lmass)
 732 !    however it might if things are not "just right" in subr activate
 733 !    the following is a safety measure to avoid negatives in explmix
 734          do k = 1, pver-1
 735          do m = 1, ntot_amode
 736             nact(k,m) = min( nact(k,m), ekkp(k) )
 737             mact(k,m) = min( mact(k,m), ekkp(k) )
 738          end do
 739          end do
 740
 741
 742 old_cloud_nsubmix_loop: &
 743          do n=1,nsubmix
 744             qncld(:)=qcld(:)
 745 !           switch nsav, nnew so that nsav is the updated aerosol
 746             ntemp=nsav
 747             nsav=nnew
 748             nnew=ntemp
 749             srcn(:)=0.0
 750             do m=1,ntot_amode
 751                lnum=numptr_amode(m)
 752                lnumcw=numptrcw_amode(m)
 753 !              update droplet source
 754 ! rce-comment- activation source in layer k involves particles from k+1
 755 !              srcn(:)=srcn(:)+nact(:,m)*(raercol(:,lnum,nsav))
 756                srcn(1:pver-1)=srcn(1:pver-1)+nact(1:pver-1,m)*(raercol(2:pver,lnum,nsav))
 757 ! rce-comment- new formulation for k=pver
 758 !              srcn(  pver  )=srcn(  pver  )+nact(  pver  ,m)*(raercol(  pver,lnum,nsav))
 759                tmpa = raercol(pver,lnum,nsav)*nact(pver,m) &
 760                     + raercol(pver,lnumcw,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
 761                srcn(pver)=srcn(pver) + max(0.0_r8,tmpa)
 762             enddo
 763             call explmix(qcld,srcn,ekkp,ekkm,overlapp,overlapm,   &
 764                          qncld,surfrate(ixndrop),zero,pver,dtmix,.false.)
 765
 766 ! rce-comment
 767 !    the interstitial particle mixratio is different in clear/cloudy portions
 768 !    of a layer, and generally higher in the clear portion.  (we have/had
 769 !    a method for diagnosing the the clear/cloudy mixratios.)  the activation
 770 !    source terms involve clear air (from below) moving into cloudy air (above).
 771 !    in theory, the clear-portion mixratio should be used when calculating
 772 !    source terms
 773             do m=1,ntot_amode
 774                lnum=numptr_amode(m)
 775                lnumcw=numptrcw_amode(m)
 776                if(lnum>0)then
 777 ! rce-comment -   activation source in layer k involves particles from k+1
 778 !                 source(:)= nact(:,m)*(raercol(:,lnum,nsav))
 779                   source(1:pver-1)= nact(1:pver-1,m)*(raercol(2:pver,lnum,nsav))
 780 ! rce-comment- new formulation for k=pver
 781 !                 source(  pver  )= nact(  pver,  m)*(raercol(  pver,lnum,nsav))
 782                   tmpa = raercol(pver,lnum,nsav)*nact(pver,m) &
 783                        + raercol(pver,lnumcw,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
 784                   source(pver) = max(0.0_r8,tmpa)
 785 !                  flxconv=cflx(i,lnum)*zn(pver)
 786                   flxconv=0.
 787                   call explmix(raercol_cw(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
 788                                raercol_cw(1,lnumcw,nsav),surfrate_cw(lnumcw),zero,pver,dtmix,&
 789                                .false.)
 790                   call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
 791                                raercol(1,lnum,nsav),surfrate(lnum),flxconv,pver,dtmix, &
 792                                .true.,raercol_cw(1,lnumcw,nsav))
 793                endif
 794
 795                do l=1,nspec_amode(m)
 796                   lmass=lmassptr_amode(l,m)
 797                   lmasscw=lmassptrcw_amode(l,m)
 798 ! rce-comment -   activation source in layer k involves particles from k+1
 799 !                 source(:)= mact(:,m)*(raercol(:,lmass,nsav))
 800                   source(1:pver-1)= mact(1:pver-1,m)*(raercol(2:pver,lmass,nsav))
 801 ! rce-comment- new formulation for k=pver
 802 !                 source(  pver  )= mact(  pver  ,m)*(raercol(  pver,lmass,nsav))
 803                   tmpa = raercol(pver,lmass,nsav)*mact(pver,m) &
 804                        + raercol(pver,lmasscw,nsav)*(mact(pver,m) - taumix_internal_pver_inv)
 805                   source(pver) = max(0.0_r8,tmpa)
 806 !                  flxconv=cflx(i,lmass)*zn(pver)
 807                   flxconv=0.
 808
 809                   call explmix(raercol_cw(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
 810                                raercol_cw(1,lmasscw,nsav),surfrate_cw(lmasscw),zero,pver,dtmix,&
 811                                .false.)
 812                   call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
 813                                raercol(1,lmass,nsav),surfrate(lmass),flxconv, pver,dtmix,&
 814                                .true.,raercol_cw(1,lmasscw,nsav))
 815                enddo   ! l
 816 !              lwater=lwaterptr_amode(m)  ! aerosol water
 817 !              source(:)=0.
 818 !              call explmix(   raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm,   &
 819 !                              raercol(1,lwater,nsav),surfrate(lwater),zero,pver,dtmix,&
 820 !                              .true.,source)
 821             enddo   ! m
 822
 823          enddo old_cloud_nsubmix_loop
 824 !          call t_stopf ('nsubmix')
 825
 826 !        evaporate particles again if no cloud
 827
 828          do k=1,pver
 829             if(cldn(i,k).eq.0.)then
 830 !              no cloud
 831                qcld(k)=0.
 832 !              convert activated aerosol to interstitial in decaying cloud
 833                do m=1,ntot_amode
 834                   lnum=numptr_amode(m)
 835                   lnumcw=numptrcw_amode(m)
 836                   if(lnum.gt.0)then
 837                         raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol_cw(k,lnumcw,nnew)
 838                         raercol_cw(k,lnumcw,nnew)=0.
 839                   endif
 840 !
 841                   do l=1,nspec_amode(m)
 842                         lmass=lmassptr_amode(l,m)
 843                         lmasscw=lmassptrcw_amode(l,m)
 844 !                       raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
 845 !                       raercol(k,lmasscw,nnew)=0.
 846                         raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol_cw(k,lmasscw,nnew)
 847                         raercol_cw(k,lmasscw,nnew)=0.
 848                   enddo
 849                 enddo
 850             endif
 851          enddo
 852
 853 !        droplet number
 854
 855          ndropcol(i)=0.
 856          do k=1,pver
 857             ndropmix(i,k)=(qcld(k)-ncldwtr(i,k))*dtinv - nsource(i,k)
 858             tendnd(i,k) = (max(qcld(k),1.e-6_r8)-ncldwtr(i,k))*dtinv
 859             ndropcol(i) = ndropcol(i) + ncldwtr(i,k)*pdel(i,k)
 860          enddo
 861          ndropcol(i) = ndropcol(i)/gravit
 862
 863
 864 !        aerosol tendency
 865          do m=1,ntot_amode
 866             lnum=numptr_amode(m)
 867             lnumcw=numptrcw_amode(m)
 868             if(lnum.gt.0)then
 869                qqcwtend(i,:,lnumcw)=(raercol_cw(:,lnumcw,nnew)-qqcw(lnumcw)%fldcw(i,:))*dtinv
 870                qqcw(lnumcw)%fldcw(i,:)=raercol_cw(:,lnumcw,nnew)    ! update cloud-borne aerosol
 871                raertend(i,:,lnum)= (raercol(:,lnum,nnew)-raer(i,:,lnum))*dtinv
 872             endif
 873             do l=1,nspec_amode(m)
 874                lmass=lmassptr_amode(l,m)
 875                lmasscw=lmassptrcw_amode(l,m)
 876                qqcwtend(i,:,lmasscw)=(raercol_cw(:,lmasscw,nnew)-qqcw(lmasscw)%fldcw(i,:))*dtinv
 877                qqcw(lmasscw)%fldcw(i,:)=raercol_cw(:,lmasscw,nnew)   ! update cloud-borne aerosol
 878                raertend(i,:,lmass)=(raercol(:,lmass,nnew)-raer(i,:,lmass))*dtinv
 879             enddo
 880 !           lwater=lwaterptr_amode(m)
 881 !           raertend(i,:,lwater)=(raercol(:,lwater,nnew)-raer(i,:,lwater))*dtinv
 882          enddo
 883
 884       enddo overall_main_i_loop
 885 ! end of main loop over i/longitude ....................................
 886
 887       call outfld('NDROPCOL', ndropcol  , pcols, lchnk   )
 888       call outfld('NDROPSRC', nsource    , pcols, lchnk   )
 889       call outfld('NDROPMIX', ndropmix    , pcols, lchnk   )
 890       call outfld('LCLOUD  ', cldn    , pcols, lchnk   )
 891       call outfld('WTKE    ', wtke    , pcols, lchnk   )
 892
 893       call ccncalc(lchnk,ncol,temp,cs,raer,qqcw,ccn,psat,supersat,alogsig,npv)
 894       do l=1,psat
 895          call outfld(ccn_name(l), ccn(1,1,l)    , pcols, lchnk   )
 896       enddo
 897
 898 ! do column tendencies
 899       do m = 1, ntot_amode
 900       do lphase = 1, 2
 901       do lspec = 0, nspec_amode(m)+1   ! loop over number + chem constituents + water
 902          if (lspec == 0) then   ! number
 903             if (lphase == 1) then
 904                l = numptr_amode(m)
 905             else
 906                l = numptrcw_amode(m)
 907             endif
 908          else if (lspec <= nspec_amode(m)) then   ! non-water mass
 909             if (lphase == 1) then
 910                l = lmassptr_amode(lspec,m)
 911             else
 912                l = lmassptrcw_amode(lspec,m)
 913             endif
 914          else   ! water mass
 915 !           if (lphase == 1) then
 916 !              l = lwaterptr_amode(m)
 917 !           else
 918                cycle
 919 !           end if
 920          end if
 921          if (l <= 0) cycle
 922
 923 ! mixnuc1 tendency for interstitial = (total tendency) - cflx
 924 ! mixnuc1 tendency for activated    = (total tendency)
 925          coltend(:) = 0.0
 926          if (lphase == 1) then
 927             tmpname = cnst_name(l)
 928             do i = 1, ncol
 929                coltend(i) = sum( pdel(i,:)*raertend(i,:,l) )/gravit
 930 !               coltend(i) = coltend(i) - cflx(i,l)
 931             end do
 932          else
 933             tmpname = cnst_name_cw(l)
 934             do i = 1, ncol
 935                coltend(i) = sum( pdel(i,:)*qqcwtend(i,:,l) )/gravit
 936             end do
 937          end if
 938          fieldname = trim(tmpname) // '_mixnuc1'
 939          call outfld( fieldname, coltend, pcols, lchnk)
 940
 941       end do   ! lspec
 942       end do   ! lphase
 943       end do   ! m
 944
 945       return
 946       end subroutine dropmixnuc
 947
 948    subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
 949                        qold, surfrate, flxconv, pver, dt, is_unact, qactold )
 950
 951 !  explicit integration of droplet/aerosol mixing
 952 !     with source due to activation/nucleation
 953
 954
 955    implicit none
 956    integer, intent(in) :: pver ! number of levels
 957    real(r8), intent(out) :: q(pver) ! mixing ratio to be updated
 958    real(r8), intent(in) :: qold(pver) ! mixing ratio from previous time step
 959    real(r8), intent(in) :: src(pver) ! source due to activation/nucleation (/s)
 960    real(r8), intent(in) :: ekkp(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
 961                       ! below layer k  (k,k+1 interface)
 962    real(r8), intent(in) :: ekkm(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
 963                       ! above layer k  (k,k+1 interface)
 964    real(r8), intent(in) :: overlapp(pver) ! cloud overlap below
 965    real(r8), intent(in) :: overlapm(pver) ! cloud overlap above
 966    real(r8), intent(in) :: surfrate ! surface exchange rate (/s)
 967    real(r8), intent(in) :: flxconv ! convergence of flux from surface
 968    real(r8), intent(in) :: dt ! time step (s)
 969    logical, intent(in) :: is_unact ! true if this is an unactivated species
 970    real(r8), intent(in),optional :: qactold(pver)
 971           ! mixing ratio of ACTIVATED species from previous step
 972           ! *** this should only be present
 973           !     if the current species is unactivated number/sfc/mass
 974
 975    integer k,kp1,km1
 976
 977    if ( is_unact ) then
 978 !     the qactold*(1-overlap) terms are resuspension of activated material
 979       do k=1,pver
 980          kp1=min(k+1,pver)
 981          km1=max(k-1,1)
 982          q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) +       &
 983                            qactold(kp1)*(1.0-overlapp(k)))               &
 984                                   + ekkm(k)*(qold(km1) - qold(k) +     &
 985                            qactold(km1)*(1.0-overlapm(k))) )
 986 !        force to non-negative
 987 !        if(q(k)<-1.e-30)then
 988 !           write(iulog,*)'q=',q(k),' in explmix'
 989             q(k)=max(q(k),0._r8)
 990 !        endif
 991       end do
 992 !     diffusion loss at base of lowest layer
 993       q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
 994 !        force to non-negative
 995 !        if(q(pver)<-1.e-30)then
 996 !           write(iulog,*)'q=',q(pver),' in explmix'
 997             q(pver)=max(q(pver),0._r8)
 998 !        endif
 999    else
1000       do k=1,pver
1001          kp1=min(k+1,pver)
1002          km1=max(k-1,1)
1003          q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) +      &
1004                                     ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
1005 !        force to non-negative
1006 !        if(q(k)<-1.e-30)then
1007 !           write(iulog,*)'q=',q(k),' in explmix'
1008             q(k)=max(q(k),0._r8)
1009 !        endif
1010       end do
1011 !     diffusion loss at base of lowest layer
1012       q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
1013 !        force to non-negative
1014 !        if(q(pver)<-1.e-30)then
1015 !           write(iulog,*)'q=',q(pver),' in explmix'
1016             q(pver)=max(q(pver),0._r8)
1017 !        endif
1018    end if
1019
1020    return
1021    end subroutine explmix
1022
1023    subroutine activate_init
1024 #ifndef WRF_PORT
1025       use ppgrid,  only: pver
1026 #else
1027       use module_cam_support, only: pver
1028 #endif
1029 #ifndef WRF_PORT
1030       use modal_aero_data
1031       use cam_history,   only: addfld, add_default, phys_decomp
1032 #else
1033       use modal_aero_data, only:sigmag_amode => sigmag_amode_mp, &
1034            numptr_amode   => numptr_amode_mp, nspec_amode => nspec_amode_mp, &
1035            lmassptr_amode => lmassptr_amode_mp, dgnum_amode => dgnum_amode_mp
1036       use module_cam_support, only: addfld, add_default, phys_decomp
1037 #endif
1038       use physconst, only: rhoh2o, mwh2o, r_universal
1039       use error_function, only: erf
1040
1041       implicit none
1042
1043       integer l,m
1044       real(r8) arg
1045       integer lnum, lmass
1046
1047       character*16 ccn_longname
1048
1049 !      mathematical constants
1050
1051       zero=0._r8
1052       third=1./3._r8
1053       twothird=2.*third
1054       sixth=1./6._r8
1055       sq2=sqrt(2._r8)
1056       pi=4._r8*atan(1.0_r8)
1057       sqpi=sqrt(pi)
1058
1059       t0=273.
1060       surften=0.076_r8
1061       aten=2.*mwh2o*surften/(r_universal*t0*rhoh2o)
1062       alogaten=log(aten)
1063       alog2=log(2._r8)
1064       alog3=log(3._r8)
1065
1066       do m=1,ntot_amode
1067 !         use only if width of size distribution is prescribed
1068           alogsig(m)=log(sigmag_amode(m))
1069           exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
1070           argfactor(m)=2./(3.*sqrt(2.)*alogsig(m))
1071           f1(m)=0.5*exp(2.5*alogsig(m)*alogsig(m))
1072           f2(m)=1.+0.25*alogsig(m)
1073           lnum=numptr_amode(m)
1074           do l=1,nspec_amode(m)
1075                lmass=lmassptr_amode(l,m)
1076                if(lmass.le.0)then
1077                   write(iulog,*)'lmassptr_amode(',l,m,')=',lmass
1078 #ifdef WRF_PORT
1079                   call wrf_message(iulog)
1080 #endif
1081                   call endrun
1082                endif
1083           enddo
1084           npv(m)=6./(pi*dgnum_amode(m)**3*exp45logsig(m))
1085       enddo
1086
1087       call addfld ('WTKE     ', 'm/s     ', pver, 'A', 'Standard deviation of updraft velocity'                  ,phys_decomp)
1088       call addfld ('LCLOUD   ', '        ', pver, 'A', 'Liquid cloud fraction'                                   ,phys_decomp)
1089       call addfld ('NDROPMIX ', '#/kg/s  ', pver, 'A', 'Droplet number mixing'                     ,phys_decomp)
1090       call addfld ('NDROPSRC ', '#/kg/s  ', pver, 'A', 'Droplet number source'                     ,phys_decomp)
1091       call addfld ('NDROPSNK ', '#/kg/s  ', pver, 'A', 'Droplet number loss by microphysics'       ,phys_decomp)
1092       call addfld ('NDROPCOL ', '#/m2    ', 1,    'A', 'Column droplet number'                     ,phys_decomp)
1093       call add_default ('WTKE    ', 1, ' ')
1094       call add_default ('LCLOUD  ', 1, ' ')
1095
1096       return
1097       end subroutine activate_init
1098
1099       subroutine activate_modal(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
1100                           na, pmode, nmode, volume, sigman, hygro, &
1101                           fn, fm, fluxn, fluxm, flux_fullact )
1102
1103 !      calculates number, surface, and mass fraction of aerosols activated as CCN
1104 !      calculates flux of cloud droplets, surface area, and aerosol mass into cloud
1105 !      assumes an internal mixture within each of up to pmode multiple aerosol modes
1106 !      a gaussiam spectrum of updrafts can be treated.
1107
1108 !      mks units
1109
1110 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1111 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1112
1113       use physconst, only: rair, epsilo, cpair, rh2o, latvap, gravit,   &
1114                                  rhoh2o, mwh2o, r_universal
1115       use wv_saturation, only: estblf, epsqs
1116       use error_function, only: erf
1117
1118       implicit none
1119
1120
1121 !      input
1122
1123       integer pmode ! dimension of modes
1124       real(r8) wbar          ! grid cell mean vertical velocity (m/s)
1125       real(r8) sigw          ! subgrid standard deviation of vertical vel (m/s)
1126       real(r8) wdiab         ! diabatic vertical velocity (0 if adiabatic)
1127       real(r8) wminf         ! minimum updraft velocity for integration (m/s)
1128       real(r8) wmaxf         ! maximum updraft velocity for integration (m/s)
1129       real(r8) tair          ! air temperature (K)
1130       real(r8) rhoair        ! air density (kg/m3)
1131       real(r8) na(pmode)           ! aerosol number concentration (/m3)
1132       integer nmode      ! number of aerosol modes
1133       real(r8) volume(pmode)     ! aerosol volume concentration (m3/m3)
1134       real(r8) sigman(pmode)  ! geometric standard deviation of aerosol size distribution
1135       real(r8) hygro(pmode)  ! hygroscopicity of aerosol mode
1136
1137 !      output
1138
1139       real(r8) fn(pmode)      ! number fraction of aerosols activated
1140       real(r8) fm(pmode)      ! mass fraction of aerosols activated
1141       real(r8) fluxn(pmode)   ! flux of activated aerosol number fraction into cloud (cm/s)
1142       real(r8) fluxm(pmode)   ! flux of activated aerosol mass fraction into cloud (cm/s)
1143       real(r8) flux_fullact   ! flux of activated aerosol fraction assuming 100% activation (cm/s)
1144                               !    rce-comment
1145                               !    used for consistency check -- this should match (ekd(k)*zs(k))
1146                               !    also, fluxm/flux_fullact gives fraction of aerosol mass flux
1147                               !       that is activated
1148
1149 !      local
1150
1151       integer, parameter:: nx=200
1152       integer iquasisect_option, isectional
1153       real(r8) integ,integf
1154       real(r8), parameter :: p0 = 1013.25e2_r8    ! reference pressure (Pa)
1155       real(r8) xmin(ntot_amode),xmax(ntot_amode) ! ln(r) at section interfaces
1156       real(r8) volmin(ntot_amode),volmax(ntot_amode) ! volume at interfaces
1157       real(r8) tmass ! total aerosol mass concentration (g/cm3)
1158       real(r8) sign(ntot_amode)    ! geometric standard deviation of size distribution
1159       real(r8) rm ! number mode radius of aerosol at max supersat (cm)
1160       real(r8) pres ! pressure (Pa)
1161       real(r8) path ! mean free path (m)
1162       real(r8) diff ! diffusivity (m2/s)
1163       real(r8) conduct ! thermal conductivity (Joule/m/sec/deg)
1164       real(r8) diff0,conduct0
1165       real(r8) es ! saturation vapor pressure
1166       real(r8) qs ! water vapor saturation mixing ratio
1167       real(r8) dqsdt ! change in qs with temperature
1168       real(r8) dqsdp ! change in qs with pressure
1169       real(r8) g ! thermodynamic function (m2/s)
1170       real(r8) zeta(ntot_amode), eta(ntot_amode)
1171       real(r8) lnsmax ! ln(smax)
1172       real(r8) alpha
1173       real(r8) gamma
1174       real(r8) beta
1175       real(r8) sqrtg(ntot_amode)
1176       real(r8) :: amcube(ntot_amode) ! cube of dry mode radius (m)
1177       real(r8) :: smcrit(ntot_amode) ! critical supersatuation for activation
1178       real(r8) :: lnsm(ntot_amode) ! ln(smcrit)
1179       real(r8) smc(ntot_amode) ! critical supersaturation for number mode radius
1180       real(r8) sumflx_fullact
1181       real(r8) sumflxn(ntot_amode)
1182       real(r8) sumflxm(ntot_amode)
1183       real(r8) sumfn(ntot_amode)
1184       real(r8) sumfm(ntot_amode)
1185       real(r8) fnold(ntot_amode)   ! number fraction activated
1186       real(r8) fmold(ntot_amode)   ! mass fraction activated
1187       real(r8) wold,gold
1188       real(r8) alogam
1189       real(r8) rlo,rhi,xint1,xint2,xint3,xint4
1190       real(r8) wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
1191       real(r8) dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar
1192       real(r8) alw,sqrtalw
1193       real(r8) smax
1194       real(r8) x,arg
1195       real(r8) xmincoeff,xcut,volcut,surfcut
1196       real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
1197       real(r8) etafactor1,etafactor2(ntot_amode),etafactor2max
1198       integer m,n
1199 !      numerical integration parameters
1200       real(r8), parameter :: eps=0.3_r8,fmax=0.99_r8,sds=3._r8
1201
1202       real(r8), parameter :: namin=1.e6_r8   ! minimum aerosol number concentration (/m3)
1203
1204       integer ndist(nx)  ! accumulates frequency distribution of integration bins required
1205       data ndist/nx*0/
1206       save ndist
1207
1208       if(ntot_amode<pmode)then
1209          write(iulog,*)'ntot_amode,pmode in activate =',ntot_amode,pmode
1210 #ifdef WRF_PORT
1211          call wrf_message(iulog)
1212 #endif
1213          call endrun('activate')
1214       endif
1215
1216       fn(:)=0._r8
1217       fm(:)=0._r8
1218       fluxn(:)=0._r8
1219       fluxm(:)=0._r8
1220       flux_fullact=0._r8
1221
1222       if(nmode.eq.1.and.na(1).lt.1.e-20_r8)return
1223
1224       if(sigw.le.1.e-5_r8.and.wbar.le.0.)return
1225
1226       pres=rair*rhoair*tair
1227       diff0=0.211e-4_r8*(p0/pres)*(tair/t0)**1.94
1228       conduct0=(5.69_r8+0.017_r8*(tair-t0))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
1229       es = estblf(tair)
1230       qs = epsilo*es/(pres-(1.0_r8 - epsqs)*es)
1231       dqsdt=latvap/(rh2o*tair*tair)*qs
1232       alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1./(rair*tair))
1233       gamma=(1+latvap/cpair*dqsdt)/(rhoair*qs)
1234       etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small.
1235
1236       do m=1,nmode
1237          if(volume(m).gt.1.e-39_r8.and.na(m).gt.1.e-39_r8)then
1238 !            number mode radius (m)
1239 !           write(iulog,*)'alogsig,volc,na=',alogsig(m),volc(m),na(m)
1240             amcube(m)=(3.*volume(m)/(4.*pi*exp45logsig(m)*na(m)))  ! only if variable size dist
1241 !           growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
1242 !           should depend on mean radius of mode to account for gas kinetic effects
1243 !           see Fountoukis and Nenes, JGR2005 and Meskhidze et al., JGR2006
1244 !           for approriate size to use for effective diffusivity.
1245             g=1._r8/(rhoh2o/(diff0*rhoair*qs)                                    &
1246               +latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1._r8))
1247             sqrtg(m)=sqrt(g)
1248             beta=2._r8*pi*rhoh2o*g*gamma
1249             etafactor2(m)=1._r8/(na(m)*beta*sqrtg(m))
1250             if(hygro(m).gt.1.e-10)then
1251                smc(m)=2.*aten*sqrt(aten/(27.*hygro(m)*amcube(m))) ! only if variable size dist
1252             else
1253                smc(m)=100.
1254             endif
1255 !           write(iulog,*)'sm,hygro,amcube=',smcrit(m),hygro(m),amcube(m)
1256          else
1257             g=1._r8/(rhoh2o/(diff0*rhoair*qs)                                    &
1258               +latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1._r8))
1259             sqrtg(m)=sqrt(g)
1260             smc(m)=1._r8
1261             etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
1262          endif
1263          lnsm(m)=log(smc(m)) ! only if variable size dist
1264 !        write(iulog,'(a,i4,4g12.2)')'m,na,amcube,hygro,sm,lnsm=', &
1265 !                   m,na(m),amcube(m),hygro(m),sm(m),lnsm(m)
1266       enddo
1267
1268       if(sigw.gt.1.e-5_r8)then ! spectrum of updrafts
1269
1270          wmax=min(wmaxf,wbar+sds*sigw)
1271          wmin=max(wminf,-wdiab)
1272          wmin=max(wmin,wbar-sds*sigw)
1273          w=wmin
1274          dwmax=eps*sigw
1275          dw=dwmax
1276          dfmax=0.2_r8
1277          dfmin=0.1_r8
1278          if(wmax.le.w)then
1279             do m=1,nmode
1280                fluxn(m)=0.
1281                fn(m)=0.
1282                fluxm(m)=0.
1283                fm(m)=0.
1284             enddo
1285             flux_fullact=0._r8
1286             return
1287          endif
1288          do m=1,nmode
1289             sumflxn(m)=0._r8
1290             sumfn(m)=0._r8
1291             fnold(m)=0._r8
1292             sumflxm(m)=0._r8
1293             sumfm(m)=0._r8
1294             fmold(m)=0._r8
1295          enddo
1296          sumflx_fullact=0._r8
1297
1298          fold=0._r8
1299          wold=0._r8
1300          gold=0._r8
1301
1302          dwmin = min( dwmax, 0.01_r8 )
1303
1304          do n=1,200
1305  100        wnuc=w+wdiab
1306 !           write(iulog,*)'wnuc=',wnuc
1307             alw=alpha*wnuc
1308             sqrtalw=sqrt(alw)
1309             etafactor1=alw*sqrtalw
1310
1311               do m=1,nmode
1312                  eta(m)=etafactor1*etafactor2(m)
1313                  zeta(m)=twothird*sqrtalw*aten/sqrtg(m)
1314               enddo
1315
1316               call maxsat(zeta,eta,nmode,smc,smax)
1317 !             write(iulog,*)'w,smax=',w,smax
1318
1319               lnsmax=log(smax)
1320
1321               x=twothird*(lnsm(nmode)-lnsmax)/(sq2*alogsig(nmode))
1322               fnew=0.5_r8*(1._r8-erf(x))
1323
1324
1325             dwnew = dw
1326             if(fnew-fold.gt.dfmax.and.n.gt.1)then
1327 !              reduce updraft increment for greater accuracy in integration
1328                if (dw .gt. 1.01*dwmin) then
1329                   dw=0.7_r8*dw
1330                   dw=max(dw,dwmin)
1331                   w=wold+dw
1332                   go to 100
1333                else
1334                   dwnew = dwmin
1335                endif
1336             endif
1337
1338             if(fnew-fold.lt.dfmin)then
1339 !              increase updraft increment to accelerate integration
1340                dwnew=min(1.5_r8*dw,dwmax)
1341             endif
1342             fold=fnew
1343
1344             z=(w-wbar)/(sigw*sq2)
1345             g=exp(-z*z)
1346             fnmin=1._r8
1347             xmincoeff=alogaten-twothird*(lnsmax-alog2)-alog3
1348
1349             do m=1,nmode
1350 !              modal
1351                x=twothird*(lnsm(m)-lnsmax)/(sq2*alogsig(m))
1352                fn(m)=0.5_r8*(1.-erf(x))
1353                fnmin=min(fn(m),fnmin)
1354 !               integration is second order accurate
1355 !               assumes linear variation of f*g with w
1356                fnbar=(fn(m)*g+fnold(m)*gold)
1357                arg=x-1.5_r8*sq2*alogsig(m)
1358                fm(m)=0.5_r8*(1._r8-erf(arg))
1359                fmbar=(fm(m)*g+fmold(m)*gold)
1360                wb=(w+wold)
1361                if(w.gt.0.)then
1362                   sumflxn(m)=sumflxn(m)+sixth*(wb*fnbar           &
1363                       +(fn(m)*g*w+fnold(m)*gold*wold))*dw
1364                   sumflxm(m)=sumflxm(m)+sixth*(wb*fmbar           &
1365                       +(fm(m)*g*w+fmold(m)*gold*wold))*dw
1366                endif
1367                sumfn(m)=sumfn(m)+0.5_r8*fnbar*dw
1368 !              write(iulog,'(a,9g10.2)')'lnsmax,lnsm(m),x,fn(m),fnold(m),g,gold,fnbar,dw=',lnsmax,lnsm(m),x,fn(m),fnold(m),g,gold,fnbar,dw
1369                fnold(m)=fn(m)
1370                sumfm(m)=sumfm(m)+0.5_r8*fmbar*dw
1371                fmold(m)=fm(m)
1372             enddo
1373 !           same form as sumflxm but replace the fm with 1.0
1374             sumflx_fullact = sumflx_fullact &
1375                            + sixth*(wb*(g+gold) + (g*w+gold*wold))*dw
1376 !            sumg=sumg+0.5_r8*(g+gold)*dw
1377             gold=g
1378             wold=w
1379             dw=dwnew
1380             if(n.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 20
1381             w=w+dw
1382          enddo
1383          write(iulog,*)'do loop is too short in activate'
1384 #ifdef WRF_PORT
1385          call wrf_message(iulog)
1386 #endif
1387          write(iulog,*)'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
1388 #ifdef WRF_PORT
1389          call wrf_message(iulog)
1390 #endif
1391          write(iulog,*)'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
1392 #ifdef WRF_PORT
1393          call wrf_message(iulog)
1394 #endif
1395          write(iulog,*)'wnuc=',wnuc
1396 #ifdef WRF_PORT
1397          call wrf_message(iulog)
1398 #endif
1399          write(iulog,*)'na=',(na(m),m=1,nmode)
1400 #ifdef WRF_PORT
1401          call wrf_message(iulog)
1402 #endif
1403          write(iulog,*)'fn=',(fn(m),m=1,nmode)
1404 #ifdef WRF_PORT
1405          call wrf_message(iulog)
1406 #endif
1407 !   dump all subr parameters to allow testing with standalone code
1408 !   (build a driver that will read input and call activate)
1409          write(iulog,*)'wbar,sigw,wdiab,tair,rhoair,nmode='
1410 #ifdef WRF_PORT
1411          call wrf_message(iulog)
1412 #endif
1413          write(iulog,*) wbar,sigw,wdiab,tair,rhoair,nmode
1414 #ifdef WRF_PORT
1415          call wrf_message(iulog)
1416 #endif
1417          write(iulog,*)'na=',na
1418 #ifdef WRF_PORT
1419          call wrf_message(iulog)
1420 #endif
1421          write(iulog,*)'volume=', (volume(m),m=1,nmode)
1422 #ifdef WRF_PORT
1423          call wrf_message(iulog)
1424 #endif
1425          write(iulog,*)'sigman=',sigman
1426 #ifdef WRF_PORT
1427          call wrf_message(iulog)
1428 #endif
1429          write(iulog,*)'hydro='
1430 #ifdef WRF_PORT
1431          call wrf_message(iulog)
1432 #endif
1433          write(iulog,*) hygro
1434 #ifdef WRF_PORT
1435          call wrf_message(iulog)
1436 #endif
1437
1438          call endrun
1439    20    continue
1440          ndist(n)=ndist(n)+1
1441          if(w.lt.wmaxf)then
1442
1443 !            contribution from all updrafts stronger than wmax
1444 !            assuming constant f (close to fmax)
1445             wnuc=w+wdiab
1446
1447             z1=(w-wbar)/(sigw*sq2)
1448             z2=(wmaxf-wbar)/(sigw*sq2)
1449             g=exp(-z1*z1)
1450             integ=sigw*0.5*sq2*sqpi*(erf(z2)-erf(z1))
1451 !            consider only upward flow into cloud base when estimating flux
1452             wf1=max(w,zero)
1453             zf1=(wf1-wbar)/(sigw*sq2)
1454             gf1=exp(-zf1*zf1)
1455             wf2=max(wmaxf,zero)
1456             zf2=(wf2-wbar)/(sigw*sq2)
1457             gf2=exp(-zf2*zf2)
1458             gf=(gf1-gf2)
1459             integf=wbar*sigw*0.5*sq2*sqpi*(erf(zf2)-erf(zf1))+sigw*sigw*gf
1460
1461             do m=1,nmode
1462                sumflxn(m)=sumflxn(m)+integf*fn(m)
1463                sumfn(m)=sumfn(m)+fn(m)*integ
1464               sumflxm(m)=sumflxm(m)+integf*fm(m)
1465                sumfm(m)=sumfm(m)+fm(m)*integ
1466             enddo
1467 !           same form as sumflxm but replace the fm with 1.0
1468             sumflx_fullact = sumflx_fullact + integf
1469 !            sumg=sumg+integ
1470          endif
1471
1472
1473          do m=1,nmode
1474             fn(m)=sumfn(m)/(sq2*sqpi*sigw)
1475 !            fn(m)=sumfn(m)/(sumg)
1476             if(fn(m).gt.1.01)then
1477                write(iulog,*)'fn=',fn(m),' > 1 in activate'
1478 #ifdef WRF_PORT
1479                call wrf_message(iulog)
1480 #endif
1481                write(iulog,*)'w,m,na,amcube=',w,m,na(m),amcube(m)
1482 #ifdef WRF_PORT
1483                call wrf_message(iulog)
1484 #endif
1485                write(iulog,*)'integ,sumfn,sigw=',integ,sumfn(m),sigw
1486 #ifdef WRF_PORT
1487                call wrf_message(iulog)
1488 #endif
1489                call endrun('activate')
1490             endif
1491             fluxn(m)=sumflxn(m)/(sq2*sqpi*sigw)
1492            fm(m)=sumfm(m)/(sq2*sqpi*sigw)
1493 !            fm(m)=sumfm(m)/(sumg)
1494             if(fm(m).gt.1.01)then
1495                write(iulog,*)'fm=',fm(m),' > 1 in activate'
1496 #ifdef WRF_PORT
1497                call wrf_message(iulog)
1498 #endif
1499             endif
1500             fluxm(m)=sumflxm(m)/(sq2*sqpi*sigw)
1501          enddo
1502 !        same form as fluxm
1503          flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
1504
1505       else
1506
1507 !        single updraft
1508          wnuc=wbar+wdiab
1509
1510          if(wnuc.gt.0._r8)then
1511
1512             w=wbar
1513             alw=alpha*wnuc
1514             sqrtalw=sqrt(alw)
1515             etafactor1=alw*sqrtalw
1516
1517             do m=1,nmode
1518                eta(m)=etafactor1*etafactor2(m)
1519                zeta(m)=twothird*sqrtalw*aten/sqrtg(m)
1520             enddo
1521
1522             call maxsat(zeta,eta,nmode,smc,smax)
1523
1524             lnsmax=log(smax)
1525             xmincoeff=alogaten-twothird*(lnsmax-alog2)-alog3
1526
1527
1528             do m=1,nmode
1529 !                 modal
1530                   x=twothird*(lnsm(m)-lnsmax)/(sq2*alogsig(m))
1531                   fn(m)=0.5_r8*(1._r8-erf(x))
1532                   arg=x-1.5_r8*sq2*alogsig(m)
1533                   fm(m)=0.5_r8*(1._r8-erf(arg))
1534                 if(wbar.gt.0._r8)then
1535                    fluxn(m)=fn(m)*w
1536                    fluxm(m)=fm(m)*w
1537                endif
1538             enddo
1539             flux_fullact = w
1540          endif
1541
1542       endif
1543
1544       return
1545       end subroutine activate_modal
1546
1547       subroutine maxsat(zeta,eta,nmode,smc,smax)
1548
1549 !      calculates maximum supersaturation for multiple
1550 !      competing aerosol modes.
1551
1552 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1553 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1554
1555       implicit none
1556
1557       integer nmode ! number of modes
1558       real(r8) smc(ntot_amode) ! critical supersaturation for number mode radius
1559       real(r8) zeta(ntot_amode), eta(ntot_amode)
1560       real(r8) smax ! maximum supersaturation
1561       integer m  ! mode index
1562       real(r8) sum, g1, g2, g1sqrt, g2sqrt
1563
1564       do m=1,nmode
1565          if(zeta(m).gt.1.e5_r8*eta(m).or.smc(m)*smc(m).gt.1.e5_r8*eta(m))then
1566 !            weak forcing. essentially none activated
1567             smax=1.e-20_r8
1568          else
1569 !            significant activation of this mode. calc activation all modes.
1570             go to 1
1571          endif
1572       enddo
1573
1574       return
1575
1576   1   continue
1577
1578       sum=0
1579       do m=1,nmode
1580          if(eta(m).gt.1.e-20_r8)then
1581             g1=zeta(m)/eta(m)
1582             g1sqrt=sqrt(g1)
1583             g1=g1sqrt*g1
1584             g2=smc(m)/sqrt(eta(m)+3._r8*zeta(m))
1585             g2sqrt=sqrt(g2)
1586             g2=g2sqrt*g2
1587             sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
1588          else
1589             sum=1.e20_r8
1590          endif
1591       enddo
1592
1593       smax=1._r8/sqrt(sum)
1594
1595       return
1596
1597       end subroutine maxsat
1598
1599       subroutine ccncalc(lchnk,ncol,tair,cs,raer,qqcw,ccn,psat,supersat,alogsig,npv)
1600
1601 !     calculates number concentration of aerosols activated as CCN at
1602 !     supersaturation supersat.
1603 !     assumes an internal mixture of a multiple externally-mixed aerosol modes
1604 !     cgs units
1605
1606 !     Ghan et al., Atmos. Res., 1993, 198-221.
1607
1608       use shr_kind_mod,  only: r8 => shr_kind_r8
1609 #ifndef WRF_PORT
1610       use ppgrid,        only: pcols, pver, pverp
1611       use constituents,  only: pcnst
1612       use modal_aero_data
1613 #else
1614       use module_cam_support, only: pcols, pver, pverp, pcnst =>pcnst_mp
1615 #endif
1616       use physconst,     only: rhoh2o, mwh2o, r_universal
1617       use modal_aero_data, only:
1618       use error_function, only: erf
1619
1620       implicit none
1621
1622 !     input
1623
1624       integer lchnk ! chunk index
1625       integer, intent(in) :: ncol ! number of columns
1626       integer, intent(in) :: psat ! number of supesaturations
1627       real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
1628       type(qqcw_type), intent(in) :: QQCW(:)
1629
1630
1631       real(r8), intent(in) :: tair(pcols,pver)          ! air temperature (K)
1632       real(r8), intent(in) :: cs(pcols,pver)    ! air density (kg/m3)
1633       real(r8), intent(in) :: supersat(psat)
1634       real(r8), intent(in) :: npv(ntot_amode) ! number per volume concentration
1635       real(r8), intent(in) :: alogsig(ntot_amode) ! natl log of geometric standard dev of aerosol
1636       real(r8), intent(out) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat (#/m3)
1637
1638 !     local
1639
1640       real(r8) naerosol(pcols) ! interstit+activated aerosol number conc (/m3)
1641       real(r8) vaerosol(pcols) ! interstit+activated aerosol volume conc (m3/m3)
1642       real(r8) exp45logsig(ntot_amode)     ! exp(4.5*alogsig**2)
1643       real(r8) amcube(pcols)
1644       real(r8) super(psat) ! supersaturation
1645       real(r8) amcubecoef(ntot_amode)
1646       real(r8) :: surften       ! surface tension of water w/respect to air (N/m)
1647       real(r8) surften_coef
1648       real(r8) a(pcols) ! surface tension parameter
1649       real(r8) hygro(pcols)  ! aerosol hygroscopicity
1650       real(r8) sm(pcols)  ! critical supersaturation at mode radius
1651       real(r8) argfactor(ntot_amode),arg(pcols)
1652 !     mathematical constants
1653       real(r8) pi
1654       real(r8) twothird,sq2
1655       integer l,m,n,i,k
1656       real(r8) log,cc
1657       real(r8) smcoefcoef,smcoef(pcols)
1658       integer phase ! phase of aerosol
1659
1660       super(:)=supersat(:)*0.01
1661       pi = 4.*atan(1.0)
1662       sq2=sqrt(2._r8)
1663       twothird=2._r8/3._r8
1664       surften=0.076_r8
1665       surften_coef=2._r8*mwh2o*surften/(r_universal*rhoh2o)
1666       smcoefcoef=2._r8/sqrt(27._r8)
1667
1668       do m=1,ntot_amode
1669           exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
1670           amcubecoef(m)=3./(4.*pi*exp45logsig(m))
1671           argfactor(m)=twothird/(sq2*alogsig(m))
1672       end do
1673
1674       do k=1,pver
1675
1676          do i=1,ncol
1677             ccn(i,k,:)=0.
1678             a(i)=surften_coef/tair(i,k)
1679             smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
1680          end do
1681
1682          do m=1,ntot_amode
1683             phase=3 ! interstitial+cloudborne
1684             call loadaer(raer,qqcw,1,ncol,k,m,cs,npv(m), phase, &
1685                          naerosol, vaerosol,  hygro )
1686             where(naerosol(:ncol)>1.e-3)
1687                amcube(:ncol)=amcubecoef(m)*vaerosol(:ncol)/naerosol(:ncol)
1688                sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol)*amcube(:ncol)) ! critical supersaturation
1689             elsewhere
1690                sm(:ncol)=1. ! value shouldn't matter much since naerosol is small
1691             endwhere
1692             do l=1,psat
1693                do i=1,ncol
1694                   arg(i)=argfactor(m)*log(sm(i)/super(l))
1695                   ccn(i,k,l)=ccn(i,k,l)+naerosol(i)*0.5*(1._r8-erf(arg(i)))
1696                enddo
1697              enddo
1698          enddo
1699       enddo
1700       ccn(:ncol,:,:)=ccn(:ncol,:,:)*1.e-6 ! convert from #/m3 to #/cm3
1701
1702       return
1703       end subroutine ccncalc
1704
1705       subroutine loadaer(raer,qqcw,istart,istop,k,m,cs,npv1, phase, &
1706                          naerosol, vaerosol,  hygro )
1707
1708       use shr_kind_mod,  only: r8 => shr_kind_r8
1709 #ifndef WRF_PORT
1710       use abortutils,    only: endrun
1711       use ppgrid,        only: pcols, pver
1712       use modal_aero_data
1713 #else
1714       use module_cam_support, only: pcols, pver, endrun, pcnst =>pcnst_mp
1715       use modal_aero_data,   only: nspec_amode => nspec_amode_mp, &
1716            lmassptr_amode    => lmassptr_amode_mp , lmassptrcw_amode  => lmassptrcw_amode_mp, &
1717            lspectype_amode   => lspectype_amode_mp, specdens_amode    => specdens_amode_mp,   &
1718            spechygro         => spechygro_mp      , numptr_amode      => numptr_amode_mp,     &
1719            numptrcw_amode    => numptrcw_amode_mp , voltonumbhi_amode => voltonumbhi_amode_mp,&
1720            voltonumblo_amode => voltonumblo_amode_mp
1721 #endif
1722
1723
1724       implicit none
1725
1726       !      load aerosol number, volume concentrations, and bulk hygroscopicity
1727
1728       type(qqcw_type), intent(in) :: QQCW(:)
1729
1730       real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
1731        integer, intent(in) :: istart, istop ! start, stop column index
1732        integer, intent(in) ::  m          ! m=mode index
1733        integer, intent(in) ::  k          ! level index
1734        real(r8), intent(in) :: cs(pcols,pver)  ! air density (kg/m3)
1735        integer, intent(in) :: phase ! phase of aerosol: 1 for interstitial, 2 for cloud-borne, 3 for sum
1736        real(r8), intent(in) :: npv1 ! number per volume
1737        real(r8), intent(out) :: naerosol(pcols)                ! interstitial number conc (/m3)
1738        real(r8), intent(out) :: vaerosol(pcols)       ! interstitial+activated volume conc (m3/m3)
1739        real(r8), intent(out) :: hygro(pcols) ! bulk hygroscopicity of mode
1740
1741 !      internal
1742
1743        real(r8) vol(pcols) ! aerosol volume mixing ratio
1744        integer i,lnum,lnumcw,l,ltype,lmass,lmasscw
1745
1746           do i=istart,istop
1747              vaerosol(i)=0._r8
1748              hygro(i)=0._r8
1749           end do
1750
1751           do l=1,nspec_amode(m)
1752              lmass=lmassptr_amode(l,m) ! interstitial
1753              lmasscw=lmassptrcw_amode(l,m) ! cloud-borne
1754              ltype=lspectype_amode(l,m)
1755              if(phase.eq.3)then
1756                 do i=istart,istop
1757 !                  vol(i)=max(raer(i,k,lmass)+raer(i,k,lmasscw),0._r8)/specdens_amode(ltype)
1758                    vol(i)=max(raer(i,k,lmass)+qqcw(lmasscw)%fldcw(i,k),0._r8)/specdens_amode(ltype)
1759                 end do
1760              elseif(phase.eq.2)then
1761                 do i=istart,istop
1762                    vol(i)=max(qqcw(lmasscw)%fldcw(i,k),0._r8)/specdens_amode(ltype)
1763                 end do
1764              elseif(phase.eq.1)then
1765                 do i=istart,istop
1766                    vol(i)=max(raer(i,k,lmass),0._r8)/specdens_amode(ltype)
1767                 end do
1768              else
1769                 write(iulog,*)'phase=',phase,' in loadaer'
1770 #ifdef WRF_PORT
1771                call wrf_message(iulog)
1772 #endif
1773                 call endrun('phase error in loadaer')
1774              endif
1775              do i=istart,istop
1776                 vaerosol(i)=vaerosol(i)+vol(i)
1777                 hygro(i)=hygro(i)+vol(i)*spechygro(ltype)
1778              end do
1779           enddo
1780           do i=istart,istop
1781              if (vaerosol(i) > 1.0e-30_r8) then   ! +++xl add 8/2/2007
1782                hygro(i)=hygro(i)/(vaerosol(i))
1783                vaerosol(i)=vaerosol(i)*cs(i,k)
1784              else
1785                hygro(i)=0.0_r8
1786                vaerosol(i)=0.0_r8
1787              endif
1788           end do
1789
1790           lnum=numptr_amode(m)
1791           lnumcw=numptrcw_amode(m)
1792           if(lnum.gt.0)then
1793 !            aerosol number predicted
1794              if(phase.eq.3)then
1795                 do i=istart,istop
1796 !                  naerosol(i)=(raer(i,k,lnum)+raer(i,k,lnumcw))*cs(i,k)
1797                    naerosol(i)=(raer(i,k,lnum)+qqcw(lnumcw)%fldcw(i,k))*cs(i,k)
1798                 end do
1799              elseif(phase.eq.2)then
1800                 do i=istart,istop
1801                    naerosol(i)=qqcw(lnumcw)%fldcw(i,k)*cs(i,k)
1802                 end do
1803              else
1804                 do i=istart,istop
1805                    naerosol(i)=raer(i,k,lnum)*cs(i,k)
1806                 end do
1807              endif
1808 !            adjust number so that dgnumlo < dgnum < dgnumhi
1809              do i=istart,istop
1810                 naerosol(i) = max( naerosol(i), vaerosol(i)*voltonumbhi_amode(m) )
1811                 naerosol(i) = min( naerosol(i), vaerosol(i)*voltonumblo_amode(m) )
1812              end do
1813           else
1814 !            aerosol number diagnosed from mass and prescribed size
1815              do i=istart,istop
1816                 naerosol(i)=vaerosol(i)*npv1
1817                 naerosol(i)=max(naerosol(i),0._r8)
1818              end do
1819           endif
1820
1821        return
1822        end subroutine loadaer
1823 #endif
1824       end module ndrop
1825
1826
1827