wrfv2_fire/phys/module_mixactivate.F

   1 !************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute,
   3 ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
   4 ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
   5 ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
   6 ! LIABILITY FOR THE USE OF THIS SOFTWARE.
   7 !
   8 ! MOSAIC module: see chem/module_mosaic_driver.F for references and terms
   9 ! of use
  10 !************************************************************************
  11
  12 MODULE module_mixactivate
  13 PRIVATE
  14 PUBLIC prescribe_aerosol_mixactivate, mixactivate
  15 CONTAINS
  16
  17
  18 !----------------------------------------------------------------------
  19 !----------------------------------------------------------------------
  20 ! 06-nov-2005 rce - grid_id & ktau added to arg list
  21 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
  22       subroutine prescribe_aerosol_mixactivate (                      &
  23                 grid_id, ktau, dtstep, naer,                          &
  24                 rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old,       &
  25                 z, dz8w, p_at_w, t_at_w, exch_h,                      &
  26         qv, qc, qi, qndrop3d,                                 &
  27         nsource,                                              &
  28                 ids,ide, jds,jde, kds,kde,                            &
  29                 ims,ime, jms,jme, kms,kme,                            &
  30                 its,ite, jts,jte, kts,kte,                            &
  31                 f_qc, f_qi                                            )
  32
  33 !        USE module_configure
  34
  35 ! wrapper to call mixactivate for mosaic description of aerosol
  36
  37         implicit none
  38
  39 !   subr arguments
  40         integer, intent(in) ::                  &
  41                 grid_id, ktau,                  &
  42                 ids, ide, jds, jde, kds, kde,   &
  43                 ims, ime, jms, jme, kms, kme,   &
  44                 its, ite, jts, jte, kts, kte
  45
  46         real, intent(in) :: dtstep
  47         real, intent(inout) :: naer ! aerosol number (/kg)
  48
  49         real, intent(in),   &
  50                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  51                 rho_phy, th_phy, pi_phy, w,  &
  52                 z, dz8w, p_at_w, t_at_w, exch_h
  53
  54         real, intent(inout),   &
  55                 dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old
  56
  57         real, intent(in),   &
  58                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  59                 qv, qc, qi
  60
  61         real, intent(inout),   &
  62                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  63                 qndrop3d
  64
  65         real, intent(out),   &
  66                 dimension( ims:ime, kms:kme, jms:jme) :: nsource
  67
  68     LOGICAL, OPTIONAL :: f_qc, f_qi
  69
  70 ! local vars
  71         integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem
  72         parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10)
  73         real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity
  74         real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol
  75         real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array
  76         integer i,j,k,l,m,n,p
  77         real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk
  78         integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer
  79         integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
  80           waterptr_aer( maxd_asize, maxd_atype ),   &
  81           numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
  82           ai_phase, cw_phase
  83         real dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
  84              dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
  85              sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
  86              dgnum_aer(maxd_asize, maxd_atype),    & ! median diameter (cm) of number distrib of mode
  87              dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
  88              mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
  89              dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
  90 ! terminology:  (pi/6) * (mean-volume diameter)**3 ==
  91 !       (volume mixing ratio of section/mode)/(number mixing ratio)
  92         real, dimension(ims:ime,kms:kme,jms:jme) :: &
  93              ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat
  94         integer idrydep_onoff
  95         real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy
  96         integer msectional
  97
  98
  99           integer ptr
 100           real maer
 101
 102       if(naer.lt.1.)then
 103              naer=1000.e6 ! #/kg default value
 104       endif
 105           ai_phase=1
 106           cw_phase=2
 107           idrydep_onoff = 0
 108           msectional = 0
 109
 110           t_phy(its:ite,kts:kte,jts:jte)=th_phy(its:ite,kts:kte,jts:jte)*pi_phy(its:ite,kts:kte,jts:jte)
 111
 112       ntype_aer=maxd_atype
 113       do n=1,ntype_aer
 114          nsize_aer(n)=maxd_asize
 115          ncomp_aer(n)=maxd_acomp
 116       end do
 117       nphase_aer=maxd_aphase
 118
 119 ! set properties for each type and size
 120        do n=1,ntype_aer
 121        do m=1,nsize_aer(n)
 122           dlo_sect( m,n )=0.01e-4    ! minimum size of section (cm)
 123           dhi_sect( m,n )=0.5e-4    ! maximum size of section (cm)
 124           sigmag_aer(m,n)=2.      ! geometric standard deviation of aerosol size dist
 125           dgnum_aer(m,n)=0.1e-4       ! median diameter (cm) of number distrib of mode
 126           dpvolmean_aer(m,n) = dgnum_aer(m,n) * exp( 1.5 * (log(sigmag_aer(m,n)))**2 )
 127           end do
 128           do l=1,ncomp_aer(n)
 129              dens_aer( l, n)=1.0   ! density (g/cm3) of material
 130              mw_aer( l, n)=132. ! molecular weight (g/mole)
 131           end do
 132       end do
 133        ptr=0
 134        do p=1,nphase_aer
 135        do n=1,ntype_aer
 136        do m=1,nsize_aer(n)
 137           ptr=ptr+1
 138           numptr_aer( m, n, p )=ptr
 139           if(p.eq.ai_phase)then
 140              chem(its:ite,kts:kte,jts:jte,ptr)=naer
 141           else
 142              chem(its:ite,kts:kte,jts:jte,ptr)=0.
 143           endif
 144         end do ! size
 145         end do ! type
 146         end do ! phase
 147        do p=1,maxd_aphase
 148        do n=1,ntype_aer
 149        do m=1,nsize_aer(n)
 150           do l=1,ncomp_aer(n)
 151           ptr=ptr+1
 152              if(ptr.gt.max_chem)then
 153                 write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
 154                 call wrf_error_fatal("1")
 155              endif
 156              massptr_aer(l, m, n, p)=ptr
 157 ! maer is ug/kg-air;  naer is #/kg-air;  dgnum is cm;  dens_aer is g/cm3
 158 ! 1.e6 factor converts g to ug
 159              maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) *   &
 160                  (dgnum_aer(m,n)**3) * exp( 4.5*((log(sigmag_aer(m,n)))**2) ) )
 161              if(p.eq.ai_phase)then
 162                 chem(its:ite,kts:kte,jts:jte,ptr)=maer
 163              else
 164                 chem(its:ite,kts:kte,jts:jte,ptr)=0.
 165              endif
 166           end do
 167         end do ! size
 168         end do ! type
 169         end do ! phase
 170        do n=1,ntype_aer
 171        do m=1,nsize_aer(n)
 172           ptr=ptr+1
 173           if(ptr.gt.max_chem)then
 174              write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
 175              call wrf_error_fatal("1")
 176           endif
 177 !wig      waterptr_aer(m, n)=ptr
 178           waterptr_aer(m, n)=-1
 179         end do ! size
 180         end do ! type
 181         ddvel(its:ite,jts:jte,:)=0.
 182     hygro(its:ite,kts:kte,jts:jte,:,:) = 0.5
 183
 184 ! 06-nov-2005 rce - grid_id & ktau added to arg list
 185       call mixactivate(  msectional,     &
 186             chem,max_chem,qv,qc,qi,qndrop3d,        &
 187             t_phy, w, ddvel, idrydep_onoff,  &
 188             maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
 189             ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
 190             numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
 191             dens_aer, mw_aer,           &
 192             waterptr_aer, hygro,  ai_phase, cw_phase,                &
 193             ids,ide, jds,jde, kds,kde,                            &
 194             ims,ime, jms,jme, kms,kme,                            &
 195             its,ite, jts,jte, kts,kte,                            &
 196             rho_phy, z, dz8w, p_at_w, t_at_w, exch_h,      &
 197             cldfra, cldfra_old, qsrflx,         &
 198             ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
 199             grid_id, ktau, dtstep, &
 200             F_QC=f_qc, F_QI=f_qi                              )
 201
 202
 203       end subroutine prescribe_aerosol_mixactivate
 204
 205 !----------------------------------------------------------------------
 206 !----------------------------------------------------------------------
 207 !   nov-04 sg ! replaced amode with aer and expanded aerosol dimension to include type and phase
 208
 209 ! 06-nov-2005 rce - grid_id & ktau added to arg list
 210 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
 211 subroutine mixactivate(  msectional,            &
 212            chem, num_chem, qv, qc, qi, qndrop3d,         &
 213            temp, w, ddvel, idrydep_onoff,  &
 214            maxd_acomp, maxd_asize, maxd_atype, maxd_aphase,   &
 215            ncomp_aer, nsize_aer, ntype_aer, nphase_aer,  &
 216            numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dpvolmean_aer,  &
 217            dens_aer, mw_aer,               &
 218            waterptr_aer, hygro, ai_phase, cw_phase,              &
 219            ids,ide, jds,jde, kds,kde,                            &
 220            ims,ime, jms,jme, kms,kme,                            &
 221            its,ite, jts,jte, kts,kte,                            &
 222            rho, zm, dz8w, p_at_w, t_at_w, kvh,      &
 223            cldfra, cldfra_old, qsrflx,          &
 224            ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource,       &
 225            grid_id, ktau, dtstep, &
 226            f_qc, f_qi                       )
 227
 228
 229 !     vertical diffusion and nucleation of cloud droplets
 230 !     assume cloud presence controlled by cloud fraction
 231 !     doesn't distinguish between warm, cold clouds
 232
 233   USE module_model_constants, only: g, rhowater, xlv, cp, rvovrd, r_d, r_v, mwdry, ep_2
 234   USE module_radiation_driver, only: cal_cldfra
 235
 236   implicit none
 237
 238 !     input
 239
 240   INTEGER, intent(in) ::         grid_id, ktau
 241   INTEGER, intent(in) ::         num_chem
 242   integer, intent(in) ::         ids,ide, jds,jde, kds,kde,    &
 243                                  ims,ime, jms,jme, kms,kme,    &
 244                                  its,ite, jts,jte, kts,kte
 245
 246   integer maxd_aphase, nphase_aer, maxd_atype, ntype_aer
 247   integer maxd_asize, maxd_acomp, nsize_aer(maxd_atype)
 248   integer, intent(in) ::   &
 249        ncomp_aer( maxd_atype  ),   &
 250        massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ),   &
 251        waterptr_aer( maxd_asize, maxd_atype ),   &
 252        numptr_aer( maxd_asize, maxd_atype, maxd_aphase), &
 253        ai_phase, cw_phase
 254   integer, intent(in) :: msectional ! 1 for sectional, 0 for modal
 255   integer, intent(in) :: idrydep_onoff
 256   real, intent(in)  ::                       &
 257        dlo_sect( maxd_asize, maxd_atype ),   & ! minimum size of section (cm)
 258        dhi_sect( maxd_asize, maxd_atype ),   & ! maximum size of section (cm)
 259        sigmag_aer(maxd_asize, maxd_atype),   & ! geometric standard deviation of aerosol size dist
 260        dens_aer( maxd_acomp, maxd_atype),    & ! density (g/cm3) of material
 261        mw_aer( maxd_acomp, maxd_atype),      & ! molecular weight (g/mole)
 262        dpvolmean_aer(maxd_asize, maxd_atype)   ! mean-volume diameter (cm) of mode
 263 ! terminology:  (pi/6) * (mean-volume diameter)**3 ==
 264 !       (volume mixing ratio of section/mode)/(number mixing ratio)
 265
 266
 267   REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) :: &
 268        chem ! aerosol molar mixing ratio (ug/kg or #/kg)
 269
 270   REAL, intent(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
 271        qv, qc, qi ! water species (vapor, cloud drops, cloud ice) mixing ratio (g/g)
 272
 273   LOGICAL, OPTIONAL :: f_qc, f_qi
 274
 275   REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
 276        qndrop3d    ! water species mixing ratio (g/g)
 277
 278   real, intent(in) :: dtstep             ! time step for microphysics (s)
 279   real, intent(in) :: temp(ims:ime, kms:kme, jms:jme)    ! temperature (K)
 280   real, intent(in) :: w(ims:ime, kms:kme, jms:jme)   ! vertical velocity (m/s)
 281   real, intent(in) :: rho(ims:ime, kms:kme, jms:jme)    ! density at mid-level  (kg/m3)
 282   REAL, intent(in) :: ddvel( its:ite, jts:jte, num_chem ) ! deposition velocity  (m/s)
 283   real, intent(in) :: zm(ims:ime, kms:kme, jms:jme)     ! geopotential height of level (m)
 284   real, intent(in) :: dz8w(ims:ime, kms:kme, jms:jme) ! layer thickness (m)
 285   real, intent(in) :: p_at_w(ims:ime, kms:kme, jms:jme) ! pressure at layer interface (Pa)
 286   real, intent(in) :: t_at_w(ims:ime, kms:kme, jms:jme) ! temperature at layer interface (K)
 287   real, intent(in) :: kvh(ims:ime, kms:kme, jms:jme)    ! vertical diffusivity (m2/s)
 288   real, intent(inout) :: cldfra_old(ims:ime, kms:kme, jms:jme)! cloud fraction on previous time step
 289   real, intent(inout) :: cldfra(ims:ime, kms:kme, jms:jme)    ! cloud fraction
 290   real, intent(in) :: hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk hygroscopicity   &
 291
 292   REAL, intent(out), DIMENSION( ims:ime, jms:jme, num_chem ) ::   qsrflx ! dry deposition rate for aerosol
 293   real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, &  ! droplet number source (#/kg/s)
 294        ccn1,ccn2,ccn3,ccn4,ccn5,ccn6  ! number conc of aerosols activated at supersat
 295
 296
 297 !--------------------Local storage-------------------------------------
 298 !
 299   real :: dgnum_aer(maxd_asize, maxd_atype) ! median diameter (cm) of number distrib of mode
 300   real :: qndrop(kms:kme)      ! cloud droplet number mixing ratio (#/kg)
 301   real :: lcldfra(kms:kme)     ! liquid cloud fraction
 302   real :: lcldfra_old(kms:kme) ! liquid cloud fraction for previous timestep
 303   real :: wtke(kms:kme)        ! turbulent vertical velocity at base of layer k (m2/s)
 304   real zn(kms:kme)             ! g/pdel (m2/g) for layer
 305   real zs(kms:kme)             ! inverse of distance between levels (m)
 306   real, parameter :: zkmin = 0.01
 307   real, parameter :: zkmax = 100.
 308   real cs(kms:kme)             ! air density (kg/m3) at layer center
 309   real csbot(kms:kme)          ! air density (kg/m3) at layer bottom
 310   real csbot_cscen(kms:kme)    ! csbot(k)/cs(k)
 311   real dz(kms:kme)             ! geometric thickness of layers (m)
 312
 313   real wdiab                   ! diabatic vertical velocity
 314 !      real, parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
 315   real, parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
 316 !      real, parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
 317   real :: qndrop_new(kms:kme)  ! droplet number nucleated on cloud boundaries
 318   real :: ekd(kms:kme)         ! diffusivity for droplets (m2/s)
 319   real :: ekk(kms:kme)         ! density*diffusivity for droplets (kg/m3 m2/s)
 320   real :: srcn(kms:kme)        ! droplet source rate (/s)
 321   real, parameter :: sq2pi = 2.5066282746
 322   real dtinv
 323
 324   integer km1,kp1
 325   real wbar,wmix,wmin,wmax
 326   real dum
 327   real tmpa, tmpb, tmpc, tmpc1, tmpc2, tmpd, tmpe, tmpf
 328   real tmpcourno
 329   real dact
 330   real fluxntot         ! (#/cm2/s)
 331   real fac_srflx
 332   real depvel_drop, depvel_tmp
 333   real, parameter :: depvel_uplimit = 1.0 ! upper limit for dep vels (m/s)
 334   real :: surfrate(num_chem) ! surface exchange rate (/s)
 335   real surfratemax      ! max surfrate for all species treated here
 336   real surfrate_drop    ! surfade exchange rate for droplelts
 337   real dtmin,tinv,dtt
 338   integer nsubmix,nsubmix_bnd
 339   integer i,j,k,m,n,nsub
 340   real dtmix
 341   real alogarg
 342   real qcld
 343   real pi
 344   integer nnew,nsav,ntemp
 345   real :: overlapp(kms:kme),overlapm(kms:kme) ! cloud overlap
 346   real ::  ekkp(kms:kme),ekkm(kms:kme) ! zn*zs*density*diffusivity
 347 !  integer, save :: count_submix(100)=0 ! wig: Note that this is a no-no for tile threads with OMP
 348
 349   integer lnum,lnumcw,l,lmass,lmasscw,lsfc,lsfccw,ltype,lsig,lwater
 350   integer :: ntype(maxd_asize)
 351
 352   real ::  naerosol(maxd_asize, maxd_atype)    ! interstitial aerosol number conc (/m3)
 353   real ::  naerosolcw(maxd_asize, maxd_atype)  ! activated number conc (/m3)
 354   real ::   maerosol(maxd_acomp,maxd_asize, maxd_atype)   ! interstit mass conc (kg/m3)
 355   real ::   maerosolcw(maxd_acomp,maxd_asize, maxd_atype) ! activated mass conc (kg/m3)
 356   real ::   maerosol_tot(maxd_asize, maxd_atype)     ! species-total interstit mass conc (kg/m3)
 357   real ::   maerosol_totcw(maxd_asize, maxd_atype)   ! species-total activated mass conc (kg/m3)
 358   real ::   vaerosol(maxd_asize, maxd_atype) ! interstit+activated aerosol volume conc (m3/m3)
 359   real ::   vaerosolcw(maxd_asize, maxd_atype) ! activated aerosol volume conc (m3/m3)
 360   real ::   raercol(kms:kme,num_chem,2) ! aerosol mass, number mixing ratios
 361   real ::   source(kms:kme) !
 362
 363   real ::   fn(maxd_asize, maxd_atype)         ! activation fraction for aerosol number
 364   real ::   fs(maxd_asize, maxd_atype)         ! activation fraction for aerosol sfcarea
 365   real ::   fm(maxd_asize, maxd_atype)         ! activation fraction for aerosol mass
 366   integer ::   ncomp(maxd_atype)
 367
 368   real ::   fluxn(maxd_asize, maxd_atype)      ! number  activation fraction flux (m/s)
 369   real ::   fluxs(maxd_asize, maxd_atype)      ! sfcarea activation fraction flux (m/s)
 370   real ::   fluxm(maxd_asize, maxd_atype)      ! mass    activation fraction flux (m/s)
 371   real ::   flux_fullact(kms:kme)              ! 100%    activation fraction flux (m/s)
 372 !     note:  activation fraction fluxes are defined as
 373 !     fluxn = [flux of activated aero. number into cloud (#/m2/s)]
 374 !           / [aero. number conc. in updraft, just below cloudbase (#/m3)]
 375
 376   real :: nact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. number  activation rate (/s)
 377   real :: mact(kms:kme,maxd_asize, maxd_atype)  ! fractional aero. mass    activation rate (/s)
 378   real :: npv(maxd_asize, maxd_atype) ! number per volume concentration (/m3)
 379   real scale
 380
 381   real :: hygro_aer(maxd_asize, maxd_atype)  ! hygroscopicity of aerosol mode
 382   real :: exp45logsig     ! exp(4.5*alogsig**2)
 383   real :: alogsig(maxd_asize, maxd_atype) ! natl log of geometric standard dev of aerosol
 384   integer, parameter :: psat=6  ! number of supersaturations to calc ccn concentration
 385   real ccn(kts:kte,psat)        ! number conc of aerosols activated at supersat
 386   real, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
 387        (/0.02,0.05,0.1,0.2,0.5,1.0/)
 388   real super(psat) ! supersaturation
 389   real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m)
 390   real :: ccnfact(psat,maxd_asize, maxd_atype)
 391   real :: amcube(maxd_asize, maxd_atype) ! cube of dry mode radius (m)
 392   real :: argfactor(maxd_asize, maxd_atype)
 393   real aten ! surface tension parameter
 394   real t0 ! reference temperature
 395   real sm ! critical supersaturation
 396   real arg
 397
 398   integer,parameter :: icheck_colmass = 0
 399            ! icheck_colmass > 0 turns on mass/number conservation checking
 400            ! values of 1, 10, 100 produce less to more diagnostics
 401   integer :: colmass_worst_ij( 2, 0:maxd_acomp, maxd_asize, maxd_atype )
 402   integer :: colmass_maxworst_i(3)
 403   real :: colmass_bgn( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
 404   real :: colmass_end( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
 405   real :: colmass_sfc( 0:maxd_acomp, maxd_asize, maxd_atype, maxd_aphase )
 406   real :: colmass_worst( 0:maxd_acomp, maxd_asize, maxd_atype )
 407   real :: colmass_maxworst_r
 408   real :: rhodz( kts:kte ), rhodzsum
 409
 410 !!$#if (defined AIX)
 411 !!$#define ERF erf
 412 !!$#define ERFC erfc
 413 !!$#else
 414 !!$#define ERF erf
 415 !!$    real erf
 416 !!$#define ERFC erfc
 417 !!$    real erfc
 418 !!$#endif
 419
 420   character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
 421
 422
 423   colmass_worst(:,:,:) = 0.0
 424   colmass_worst_ij(:,:,:,:) = -1
 425
 426
 427   arg = 1.0
 428   if (abs(0.8427-ERF_ALT(arg))/0.8427>0.001) then
 429      write (6,*) 'erf_alt(1.0) = ',ERF_ALT(arg)
 430      call wrf_error_fatal('dropmixnuc: Error function error')
 431   endif
 432   arg = 0.0
 433   if (ERF_ALT(arg) /= 0.0) then
 434      write (6,*) 'erf_alt(0.0) = ',ERF_ALT(arg)
 435      call wrf_error_fatal('dropmixnuc: Error function error')
 436   endif
 437
 438   pi = 4.*atan(1.0)
 439   dtinv=1./dtstep
 440
 441   depvel_drop =  0.1 ! prescribed here rather than getting it from dry_dep_driver
 442   if (idrydep_onoff .le. 0) depvel_drop =  0.0
 443   depvel_drop =  min(depvel_drop,depvel_uplimit)
 444
 445   do n=1,ntype_aer
 446      do m=1,nsize_aer(n)
 447         ncomp(n)=ncomp_aer(n)
 448         alogsig(m,n)=alog(sigmag_aer(m,n))
 449         dgnum_aer(m,n) = dpvolmean_aer(m,n) * exp( -1.5*alogsig(m,n)*alogsig(m,n) )
 450 !    print *,'sigmag_aer,dgnum_aer=',sigmag_aer(m,n),dgnum_aer(m,n)
 451         ! npv is used only if number is diagnosed from volume
 452         npv(m,n)=6./(pi*(0.01*dgnum_aer(m,n))**3*exp(4.5*alogsig(m,n)*alogsig(m,n)))
 453      end do
 454   end do
 455   t0=273.15   !wig, 1-Mar-2009: Added .15
 456   aten=2.*surften/(r_v*t0*rhowater)
 457   super(:)=0.01*supersat(:)
 458   do n=1,ntype_aer
 459      do m=1,nsize_aer(n)
 460         exp45logsig=exp(4.5*alogsig(m,n)*alogsig(m,n))
 461         argfactor(m,n)=2./(3.*sqrt(2.)*alogsig(m,n))
 462         amcube(m,n)=3./(4.*pi*exp45logsig*npv(m,n))
 463      enddo
 464   enddo
 465
 466   IF( PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
 467      CALL cal_cldfra(CLDFRA,qc,qi,f_qc,f_qi,      &
 468           ids,ide, jds,jde, kds,kde,              &
 469           ims,ime, jms,jme, kms,kme,              &
 470           its,ite, jts,jte, kts,kte               )
 471   END IF
 472
 473   qsrflx(its:ite,jts:jte,:) = 0.
 474
 475 !     start loop over columns
 476
 477 OVERALL_MAIN_J_LOOP: do j=jts,jte
 478 OVERALL_MAIN_I_LOOP: do i=its,ite
 479
 480 !        load number nucleated into qndrop on cloud boundaries
 481
 482 ! initialization for current i .........................................
 483
 484      do k=kts+1,kte
 485             zs(k)=1./(zm(i,k,j)-zm(i,k-1,j))
 486          enddo
 487          zs(kts)=zs(kts+1)
 488      zs(kte+1)=0.
 489
 490      do k=kts,kte
 491 !!$         if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
 492 !!$!           call wrf_error_fatal("1")
 493 !!$         endif
 494         if(f_qi)then
 495            qcld=qc(i,k,j)+qi(i,k,j)
 496         else
 497            qcld=qc(i,k,j)
 498         endif
 499         if(qcld.lt.-1..or.qcld.gt.1.)then
 500            write(6,'(a,g12.2,a,3i5)')'qcld=',qcld,' for i,k,j=',i,k,j
 501            call wrf_error_fatal("1")
 502         endif
 503         if(qcld.gt.1.e-20)then
 504            lcldfra(k)=cldfra(i,k,j)*qc(i,k,j)/qcld
 505            lcldfra_old(k)=cldfra_old(i,k,j)*qc(i,k,j)/qcld
 506         else
 507            lcldfra(k)=0.
 508            lcldfra_old(k)=0.
 509         endif
 510         qndrop(k)=qndrop3d(i,k,j)
 511 !           qndrop(k)=1.e5
 512         cs(k)=rho(i,k,j) ! air density (kg/m3)
 513         dz(k)=dz8w(i,k,j)
 514         do n=1,ntype_aer
 515            do m=1,nsize_aer(n)
 516               nact(k,m,n)=0.
 517               mact(k,m,n)=0.
 518            enddo
 519         enddo
 520         zn(k)=1./(cs(k)*dz(k))
 521         if(k>kts)then
 522            ekd(k)=kvh(i,k,j)
 523            ekd(k)=max(ekd(k),zkmin)
 524            ekd(k)=min(ekd(k),zkmax)
 525         else
 526            ekd(k)=0
 527         endif
 528 !           diagnose subgrid vertical velocity from diffusivity
 529         if(k.eq.kts)then
 530            wtke(k)=sq2pi*depvel_drop
 531 !               wtke(k)=sq2pi*kvh(i,k,j)
 532 !               wtke(k)=max(wtke(k),wmixmin)
 533         else
 534            wtke(k)=sq2pi*ekd(k)/dz(k)
 535         endif
 536         wtke(k)=max(wtke(k),wmixmin)
 537         nsource(i,k,j)=0.
 538      enddo
 539      nsource(i,kte+1,j) = 0.
 540      qndrop(kte+1)      = 0.
 541      zn(kte+1)          = 0.
 542
 543      do k = kts+1, kte
 544         tmpa = dz(k-1) ; tmpb = dz(k)
 545         tmpc = tmpa/(tmpa + tmpb)
 546         csbot(k) = cs(k-1)*(1.0-tmpc) + cs(k)*tmpc
 547         csbot_cscen(k) = csbot(k)/cs(k)
 548      end do
 549      csbot(kts) = cs(kts)
 550      csbot_cscen(kts) = 1.0
 551      csbot(kte+1) = cs(kte)
 552      csbot_cscen(kte+1) = 1.0
 553
 554     !  calculate surface rate and mass mixing ratio for aerosol
 555
 556      surfratemax = 0.0
 557      nsav=1
 558      nnew=2
 559      surfrate_drop=depvel_drop/dz(kts)
 560      surfratemax = max( surfratemax, surfrate_drop )
 561      do n=1,ntype_aer
 562         do m=1,nsize_aer(n)
 563            lnum=numptr_aer(m,n,ai_phase)
 564            lnumcw=numptr_aer(m,n,cw_phase)
 565            if(lnum>0)then
 566               depvel_tmp = max( 0.0, min( ddvel(i,j,lnum), depvel_uplimit ) )
 567               surfrate(lnum)=depvel_tmp/dz(kts)
 568               surfrate(lnumcw)=surfrate_drop
 569               surfratemax = max( surfratemax, surfrate(lnum) )
 570 !             scale = 1000./mwdry ! moles/kg
 571               scale = 1.
 572               raercol(kts:kte,lnumcw,nsav)=chem(i,kts:kte,j,lnumcw)*scale ! #/kg
 573               raercol(kts:kte,lnum,nsav)=chem(i,kts:kte,j,lnum)*scale
 574            endif
 575            do l=1,ncomp(n)
 576               lmass=massptr_aer(l,m,n,ai_phase)
 577               lmasscw=massptr_aer(l,m,n,cw_phase)
 578 !             scale = mw_aer(l,n)/mwdry
 579               scale = 1.e-9 ! kg/ug
 580               depvel_tmp = max( 0.0, min( ddvel(i,j,lmass), depvel_uplimit ) )
 581               surfrate(lmass)=depvel_tmp/dz(kts)
 582               surfrate(lmasscw)=surfrate_drop
 583               surfratemax = max( surfratemax, surfrate(lmass) )
 584               raercol(kts:kte,lmasscw,nsav)=chem(i,kts:kte,j,lmasscw)*scale ! kg/kg
 585               raercol(kts:kte,lmass,nsav)=chem(i,kts:kte,j,lmass)*scale ! kg/kg
 586            enddo
 587            lwater=waterptr_aer(m,n)
 588            if(lwater>0)then
 589               depvel_tmp = max( 0.0, min( ddvel(i,j,lwater), depvel_uplimit ) )
 590               surfrate(lwater)=depvel_tmp/dz(kts)
 591               surfratemax = max( surfratemax, surfrate(lwater) )
 592               raercol(kts:kte,lwater,nsav)=chem(i,kts:kte,j,lwater) ! don't bother to convert units,
 593              ! because it doesn't contribute to aerosol mass
 594            endif
 595         enddo ! size
 596      enddo ! type
 597
 598
 599 ! mass conservation checking
 600      if (icheck_colmass > 0) then
 601
 602 ! calc initial column burdens
 603         colmass_bgn(:,:,:,:) = 0.0
 604         colmass_end(:,:,:,:) = 0.0
 605         colmass_sfc(:,:,:,:) = 0.0
 606         rhodz(kts:kte) = 1.0/zn(kts:kte)
 607         rhodzsum = sum( rhodz(kts:kte) )
 608         do n=1,ntype_aer
 609            do m=1,nsize_aer(n)
 610               lnum=numptr_aer(m,n,ai_phase)
 611               lnumcw=numptr_aer(m,n,cw_phase)
 612               if(lnum>0)then
 613                  colmass_bgn(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
 614                  colmass_bgn(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
 615               endif
 616               do l=1,ncomp(n)
 617                  lmass=massptr_aer(l,m,n,ai_phase)
 618                  lmasscw=massptr_aer(l,m,n,cw_phase)
 619                  colmass_bgn(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
 620                  colmass_bgn(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
 621               enddo
 622            enddo ! size
 623         enddo ! type
 624      endif ! (icheck_colmass > 0)
 625
 626
 627 !        droplet nucleation/aerosol activation
 628
 629 ! k-loop for growing/shrinking cloud calcs .............................
 630 GROW_SHRINK_MAIN_K_LOOP: do k=kts,kte
 631         km1=max0(k-1,1)
 632         kp1=min0(k+1,kde-1)
 633
 634
 635 !       if(lcldfra(k)-lcldfra_old(k).gt.0.01)then   ! this line is the "old" criterion
 636 !       go to 10
 637
 638 !       growing cloud PLUS
 639 !       upwards vertical advection when lcldfra(k-1) < lcldfra(k)
 640 !
 641 ! tmpc1 = cloud fraction increase from previous time step
 642         tmpc1 = max( (lcldfra(k)-lcldfra_old(k)), 0.0 )
 643         if (k > kts) then
 644 ! tmpc2 = fraction of layer for which vertical advection from below
 645 !             (over dtstep) displaces cloudy air with clear air
 646 !       = (courant number using upwards w at layer bottom)*(difference in cloud fraction)
 647            tmpcourno = dtstep*max(w(i,k,j),0.0)/dz(k)
 648            tmpc2 = max( (lcldfra(k)-lcldfra(km1)), 0.0 ) * tmpcourno
 649            tmpc2 = min( tmpc2, 1.0 )
 650 !          tmpc2 = 0.0   ! this turns off the vertical advect part
 651         else
 652            tmpc2 = 0.0
 653         endif
 654
 655         if ((tmpc1 > 0.001) .or. (tmpc2 > 0.001)) then
 656
 657 !                wmix=wtke(k)
 658            wbar=w(i,k,j)+wtke(k)
 659            wmix=0.
 660            wmin=0.
 661 ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
 662            wmax=50.
 663            wdiab=0
 664
 665 !                load aerosol properties, assuming external mixtures
 666            do n=1,ntype_aer
 667               do m=1,nsize_aer(n)
 668                  call loadaer(raercol(1,1,nsav),k,kms,kme,num_chem,    &
 669                       cs(k), npv(m,n), dlo_sect(m,n),dhi_sect(m,n),             &
 670                       maxd_acomp, ncomp(n), &
 671                       grid_id, ktau, i, j, m, n,   &
 672                       numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
 673                       dens_aer(1,n),    &
 674                       massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
 675                       maerosol(1,m,n), maerosolcw(1,m,n),          &
 676                       maerosol_tot(m,n), maerosol_totcw(m,n),      &
 677                       naerosol(m,n), naerosolcw(m,n),                  &
 678                       vaerosol(m,n), vaerosolcw(m,n) )
 679
 680                  hygro_aer(m,n)=hygro(i,k,j,m,n)
 681               enddo
 682            enddo
 683
 684 ! 06-nov-2005 rce - grid_id & ktau added to arg list
 685            call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
 686                 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
 687                 naerosol, vaerosol,  &
 688                 dlo_sect,dhi_sect,sigmag_aer,hygro_aer,              &
 689                 fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )
 690
 691            do n = 1,ntype_aer
 692               do m = 1,nsize_aer(n)
 693                  lnum   = numptr_aer(m,n,ai_phase)
 694                  lnumcw = numptr_aer(m,n,cw_phase)
 695                  if (tmpc1 > 0.0) then
 696                     dact = tmpc1*fn(m,n)*raercol(k,lnum,nsav) ! interstitial only
 697                  else
 698                     dact = 0.0
 699                  endif
 700                  if (tmpc2 > 0.0) then
 701                     dact = dact + tmpc2*fn(m,n)*raercol(km1,lnum,nsav) ! interstitial only
 702                  endif
 703                  dact = min( dact, 0.99*raercol(k,lnum,nsav) )
 704                  raercol(k,lnumcw,nsav) = raercol(k,lnumcw,nsav)+dact
 705                  raercol(k,lnum,  nsav) = raercol(k,lnum,  nsav)-dact
 706                  qndrop(k) = qndrop(k)+dact
 707                  nsource(i,k,j) = nsource(i,k,j)+dact*dtinv
 708                  do l = 1,ncomp(n)
 709                     lmass   = massptr_aer(l,m,n,ai_phase)
 710                     lmasscw = massptr_aer(l,m,n,cw_phase)
 711                     if (tmpc1 > 0.0) then
 712                        dact = tmpc1*fm(m,n)*raercol(k,lmass,nsav) ! interstitial only
 713                     else
 714                        dact = 0.0
 715                     endif
 716                     if (tmpc2 > 0.0) then
 717                        dact = dact + tmpc2*fm(m,n)*raercol(km1,lmass,nsav) ! interstitial only
 718                     endif
 719                     dact = min( dact, 0.99*raercol(k,lmass,nsav) )
 720                     raercol(k,lmasscw,nsav)  =  raercol(k,lmasscw,nsav)+dact
 721                     raercol(k,lmass,  nsav)  =  raercol(k,lmass,  nsav)-dact
 722                  enddo
 723               enddo
 724            enddo
 725 !   10 continue
 726         endif   ! ((tmpc1 > 0.001) .or. (tmpc2 > 0.001))
 727
 728
 729         if(lcldfra(k) < lcldfra_old(k) .and. lcldfra_old(k) > 1.e-20)then   ! this line is the "old" criterion
 730 !         go to 20
 731
 732 !       shrinking cloud ......................................................
 733
 734 !                droplet loss in decaying cloud
 735            nsource(i,k,j)=nsource(i,k,j)+qndrop(k)*(lcldfra(k)-lcldfra_old(k))*dtinv
 736            qndrop(k)=qndrop(k)*(1.+lcldfra(k)-lcldfra_old(k))
 737 !                 convert activated aerosol to interstitial in decaying cloud
 738
 739            tmpc = (lcldfra(k)-lcldfra_old(k))/lcldfra_old(k)
 740            do n=1,ntype_aer
 741               do m=1,nsize_aer(n)
 742                  lnum=numptr_aer(m,n,ai_phase)
 743                  lnumcw=numptr_aer(m,n,cw_phase)
 744                  if(lnum.gt.0)then
 745                     dact=raercol(k,lnumcw,nsav)*tmpc
 746                     raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
 747                     raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
 748                  endif
 749                  do l=1,ncomp(n)
 750                     lmass=massptr_aer(l,m,n,ai_phase)
 751                     lmasscw=massptr_aer(l,m,n,cw_phase)
 752                     dact=raercol(k,lmasscw,nsav)*tmpc
 753                     raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
 754                     raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
 755                  enddo
 756               enddo
 757            enddo
 758 !             20 continue
 759         endif
 760
 761      enddo GROW_SHRINK_MAIN_K_LOOP
 762 ! end of k-loop for growing/shrinking cloud calcs ......................
 763
 764
 765
 766 ! ......................................................................
 767 ! start of main k-loop for calc of old cloud activation tendencies ..........
 768 ! this loop does "set up" for the nsubmix loop
 769 !
 770 ! rce-comment
 771 !    changed this part of code to use current cloud fraction (lcldfra) exclusively
 772
 773 OLD_CLOUD_MAIN_K_LOOP: do k=kts,kte
 774         km1=max0(k-1,kts)
 775         kp1=min0(k+1,kde-1)
 776         flux_fullact(k) = 0.0
 777         if(lcldfra(k).gt.0.01)then
 778 !          go to 30
 779
 780 !               old cloud
 781               if(lcldfra(k)-lcldfra(km1).gt.0.01.or.k.eq.kts)then
 782
 783 !                   interior cloud
 784 !                   cloud base
 785
 786                  wdiab=0
 787                  wmix=wtke(k) ! spectrum of updrafts
 788                  wbar=w(i,k,j) ! spectrum of updrafts
 789 !                    wmix=0. ! single updraft
 790 !               wbar=wtke(k) ! single updraft
 791 ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
 792                  wmax=50.
 793                  ekd(k)=wtke(k)*dz(k)/sq2pi
 794                  alogarg=max(1.e-20,1/lcldfra(k)-1.)
 795                  wmin=wbar+wmix*0.25*sq2pi*alog(alogarg)
 796
 797                  do n=1,ntype_aer
 798                     do m=1,nsize_aer(n)
 799                        call loadaer(raercol(1,1,nsav),km1,kms,kme,num_chem,    &
 800                             cs(k), npv(m,n),dlo_sect(m,n),dhi_sect(m,n),               &
 801                             maxd_acomp, ncomp(n), &
 802                             grid_id, ktau, i, j, m, n,   &
 803                             numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase),  &
 804                             dens_aer(1,n),   &
 805                             massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase),  &
 806                             maerosol(1,m,n), maerosolcw(1,m,n),          &
 807                             maerosol_tot(m,n), maerosol_totcw(m,n),      &
 808                             naerosol(m,n), naerosolcw(m,n),                  &
 809                             vaerosol(m,n), vaerosolcw(m,n) )
 810
 811                        hygro_aer(m,n)=hygro(i,k,j,m,n)
 812
 813                     enddo
 814                  enddo
 815 !          print *,'old cloud wbar,wmix=',wbar,wmix
 816
 817                  call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
 818                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
 819                       naerosol, vaerosol,  &
 820                       dlo_sect,dhi_sect, sigmag_aer,hygro_aer,                    &
 821                       fn,fs,fm,fluxn,fluxs,fluxm,flux_fullact(k), grid_id, ktau, i, j, k )
 822
 823 ! rce-comment
 824 !    the activation-fraction fluxes (fluxn, fluxm) from subr activate assume that
 825 !       wbar << wmix, which is valid for global-model scale but not mesoscale
 826 !    for wrf-chem application, divide these by flux_fullact to get a
 827 !       "flux-weighted-average" activation fraction, then multiply by (ekd(k)*zs(k))
 828 !       which is the local "turbulent vertical-mixing velocity"
 829                  if (k > kts) then
 830                     if (flux_fullact(k) > 1.0e-20) then
 831                        tmpa = ekd(k)*zs(k)
 832                        tmpf = flux_fullact(k)
 833                        do n=1,ntype_aer
 834                        do m=1,nsize_aer(n)
 835                           tmpb = max( fluxn(m,n), 0.0 ) / max( fluxn(m,n), tmpf )
 836                           fluxn(m,n) = tmpa*tmpb
 837                           tmpb = max( fluxm(m,n), 0.0 ) / max( fluxm(m,n), tmpf )
 838                           fluxm(m,n) = tmpa*tmpb
 839                        enddo
 840                        enddo
 841                     else
 842                        fluxn(:,:) = 0.0
 843                        fluxm(:,:) = 0.0
 844                     endif
 845                  endif
 846
 847                  if(k.gt.kts)then
 848                     tmpc = lcldfra(k)-lcldfra(km1)
 849                  else
 850                     tmpc=lcldfra(k)
 851                  endif
 852 ! rce-comment
 853 !    flux of activated mass into layer k (in kg/m2/s)
 854 !       = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
 855 !    source of activated mass (in kg/kg/s) = flux divergence
 856 !       = actmassflux/(cs(i,k)*dz(i,k))
 857 !    so need factor of csbot_cscen = csbot(k)/cs(i,k)
 858 !                tmpe=1./(dz(k))
 859                  tmpe = csbot_cscen(k)/(dz(k))
 860                  fluxntot=0.
 861                  do n=1,ntype_aer
 862                  do m=1,nsize_aer(n)
 863                     fluxn(m,n)=fluxn(m,n)*tmpc
 864 !                   fluxs(m,n)=fluxs(m,n)*tmpc
 865                     fluxm(m,n)=fluxm(m,n)*tmpc
 866                     lnum=numptr_aer(m,n,ai_phase)
 867                     fluxntot=fluxntot+fluxn(m,n)*raercol(km1,lnum,nsav)
 868 !             print *,'fn=',fn(m,n),' for m,n=',m,n
 869 !             print *,'old cloud tmpc=',tmpc,' fn=',fn(m,n),' for m,n=',m,n
 870                     nact(k,m,n)=nact(k,m,n)+fluxn(m,n)*tmpe
 871                     mact(k,m,n)=mact(k,m,n)+fluxm(m,n)*tmpe
 872                  enddo
 873                  enddo
 874                  flux_fullact(k) = flux_fullact(k)*tmpe
 875                  nsource(i,k,j)=nsource(i,k,j)+fluxntot*zs(k)
 876                  fluxntot=fluxntot*cs(k)
 877               endif
 878 !       30 continue
 879
 880         else
 881 !       go to 40
 882 !              no cloud
 883            if(qndrop(k).gt.10000.e6)then
 884               print *,'i,k,j,lcldfra,qndrop=',i,k,j,lcldfra(k),qndrop(k)
 885               print *,'cldfra,ql,qi',cldfra(i,k,j),qc(i,k,j),qi(i,k,j)
 886            endif
 887            nsource(i,k,j)=nsource(i,k,j)-qndrop(k)*dtinv
 888            qndrop(k)=0.
 889 !              convert activated aerosol to interstitial in decaying cloud
 890            do n=1,ntype_aer
 891               do m=1,nsize_aer(n)
 892                  lnum=numptr_aer(m,n,ai_phase)
 893                  lnumcw=numptr_aer(m,n,cw_phase)
 894                  if(lnum.gt.0)then
 895                     raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol(k,lnumcw,nsav)
 896                     raercol(k,lnumcw,nsav)=0.
 897                  endif
 898                  do l=1,ncomp(n)
 899                     lmass=massptr_aer(l,m,n,ai_phase)
 900                     lmasscw=massptr_aer(l,m,n,cw_phase)
 901                     raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol(k,lmasscw,nsav)
 902                     raercol(k,lmasscw,nsav)=0.
 903                  enddo
 904               enddo
 905            enddo
 906 !      40 continue
 907         endif
 908
 909      enddo OLD_CLOUD_MAIN_K_LOOP
 910
 911 !    cycle OVERALL_MAIN_I_LOOP
 912
 913
 914 !    switch nsav, nnew so that nnew is the updated aerosol
 915
 916      ntemp=nsav
 917      nsav=nnew
 918      nnew=ntemp
 919
 920 !    load new droplets in layers above, below clouds
 921
 922      dtmin=dtstep
 923      ekk(kts)=0.0
 924 ! rce-comment -- ekd(k) is eddy-diffusivity at k/k-1 interface
 925 !   want ekk(k) = ekd(k) * (density at k/k-1 interface)
 926      do k=kts+1,kte
 927         ekk(k)=ekd(k)*csbot(k)
 928      enddo
 929      ekk(kte+1)=0.0
 930      do k=kts,kte
 931         ekkp(k)=zn(k)*ekk(k+1)*zs(k+1)
 932         ekkm(k)=zn(k)*ekk(k)*zs(k)
 933         tinv=ekkp(k)+ekkm(k)
 934         if(k.eq.kts)tinv=tinv+surfratemax
 935         if(tinv.gt.1.e-6)then
 936            dtt=1./tinv
 937            dtmin=min(dtmin,dtt)
 938         endif
 939      enddo
 940      dtmix=0.9*dtmin
 941      nsubmix=dtstep/dtmix+1
 942      if(nsubmix>100)then
 943         nsubmix_bnd=100
 944      else
 945         nsubmix_bnd=nsubmix
 946      endif
 947 !     count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
 948      dtmix=dtstep/nsubmix
 949      fac_srflx = -1.0/(zn(1)*nsubmix)
 950
 951      do k=kts,kte
 952         kp1=min(k+1,kde-1)
 953         km1=max(k-1,1)
 954         if(lcldfra(kp1).gt.0)then
 955            overlapp(k)=min(lcldfra(k)/lcldfra(kp1),1.)
 956         else
 957            overlapp(k)=1.
 958         endif
 959         if(lcldfra(km1).gt.0)then
 960            overlapm(k)=min(lcldfra(k)/lcldfra(km1),1.)
 961         else
 962            overlapm(k)=1.
 963         endif
 964      enddo
 965
 966
 967
 968 ! ......................................................................
 969 ! start of nsubmix-loop for calc of old cloud activation tendencies ....
 970 OLD_CLOUD_NSUBMIX_LOOP: do nsub=1,nsubmix
 971         qndrop_new(kts:kte)=qndrop(kts:kte)
 972 !           switch nsav, nnew so that nsav is the updated aerosol
 973         ntemp=nsav
 974         nsav=nnew
 975         nnew=ntemp
 976         srcn(:)=0.0
 977         do n=1,ntype_aer
 978            do m=1,nsize_aer(n)
 979               lnum=numptr_aer(m,n,ai_phase)
 980 !              update droplet source
 981 ! rce-comment - activation source in layer k involves particles from k-1
 982 !             srcn(kts  :kte)=srcn(kts  :kte)+nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
 983               srcn(kts+1:kte)=srcn(kts+1:kte)+nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
 984 ! rce-comment - new formulation for k=kts should be implemented
 985               srcn(kts      )=srcn(kts      )+nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
 986            enddo
 987         enddo
 988         call explmix(qndrop,srcn,ekkp,ekkm,overlapp,overlapm,   &
 989              qndrop_new,surfrate_drop,kms,kme,kts,kte,dtmix,.false.)
 990         do n=1,ntype_aer
 991            do m=1,nsize_aer(n)
 992               lnum=numptr_aer(m,n,ai_phase)
 993               lnumcw=numptr_aer(m,n,cw_phase)
 994               if(lnum>0)then
 995 ! rce-comment - activation source in layer k involves particles from k-1
 996 !                source(kts  :kte)= nact(kts  :kte,m,n)*(raercol(kts:kte  ,lnum,nsav))
 997                  source(kts+1:kte)= nact(kts+1:kte,m,n)*(raercol(kts:kte-1,lnum,nsav))
 998 ! rce-comment - new formulation for k=kts should be implemented
 999                  source(kts      )= nact(kts      ,m,n)*(raercol(kts      ,lnum,nsav))
1000                  call explmix(raercol(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
1001                       raercol(1,lnumcw,nsav),surfrate(lnumcw),kms,kme,kts,kte,dtmix,&
1002                       .false.)
1003                  call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
1004                       raercol(1,lnum,nsav),surfrate(lnum),kms,kme,kts,kte,dtmix, &
1005                       .true.,raercol(1,lnumcw,nsav))
1006                  qsrflx(i,j,lnum) = qsrflx(i,j,lnum) + fac_srflx*            &
1007                       raercol(kts,lnum,nsav)*surfrate(lnum)
1008                  qsrflx(i,j,lnumcw) = qsrflx(i,j,lnumcw) + fac_srflx*        &
1009                       raercol(kts,lnumcw,nsav)*surfrate(lnumcw)
1010                  if (icheck_colmass > 0) then
1011                     tmpf = dtmix*rhodz(kts)
1012                     colmass_sfc(0,m,n,1) = colmass_sfc(0,m,n,1) &
1013                           + raercol(kts,lnum  ,nsav)*surfrate(lnum  )*tmpf
1014                     colmass_sfc(0,m,n,2) = colmass_sfc(0,m,n,2) &
1015                           + raercol(kts,lnumcw,nsav)*surfrate(lnumcw)*tmpf
1016                  endif
1017               endif
1018               do l=1,ncomp(n)
1019                  lmass=massptr_aer(l,m,n,ai_phase)
1020                  lmasscw=massptr_aer(l,m,n,cw_phase)
1021 ! rce-comment - activation source in layer k involves particles from k-1
1022 !                source(kts  :kte)= mact(kts  :kte,m,n)*(raercol(kts:kte  ,lmass,nsav))
1023                  source(kts+1:kte)= mact(kts+1:kte,m,n)*(raercol(kts:kte-1,lmass,nsav))
1024 ! rce-comment - new formulation for k=kts should be implemented
1025                  source(kts      )= mact(kts      ,m,n)*(raercol(kts      ,lmass,nsav))
1026                  call explmix(raercol(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
1027                       raercol(1,lmasscw,nsav),surfrate(lmasscw),kms,kme,kts,kte,dtmix,  &
1028                       .false.)
1029                  call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm,  &
1030                       raercol(1,lmass,nsav),surfrate(lmass),kms,kme,kts,kte,dtmix,  &
1031                       .true.,raercol(1,lmasscw,nsav))
1032                  qsrflx(i,j,lmass) = qsrflx(i,j,lmass) + fac_srflx*          &
1033                       raercol(kts,lmass,nsav)*surfrate(lmass)
1034                  qsrflx(i,j,lmasscw) = qsrflx(i,j,lmasscw) + fac_srflx*      &
1035                       raercol(kts,lmasscw,nsav)*surfrate(lmasscw)
1036                  if (icheck_colmass > 0) then
1037                     ! colmass_sfc calculation
1038                     !    colmass_bgn/end = bgn/end column burden = sum.over.k.of{ rho(k)*dz(k)*chem(k,l) }
1039                     !    colmass_sfc = surface loss over dtstep
1040                     !       = sum.over.nsubmix.substeps{ depvel(l)*rho(kts)*chem(kts,l)*dtmix }
1041                     !    surfrate(l) = depvel(l)/dz(kts) so need to multiply by dz(kts)
1042                     !    for mass, raercol(k,l) = chem(k,l)*1.0e-9, so need to multiply by 1.0e9
1043                     tmpf = dtmix*rhodz(kts)*1.0e9
1044                     colmass_sfc(l,m,n,1) = colmass_sfc(l,m,n,1) &
1045                           + raercol(kts,lmass  ,nsav)*surfrate(lmass  )*tmpf
1046                     colmass_sfc(l,m,n,2) = colmass_sfc(l,m,n,2) &
1047                           + raercol(kts,lmasscw,nsav)*surfrate(lmasscw)*tmpf
1048                  endif
1049               enddo
1050               lwater=waterptr_aer(m,n)  ! aerosol water
1051               if(lwater>0)then
1052                  source(:)=0.
1053                  call explmix(   raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm,   &
1054                       raercol(1,lwater,nsav),surfrate(lwater),kms,kme,kts,kte,dtmix,  &
1055                       .true.,source)
1056               endif
1057            enddo ! size
1058         enddo ! type
1059
1060      enddo OLD_CLOUD_NSUBMIX_LOOP
1061
1062 !    cycle OVERALL_MAIN_I_LOOP
1063
1064 !        evaporate particles again if no cloud
1065
1066      do k=kts,kte
1067         if(lcldfra(k).eq.0.)then
1068
1069 !              no cloud
1070
1071            qndrop(k)=0.
1072 !              convert activated aerosol to interstitial in decaying cloud
1073            do n=1,ntype_aer
1074               do m=1,nsize_aer(n)
1075                  lnum=numptr_aer(m,n,ai_phase)
1076                  lnumcw=numptr_aer(m,n,cw_phase)
1077                  if(lnum.gt.0)then
1078                     raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew)
1079                     raercol(k,lnumcw,nnew)=0.
1080                  endif
1081                  do l=1,ncomp(n)
1082                     lmass=massptr_aer(l,m,n,ai_phase)
1083                     lmasscw=massptr_aer(l,m,n,cw_phase)
1084                     raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
1085                     raercol(k,lmasscw,nnew)=0.
1086                  enddo
1087               enddo
1088            enddo
1089         endif
1090      enddo
1091
1092 !         cycle OVERALL_MAIN_I_LOOP
1093
1094 !        droplet number
1095
1096      do k=kts,kte
1097 !       if(lcldfra(k).gt.0.1)then
1098 !           write(6,'(a,3i5,f12.1)')'i,j,k,qndrop=',i,j,k,qndrop(k)
1099 !       endif
1100         if(qndrop(k).lt.-10.e6.or.qndrop(k).gt.1.e12)then
1101            write(6,'(a,g12.2,a,3i5)')'after qndrop=',qndrop(k),' for i,k,j=',i,k,j
1102         endif
1103
1104         qndrop3d(i,k,j) = max(qndrop(k),1.e-6)
1105
1106         if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
1107            write(6,'(a,g12.2,a,3i5)')'after qndrop3d=',qndrop3d(i,k,j),' for i,k,j=',i,k,j
1108         endif
1109         if(qc(i,k,j).lt.-1..or.qc(i,k,j).gt.1.)then
1110            write(6,'(a,g12.2,a,3i5)')'qc=',qc(i,k,j),' for i,k,j=',i,k,j
1111            call wrf_error_fatal("1")
1112         endif
1113         if(qi(i,k,j).lt.-1..or.qi(i,k,j).gt.1.)then
1114            write(6,'(a,g12.2,a,3i5)')'qi=',qi(i,k,j),' for i,k,j=',i,k,j
1115            call wrf_error_fatal("1")
1116         endif
1117         if(qv(i,k,j).lt.-1..or.qv(i,k,j).gt.1.)then
1118            write(6,'(a,g12.2,a,3i5)')'qv=',qv(i,k,j),' for i,k,j=',i,k,j
1119            call wrf_error_fatal("1")
1120         endif
1121         cldfra_old(i,k,j) = cldfra(i,k,j)
1122 !       if(k.gt.6.and.k.lt.11)cldfra_old(i,k,j)=1.
1123      enddo
1124
1125 !    cycle OVERALL_MAIN_I_LOOP
1126
1127 !        update chem and convert back to mole/mole
1128
1129      ccn(:,:) = 0.
1130      do n=1,ntype_aer
1131         do m=1,nsize_aer(n)
1132            lnum=numptr_aer(m,n,ai_phase)
1133            lnumcw=numptr_aer(m,n,cw_phase)
1134            if(lnum.gt.0)then
1135               !          scale=mwdry*0.001
1136               scale = 1.
1137               chem(i,kts:kte,j,lnumcw)= raercol(kts:kte,lnumcw,nnew)*scale
1138               chem(i,kts:kte,j,lnum)= raercol(kts:kte,lnum,nnew)*scale
1139            endif
1140            do l=1,ncomp(n)
1141               lmass=massptr_aer(l,m,n,ai_phase)
1142               lmasscw=massptr_aer(l,m,n,cw_phase)
1143 !          scale = mwdry/mw_aer(l,n)
1144               scale = 1.e9
1145               chem(i,kts:kte,j,lmasscw)=raercol(kts:kte,lmasscw,nnew)*scale ! ug/kg
1146               chem(i,kts:kte,j,lmass)=raercol(kts:kte,lmass,nnew)*scale ! ug/kg
1147            enddo
1148            lwater=waterptr_aer(m,n)
1149            if(lwater>0)chem(i,kts:kte,j,lwater)=raercol(kts:kte,lwater,nnew) ! don't convert units
1150            do k=kts,kte
1151               sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n)))
1152               do l=1,psat
1153                  arg=argfactor(m,n)*log(sm/super(l))
1154                  if(arg<2)then
1155                     if(arg<-2)then
1156                        ccnfact(l,m,n)=1.e-6 ! convert from #/m3 to #/cm3
1157                     else
1158                        ccnfact(l,m,n)=1.e-6*0.5*ERFC_NUM_RECIPES(arg)
1159                     endif
1160                  else
1161                     ccnfact(l,m,n) = 0.
1162                  endif
1163 !                 ccn concentration as diagnostic
1164 !                 assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles
1165                  ccn(k,l)=ccn(k,l)+(raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew))*cs(k)*ccnfact(l,m,n)
1166               enddo
1167            enddo
1168         enddo
1169      enddo
1170      do l=1,psat
1171         !wig, 22-Nov-2006: added vertical bounds to prevent out-of-bounds at top
1172         if(l.eq.1)ccn1(i,kts:kte,j)=ccn(:,l)
1173         if(l.eq.2)ccn2(i,kts:kte,j)=ccn(:,l)
1174         if(l.eq.3)ccn3(i,kts:kte,j)=ccn(:,l)
1175         if(l.eq.4)ccn4(i,kts:kte,j)=ccn(:,l)
1176         if(l.eq.5)ccn5(i,kts:kte,j)=ccn(:,l)
1177         if(l.eq.6)ccn6(i,kts:kte,j)=ccn(:,l)
1178      end do
1179
1180 ! mass conservation checking
1181      if (icheck_colmass > 0) then
1182 ! calc final column burdens
1183         do n=1,ntype_aer
1184         do m=1,nsize_aer(n)
1185            lnum=numptr_aer(m,n,ai_phase)
1186            lnumcw=numptr_aer(m,n,cw_phase)
1187            if(lnum>0)then
1188               colmass_end(0,m,n,1) = sum( chem(i,kts:kte,j,lnum  )*rhodz(kts:kte) )
1189               colmass_end(0,m,n,2) = sum( chem(i,kts:kte,j,lnumcw)*rhodz(kts:kte) )
1190            endif
1191            do l=1,ncomp(n)
1192               lmass=massptr_aer(l,m,n,ai_phase)
1193               lmasscw=massptr_aer(l,m,n,cw_phase)
1194               colmass_end(l,m,n,1) = sum( chem(i,kts:kte,j,lmass  )*rhodz(kts:kte) )
1195               colmass_end(l,m,n,2) = sum( chem(i,kts:kte,j,lmasscw)*rhodz(kts:kte) )
1196            enddo
1197         enddo ! size
1198         enddo ! type
1199 ! calc burden change errors for each interstitial/activated pair
1200         do n=1,ntype_aer
1201         do m=1,nsize_aer(n)
1202            do l=0,ncomp(n)
1203               ! tmpa & tmpb = beginning & ending column burden divided by rhodzsum,
1204               !             = beginning & ending column-mean mixing ratios
1205               ! tmpc = loss to surface divided by rhodzsum,
1206               tmpa = ( colmass_bgn(l,m,n,1) + colmass_bgn(l,m,n,2) )/rhodzsum
1207               tmpb = ( colmass_end(l,m,n,1) + colmass_end(l,m,n,2) )/rhodzsum
1208               tmpc = ( colmass_sfc(l,m,n,1) + colmass_sfc(l,m,n,2) )/rhodzsum
1209
1210               ! tmpd = ((final burden) + (sfc loss)) - (initial burden)
1211               !      = burden change error
1212               tmpd = (tmpb + tmpc) - tmpa
1213               tmpe = max( tmpa, 1.0e-20 )
1214
1215               ! tmpf = (burden change error) / (initial burden)
1216               if (abs(tmpd) < 1.0e5*tmpe) then
1217                  tmpf = tmpd/tmpe
1218               else if (tmpf < 0.0) then
1219                  tmpf = -1.0e5
1220               else
1221                  tmpf = 1.0e5
1222               end if
1223               if (abs(tmpf) > abs(colmass_worst(l,m,n))) then
1224                  colmass_worst(l,m,n) = tmpf
1225                  colmass_worst_ij(1,l,m,n) = i
1226                  colmass_worst_ij(2,l,m,n) = j
1227               endif
1228            enddo
1229         enddo ! size
1230         enddo ! type
1231      endif ! (icheck_colmass > 0)
1232
1233
1234      enddo OVERALL_MAIN_I_LOOP ! end of main loop over i
1235      enddo OVERALL_MAIN_J_LOOP ! end of main loop over j
1236
1237
1238 ! mass conservation checking
1239      if (icheck_colmass > 0) then
1240         if (icheck_colmass >= 100) write(*,'(a)') &
1241              'mixactivate colmass worst errors bgn - type, size, comp, err, i, j'
1242         colmass_maxworst_r = 0.0
1243         colmass_maxworst_i(:) = -1
1244         do n=1,ntype_aer
1245         do m=1,nsize_aer(n)
1246            do l=0,ncomp(n)
1247               if (icheck_colmass >= 100) &
1248                  write(*,'(3i3,1p,e10.2,2i4)') n, m, l, &
1249                  colmass_worst(l,m,n), colmass_worst_ij(1:2,l,m,n)
1250               if (abs(colmass_worst(l,m,n)) > abs(colmass_maxworst_r)) then
1251                  colmass_maxworst_r = colmass_worst(l,m,n)
1252                  colmass_maxworst_i(1) = n
1253                  colmass_maxworst_i(2) = m
1254                  colmass_maxworst_i(3) = l
1255               end if
1256            enddo
1257         enddo ! size
1258         enddo ! type
1259         if ((icheck_colmass >= 10) .or. (abs(colmass_maxworst_r) >= 1.0e-6)) &
1260              write(*,'(a,3i3,1p,e10.2)') 'mixactivate colmass maxworst', &
1261              colmass_maxworst_i(1:3), colmass_maxworst_r
1262      endif ! (icheck_colmass > 0)
1263
1264
1265      return
1266    end subroutine mixactivate
1267
1268
1269 !----------------------------------------------------------------------
1270 !----------------------------------------------------------------------
1271    subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
1272                        qold, surfrate, kms, kme, kts, kte, dt, &
1273                        is_unact, qactold )
1274
1275 !  explicit integration of droplet/aerosol mixing
1276 !     with source due to activation/nucleation
1277
1278
1279    implicit none
1280    integer, intent(in) :: kms,kme ! number of levels for array definition
1281    integer, intent(in) :: kts,kte ! number of levels for looping
1282    real, intent(inout) :: q(kms:kme) ! mixing ratio to be updated
1283    real, intent(in) :: qold(kms:kme) ! mixing ratio from previous time step
1284    real, intent(in) :: src(kms:kme) ! source due to activation/nucleation (/s)
1285    real, intent(in) :: ekkp(kms:kme) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1286                       ! below layer k  (k,k+1 interface)
1287    real, intent(in) :: ekkm(kms:kme) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1288                       ! above layer k  (k,k+1 interface)
1289    real, intent(in) :: overlapp(kms:kme) ! cloud overlap below
1290    real, intent(in) :: overlapm(kms:kme) ! cloud overlap above
1291    real, intent(in) :: surfrate ! surface exchange rate (/s)
1292    real, intent(in) :: dt ! time step (s)
1293    logical, intent(in) :: is_unact ! true if this is an unactivated species
1294    real, intent(in),optional :: qactold(kms:kme)
1295           ! mixing ratio of ACTIVATED species from previous step
1296           ! *** this should only be present
1297           !     if the current species is unactivated number/sfc/mass
1298
1299    integer k,kp1,km1
1300
1301    if ( is_unact ) then
1302 !     the qactold*(1-overlap) terms are resuspension of activated material
1303       do k=kts,kte
1304          kp1=min(k+1,kte)
1305          km1=max(k-1,kts)
1306          q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) +  &
1307                            qactold(kp1)*(1.0-overlapp(k)))               &
1308                                   + ekkm(k)*(qold(km1) - qold(k) +     &
1309                            qactold(km1)*(1.0-overlapm(k))) )
1310 !         if(q(k)<-1.e-30)then ! force to non-negative
1311 !            print *,'q=',q(k),' in explmix'
1312              q(k)=max(q(k),0.)
1313 !         endif
1314       end do
1315
1316    else
1317       do k=kts,kte
1318          kp1=min(k+1,kte)
1319          km1=max(k-1,kts)
1320          q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) +  &
1321                                     ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
1322 !        if(q(k)<-1.e-30)then ! force to non-negative
1323 !           print *,'q=',q(k),' in explmix'
1324             q(k)=max(q(k),0.)
1325 !        endif
1326       end do
1327    end if
1328
1329 !  dry deposition loss at base of lowest layer
1330    q(kts)=q(kts)-surfrate*qold(kts)*dt
1331 !  if(q(kts)<-1.e-30)then ! force to non-negative
1332 !     print *,'q=',q(kts),' in explmix'
1333       q(kts)=max(q(kts),0.)
1334 !  endif
1335
1336    return
1337    end subroutine explmix
1338
1339 !----------------------------------------------------------------------
1340 !----------------------------------------------------------------------
1341 ! 06-nov-2005 rce - grid_id & ktau added to arg list
1342       subroutine activate(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair,  &
1343                       msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer,    &
1344                       na, volc, dlo_sect,dhi_sect,sigman, hygro, &
1345                       fn, fs, fm, fluxn, fluxs, fluxm, flux_fullact, &
1346                       grid_id, ktau, ii, jj, kk )
1347
1348 !      calculates number, surface, and mass fraction of aerosols activated as CCN
1349 !      calculates flux of cloud droplets, surface area, and aerosol mass into cloud
1350 !      assumes an internal mixture within each of aerosol mode.
1351 !      A sectional treatment within each type is assumed if ntype_aer >7.
1352 !      A gaussiam spectrum of updrafts can be treated.
1353
1354 !      mks units
1355
1356 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1357 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1358
1359       USE module_model_constants, only: g,rhowater, xlv, cp, rvovrd, r_d, r_v, &
1360               mwdry,svp1,svp2,svp3,ep_2
1361
1362       implicit none
1363
1364
1365 !      input
1366
1367       integer,intent(in) :: maxd_atype      ! dimension of types
1368       integer,intent(in) :: maxd_asize      ! dimension of sizes
1369       integer,intent(in) :: ntype_aer       ! number of types
1370       integer,intent(in) :: nsize_aer(maxd_atype) ! number of sizes for type
1371       integer,intent(in) :: msectional      ! 1 for sectional, 0 for modal
1372       integer,intent(in) :: grid_id         ! WRF grid%id
1373       integer,intent(in) :: ktau            ! WRF time step count
1374       integer,intent(in) :: ii, jj, kk      ! i,j,k of current grid cell
1375       real,intent(in) :: wbar          ! grid cell mean vertical velocity (m/s)
1376       real,intent(in) :: sigw          ! subgrid standard deviation of vertical vel (m/s)
1377       real,intent(in) :: wdiab         ! diabatic vertical velocity (0 if adiabatic)
1378       real,intent(in) :: wminf         ! minimum updraft velocity for integration (m/s)
1379       real,intent(in) :: wmaxf         ! maximum updraft velocity for integration (m/s)
1380       real,intent(in) :: tair          ! air temperature (K)
1381       real,intent(in) :: rhoair        ! air density (kg/m3)
1382       real,intent(in) :: na(maxd_asize,maxd_atype)     ! aerosol number concentration (/m3)
1383       real,intent(in) :: sigman(maxd_asize,maxd_atype) ! geometric standard deviation of aerosol size distribution
1384       real,intent(in) :: hygro(maxd_asize,maxd_atype)  ! bulk hygroscopicity of aerosol mode
1385       real,intent(in) :: volc(maxd_asize,maxd_atype)   ! total aerosol volume  concentration (m3/m3)
1386       real,intent(in) :: dlo_sect( maxd_asize, maxd_atype ), &  ! minimum size of section (cm)
1387            dhi_sect( maxd_asize, maxd_atype )     ! maximum size of section (cm)
1388
1389 !      output
1390
1391       real,intent(inout) :: fn(maxd_asize,maxd_atype)    ! number fraction of aerosols activated
1392       real,intent(inout) :: fs(maxd_asize,maxd_atype)    ! surface fraction of aerosols activated
1393       real,intent(inout) :: fm(maxd_asize,maxd_atype)    ! mass fraction of aerosols activated
1394       real,intent(inout) :: fluxn(maxd_asize,maxd_atype) ! flux of activated aerosol number fraction into cloud (m/s)
1395       real,intent(inout) :: fluxs(maxd_asize,maxd_atype) ! flux of activated aerosol surface fraction (m/s)
1396       real,intent(inout) :: fluxm(maxd_asize,maxd_atype) ! flux of activated aerosol mass fraction into cloud (m/s)
1397       real,intent(inout) :: flux_fullact                 ! flux when activation fraction = 100% (m/s)
1398
1399 !      local
1400
1401 !!$      external erf,erfc
1402 !!$      real erf,erfc
1403 !      external qsat_water
1404       integer, parameter:: nx=200
1405       integer iquasisect_option, isectional
1406       real integ,integf
1407       real, parameter :: surften = 0.076 ! surface tension of water w/respect to air (N/m)
1408       real, parameter :: p0 = 1013.25e2  ! reference pressure (Pa)
1409       real, parameter :: t0 = 273.15     ! reference temperature (K)
1410       real ylo(maxd_asize,maxd_atype),yhi(maxd_asize,maxd_atype) ! 1-particle volume at section interfaces
1411       real ymean(maxd_asize,maxd_atype) ! 1-particle volume at r=rmean
1412       real ycut, lnycut, betayy, betayy2, gammayy, phiyy
1413       real surfc(maxd_asize,maxd_atype) ! surface concentration (m2/m3)
1414       real sign(maxd_asize,maxd_atype)    ! geometric standard deviation of size distribution
1415       real alnsign(maxd_asize,maxd_atype) ! natl log of geometric standard dev of aerosol
1416       real am(maxd_asize,maxd_atype) ! number mode radius of dry aerosol (m)
1417       real lnhygro(maxd_asize,maxd_atype) ! ln(b)
1418       real f1(maxd_asize,maxd_atype) ! array to hold parameter for maxsat
1419       real pres ! pressure (Pa)
1420       real path ! mean free path (m)
1421       real diff ! diffusivity (m2/s)
1422       real conduct ! thermal conductivity (Joule/m/sec/deg)
1423       real diff0,conduct0
1424       real es ! saturation vapor pressure
1425       real qs ! water vapor saturation mixing ratio
1426       real dqsdt ! change in qs with temperature
1427       real dqsdp ! change in qs with pressure
1428       real gg ! thermodynamic function (m2/s)
1429       real sqrtg ! sqrt(gg)
1430       real sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
1431       real lnsm(maxd_asize,maxd_atype) ! ln( sm )
1432       real zeta, eta(maxd_asize,maxd_atype)
1433       real lnsmax ! ln(smax)
1434       real alpha
1435       real gamma
1436       real beta
1437       real gaus
1438       logical :: top        ! true if cloud top, false if cloud base or new cloud
1439       real asub(maxd_asize,maxd_atype),bsub(maxd_asize,maxd_atype) ! coefficients of submode size distribution N=a+bx
1440       real totn(maxd_atype) ! total aerosol number concentration
1441       real aten ! surface tension parameter
1442       real gmrad(maxd_atype) ! geometric mean radius
1443       real gmradsq(maxd_atype) ! geometric mean of radius squared
1444       real gmlnsig(maxd_atype) ! geometric standard deviation
1445       real gmsm(maxd_atype) ! critical supersaturation at radius gmrad
1446       real sumflxn(maxd_asize,maxd_atype)
1447       real sumflxs(maxd_asize,maxd_atype)
1448       real sumflxm(maxd_asize,maxd_atype)
1449       real sumflx_fullact
1450       real sumfn(maxd_asize,maxd_atype)
1451       real sumfs(maxd_asize,maxd_atype)
1452       real sumfm(maxd_asize,maxd_atype)
1453       real sumns(maxd_atype)
1454       real fnold(maxd_asize,maxd_atype)   ! number fraction activated
1455       real fsold(maxd_asize,maxd_atype)   ! surface fraction activated
1456       real fmold(maxd_asize,maxd_atype)   ! mass fraction activated
1457       real wold,gold
1458       real alogten,alog2,alog3,alogaten
1459       real alogam
1460       real rlo(maxd_asize,maxd_atype), rhi(maxd_asize,maxd_atype)
1461       real rmean(maxd_asize,maxd_atype)
1462                   ! mean radius (m) for the section (not used with modal)
1463                   ! calculated from current volume & number
1464       real ccc
1465       real dumaa,dumbb
1466       real wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
1467       real dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar
1468       real alw,sqrtalw
1469       real smax
1470       real x,arg
1471       real xmincoeff,xcut
1472       real z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
1473       real etafactor1,etafactor2(maxd_asize,maxd_atype),etafactor2max
1474       integer m,n,nw,nwmax
1475
1476 !      numerical integration parameters
1477       real, parameter :: eps = 0.3
1478       real, parameter :: fmax = 0.99
1479       real, parameter :: sds = 3.
1480
1481 !      mathematical constants
1482       real third, twothird, sixth, zero, one, two, three
1483
1484       real, parameter :: sq2  = 1.4142135624
1485       real, parameter :: sqpi = 1.7724538509
1486       real, parameter :: pi   = 3.1415926536
1487
1488 !      integer, save :: ndist(nx)  ! accumulates frequency distribution of integration bins required
1489 !      data ndist/nx*0/
1490
1491 !     for nsize_aer>7, a sectional approach is used and isectional = iquasisect_option
1492 !     activation fractions (fn,fs,fm) are computed as follows
1493 !     iquasisect_option = 1,3 - each section treated as a narrow lognormal
1494 !     iquasisect_option = 2,4 - within-section dn/dx = a + b*x,  x = ln(r)
1495 !     smax is computed as follows (when explicit activation is OFF)
1496 !     iquasisect_option = 1,2 - razzak-ghan modal parameterization with
1497 !     single mode having same ntot, dgnum, sigmag as the combined sections
1498 !     iquasisect_option = 3,4 - razzak-ghan sectional parameterization
1499 !     for nsize_aer=<9, a modal approach is used and isectional = 0
1500
1501 ! rce 08-jul-2005
1502 ! if either (na(n,m) < nsmall) or (volc(n,m) < vsmall)
1503 ! then treat bin/mode (n,m) as being empty, and set its fn/fs/fm=0.0
1504 !     (for single precision, gradual underflow starts around 1.0e-38,
1505 !      and strange things can happen when in that region)
1506       real, parameter :: nsmall = 1.0e-20    ! aer number conc in #/m3
1507       real, parameter :: vsmall = 1.0e-37    ! aer volume conc in m3/m3
1508       logical bin_is_empty(maxd_asize,maxd_atype), all_bins_empty
1509       logical bin_is_narrow(maxd_asize,maxd_atype)
1510
1511       integer idiagaa, ipass_nwloop
1512       integer idiag_dndy_neg, idiag_fnsm_prob
1513
1514 ! The flag for cloud top is no longer used so set it to false. This is an
1515 ! antiquated feature related to radiation enhancing mass fluxes at cloud
1516 ! top. It is currently, as of Feb. 2009, set to false in the CAM version
1517 ! as well.
1518       top = .false.
1519
1520 !.......................................................................
1521 !
1522 !   start calc. of modal or sectional activation properties (start of section 1)
1523 !
1524 !.......................................................................
1525       idiag_dndy_neg = 1      ! set this to 0 to turn off
1526                               !     warnings about dn/dy < 0
1527       idiag_fnsm_prob = 1     ! set this to 0 to turn off
1528                               !     warnings about fn/fs/fm misbehavior
1529
1530       iquasisect_option = 2
1531       if(msectional.gt.0)then
1532          isectional = iquasisect_option
1533       else
1534          isectional = 0
1535       endif
1536
1537       do n=1,ntype_aer
1538 !         print *,'ntype_aer,n,nsize_aer(n)=',ntype_aer,n,nsize_aer(n)
1539
1540         if(ntype_aer.eq.1.and.nsize_aer(n).eq.1.and.na(1,1).lt.1.e-20)then
1541          fn(1,1)=0.
1542          fs(1,1)=0.
1543          fm(1,1)=0.
1544          fluxn(1,1)=0.
1545          fluxs(1,1)=0.
1546          fluxm(1,1)=0.
1547          flux_fullact=0.
1548          return
1549         endif
1550       enddo
1551
1552       zero = 0.0
1553       one = 1.0
1554       two = 2.0
1555       three = 3.0
1556       third = 1.0/3.0
1557       twothird = 2.0/3.0 !wig, 1-Mar-2009: Corrected value from 2/6
1558       sixth = 1.0/6.0
1559
1560       pres=r_d*rhoair*tair
1561       diff0=0.211e-4*(p0/pres)*(tair/t0)**1.94
1562       conduct0=(5.69+0.017*(tair-t0))*4.186e2*1.e-5 ! convert to J/m/s/deg
1563       es=1000.*svp1*exp( svp2*(tair-t0)/(tair-svp3) )
1564       qs=ep_2*es/(pres-es)
1565       dqsdt=xlv/(r_v*tair*tair)*qs
1566       alpha=g*(xlv/(cp*r_v*tair*tair)-1./(r_d*tair))
1567       gamma=(1+xlv/cp*dqsdt)/(rhoair*qs)
1568       gg=1./(rhowater/(diff0*rhoair*qs)+xlv*rhowater/(conduct0*tair)*(xlv/(r_v*tair)-1.))
1569       sqrtg=sqrt(gg)
1570       beta=4.*pi*rhowater*gg*gamma
1571       aten=2.*surften/(r_v*tair*rhowater)
1572       alogaten=log(aten)
1573       alog2=log(two)
1574       alog3=log(three)
1575       ccc=4.*pi*third
1576       etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small.
1577
1578       all_bins_empty = .true.
1579       do n=1,ntype_aer
1580       totn(n)=0.
1581       gmrad(n)=0.
1582       gmradsq(n)=0.
1583       sumns(n)=0.
1584       do m=1,nsize_aer(n)
1585          alnsign(m,n)=log(sigman(m,n))
1586 !         internal mixture of aerosols
1587
1588          bin_is_empty(m,n) = .true.
1589          if (volc(m,n).gt.vsmall .and. na(m,n).gt.nsmall) then
1590             bin_is_empty(m,n) = .false.
1591             all_bins_empty = .false.
1592             lnhygro(m,n)=log(hygro(m,n))
1593 !            number mode radius (m,n)
1594 !           write(6,*)'alnsign,volc,na=',alnsign(m,n),volc(m,n),na(m,n)
1595             am(m,n)=exp(-1.5*alnsign(m,n)*alnsign(m,n))*              &
1596               (3.*volc(m,n)/(4.*pi*na(m,n)))**third
1597
1598             if (isectional .gt. 0) then
1599 !               sectional model.
1600 !               need to use bulk properties because parameterization doesn't
1601 !               work well for narrow bins.
1602                totn(n)=totn(n)+na(m,n)
1603                alogam=log(am(m,n))
1604                gmrad(n)=gmrad(n)+na(m,n)*alogam
1605                gmradsq(n)=gmradsq(n)+na(m,n)*alogam*alogam
1606             endif
1607             etafactor2(m,n)=1./(na(m,n)*beta*sqrtg)
1608
1609             if(hygro(m,n).gt.1.e-10)then
1610                sm(m,n)=2.*aten/(3.*am(m,n))*sqrt(aten/(3.*hygro(m,n)*am(m,n)))
1611             else
1612                sm(m,n)=100.
1613             endif
1614 !           write(6,*)'sm,hygro,am=',sm(m,n),hygro(m,n),am(m,n)
1615          else
1616             sm(m,n)=1.
1617             etafactor2(m,n)=etafactor2max ! this should make eta big if na is very small.
1618
1619          endif
1620          lnsm(m,n)=log(sm(m,n))
1621          if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
1622             sumns(n)=sumns(n)+na(m,n)/sm(m,n)**twothird
1623          endif
1624 !        write(6,'(a,i4,6g12.2)')'m,na,am,hygro,lnhygro,sm,lnsm=',m,na(m,n),am(m,n),hygro(m,n),lnhygro(m,n),sm(m,n),lnsm(m,n)
1625       end do ! size
1626       end do ! type
1627
1628 !  if all bins are empty, set all activation fractions to zero and exit
1629          if ( all_bins_empty ) then
1630             do n=1,ntype_aer
1631             do m=1,nsize_aer(n)
1632                fluxn(m,n)=0.
1633                fn(m,n)=0.
1634                fluxs(m,n)=0.
1635                fs(m,n)=0.
1636                fluxm(m,n)=0.
1637                fm(m,n)=0.
1638             end do
1639             end do
1640             flux_fullact=0.
1641             return
1642          endif
1643
1644
1645
1646          if (isectional .le. 0) then
1647             ! Initialize maxsat at this cell and timestep for the
1648             ! modal setup (the sectional case is handled below).
1649             call maxsat_init(maxd_atype, ntype_aer, &
1650                  maxd_asize, nsize_aer, alnsign, f1)
1651
1652             goto 30000
1653          end if
1654
1655          do n=1,ntype_aer
1656             !wig 19-Oct-2006: Add zero trap based May 2006 e-mail from
1657             !Ghan. Transport can clear out a cell leading to
1658             !inconsistencies with the mass.
1659             gmrad(n)=gmrad(n)/max(totn(n),1e-20)
1660             gmlnsig=gmradsq(n)/totn(n)-gmrad(n)*gmrad(n)    ! [ln(sigmag)]**2
1661             gmlnsig(n)=sqrt( max( 1.e-4, gmlnsig(n) ) )
1662             gmrad(n)=exp(gmrad(n))
1663             if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
1664                gmsm(n)=totn(n)/sumns(n)
1665                gmsm(n)=gmsm(n)*gmsm(n)*gmsm(n)
1666                gmsm(n)=sqrt(gmsm(n))
1667             else
1668 !              gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(1,n)*gmrad(n)))
1669                gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(nsize_aer(n),n)*gmrad(n)))
1670             endif
1671          enddo
1672
1673          ! Initialize maxsat at this cell and timestep for the
1674          ! sectional setup (the modal case is handled above)...
1675          call maxsat_init(maxd_atype, ntype_aer, &
1676               maxd_asize, (/1/), gmlnsig, f1)
1677
1678 !.......................................................................
1679 !   calculate sectional "sub-bin" size distribution
1680 !
1681 !   dn/dy = nt*( a + b*y )   for  ylo < y < yhi
1682 !
1683 !   nt = na(m,n) = number mixing ratio of the bin
1684 !   y = v/vhi
1685 !       v = (4pi/3)*r**3 = particle volume
1686 !       vhi = v at r=rhi (upper bin boundary)
1687 !   ylo = y at lower bin boundary = vlo/vhi = (rlo/rhi)**3
1688 !   yhi = y at upper bin boundary = 1.0
1689 !
1690 !   dv/dy = v * dn/dy = nt*vhi*( a*y + b*y*y )
1691 !
1692 !.......................................................................
1693 ! 02-may-2006 - this dn/dy replaces the previous
1694 !       dn/dx = a + b*x   where l = ln(r)
1695 !    the old dn/dx was overly complicated for cases of rmean near rlo or rhi
1696 !    the new dn/dy is consistent with that used in the movesect routine,
1697 !       which does continuous growth by condensation and aqueous chemistry
1698 !.......................................................................
1699              do 25002 n = 1,ntype_aer
1700              do 25000 m = 1,nsize_aer(n)
1701
1702 ! convert from diameter in cm to radius in m
1703                 rlo(m,n) = 0.5*0.01*dlo_sect(m,n)
1704                 rhi(m,n) = 0.5*0.01*dhi_sect(m,n)
1705                 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1706                 yhi(m,n) = 1.0
1707
1708 ! 04-nov-2005 - extremely narrow bins will be treated using 0/1 activation
1709 !    this is to avoid potential numerical problems
1710                 bin_is_narrow(m,n) = .false.
1711                 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) bin_is_narrow(m,n) = .true.
1712
1713 ! rmean is mass mean radius for the bin; xmean = log(rmean)
1714 ! just use section midpoint if bin is empty
1715                 if ( bin_is_empty(m,n) ) then
1716                    rmean(m,n) = sqrt(rlo(m,n)*rhi(m,n))
1717                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1718                    goto 25000
1719                 end if
1720
1721                 rmean(m,n) = (volc(m,n)/(ccc*na(m,n)))**third
1722                 rmean(m,n) = max( rlo(m,n), min( rhi(m,n), rmean(m,n) ) )
1723                 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1724                 if ( bin_is_narrow(m,n) ) goto 25000
1725
1726 ! if rmean is extremely close to either rlo or rhi,
1727 ! treat the bin as extremely narrow
1728                 if ((rhi(m,n)/rmean(m,n)) .le. 1.01) then
1729                    bin_is_narrow(m,n) = .true.
1730                    rlo(m,n) = min( rmean(m,n), (rhi(m,n)/1.01) )
1731                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1732                    goto 25000
1733                 else if ((rmean(m,n)/rlo(m,n)) .le. 1.01) then
1734                    bin_is_narrow(m,n) = .true.
1735                    rhi(m,n) = max( rmean(m,n), (rlo(m,n)*1.01) )
1736                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1737                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1738                    goto 25000
1739                 endif
1740
1741 ! if rmean is somewhat close to either rlo or rhi, then dn/dy will be
1742 !    negative near the upper or lower bin boundary
1743 ! in these cases, assume that all the particles are in a subset of the full bin,
1744 !    and adjust rlo or rhi so that rmean will be near the center of this subset
1745 ! note that the bin is made narrower LOCALLY/TEMPORARILY,
1746 !    just for the purposes of the activation calculation
1747                 gammayy = (ymean(m,n)-ylo(m,n)) / (yhi(m,n)-ylo(m,n))
1748                 if (gammayy .lt. 0.34) then
1749                    dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*(gammayy/0.34)
1750                    rhi(m,n) = rhi(m,n)*(dumaa**third)
1751                    ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1752                    ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1753                 else if (gammayy .ge. 0.66) then
1754                    dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*((gammayy-0.66)/0.34)
1755                    ylo(m,n) = dumaa
1756                    rlo(m,n) = rhi(m,n)*(dumaa**third)
1757                 end if
1758                 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) then
1759                    bin_is_narrow(m,n) = .true.
1760                    goto 25000
1761                 end if
1762
1763                 betayy = ylo(m,n)/yhi(m,n)
1764                 betayy2 = betayy*betayy
1765                 bsub(m,n) = (12.0*ymean(m,n) - 6.0*(1.0+betayy)) /   &
1766                    (4.0*(1.0-betayy2*betayy) - 3.0*(1.0-betayy2)*(1.0+betayy))
1767                 asub(m,n) = (1.0 - bsub(m,n)*(1.0-betayy2)*0.5) / (1.0-betayy)
1768
1769                 if ( asub(m,n)+bsub(m,n)*ylo(m,n) .lt. 0. ) then
1770                   if (idiag_dndy_neg .gt. 0) then
1771                     print *,'dndy<0 at lower boundary'
1772                     print *,'n,m=',n,m
1773                     print *,'na=',na(m,n),' volc=',volc(m,n)
1774                     print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
1775                     print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
1776                     print *,'dlo_sect/2,dhi_sect/2=',   &
1777                              (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
1778                     print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
1779                     print *,'asub+bsub*ylo=',   &
1780                              (asub(m,n)+bsub(m,n)*ylo(m,n))
1781                     print *,'subr activate error 11 - i,j,k =', ii, jj, kk
1782                   endif
1783                 endif
1784                 if ( asub(m,n)+bsub(m,n)*yhi(m,n) .lt. 0. ) then
1785                   if (idiag_dndy_neg .gt. 0) then
1786                     print *,'dndy<0 at upper boundary'
1787                     print *,'n,m=',n,m
1788                     print *,'na=',na(m,n),' volc=',volc(m,n)
1789                     print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
1790                     print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
1791                     print *,'dlo_sect/2,dhi_sect/2=',   &
1792                              (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
1793                     print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
1794                     print *,'asub+bsub*yhi=',   &
1795                              (asub(m,n)+bsub(m,n)*yhi(m,n))
1796                     print *,'subr activate error 12 - i,j,k =', ii, jj, kk
1797                   endif
1798                 endif
1799
1800 25000        continue      ! m=1,nsize_aer(n)
1801 25002        continue      ! n=1,ntype_aer
1802
1803
1804 30000    continue
1805 !.......................................................................
1806 !
1807 !   end calc. of modal or sectional activation properties (end of section 1)
1808 !
1809 !.......................................................................
1810
1811
1812
1813 !      sjg 7-16-98  upward
1814 !      print *,'wbar,sigw=',wbar,sigw
1815
1816       if(sigw.le.1.e-5) goto 50000
1817
1818 !.......................................................................
1819 !
1820 !   start calc. of activation fractions/fluxes
1821 !   for spectrum of updrafts (start of section 2)
1822 !
1823 !.......................................................................
1824          ipass_nwloop = 1
1825          idiagaa = 0
1826 ! 06-nov-2005 rce - set idiagaa=1 for testing/debugging
1827 !        if ((grid_id.eq.1) .and. (ktau.eq.167) .and.   &
1828 !            (ii.eq.24) .and. (jj.eq. 1) .and. (kk.eq.14)) idiagaa = 1
1829
1830 40000    continue
1831          if(top)then
1832            wmax=0.
1833            wmin=min(zero,-wdiab)
1834          else
1835            wmax=min(wmaxf,wbar+sds*sigw)
1836            wmin=max(wminf,-wdiab)
1837          endif
1838          wmin=max(wmin,wbar-sds*sigw)
1839          w=wmin
1840          dwmax=eps*sigw
1841          dw=dwmax
1842          dfmax=0.2
1843          dfmin=0.1
1844          if(wmax.le.w)then
1845             do n=1,ntype_aer
1846             do m=1,nsize_aer(n)
1847                fluxn(m,n)=0.
1848                fn(m,n)=0.
1849                fluxs(m,n)=0.
1850                fs(m,n)=0.
1851                fluxm(m,n)=0.
1852                fm(m,n)=0.
1853             end do
1854             end do
1855             flux_fullact=0.
1856             return
1857          endif
1858          do n=1,ntype_aer
1859          do m=1,nsize_aer(n)
1860             sumflxn(m,n)=0.
1861             sumfn(m,n)=0.
1862             fnold(m,n)=0.
1863             sumflxs(m,n)=0.
1864             sumfs(m,n)=0.
1865             fsold(m,n)=0.
1866             sumflxm(m,n)=0.
1867             sumfm(m,n)=0.
1868             fmold(m,n)=0.
1869          enddo
1870          enddo
1871          sumflx_fullact=0.
1872
1873          fold=0
1874          gold=0
1875 ! 06-nov-2005 rce - set wold=w here
1876 !        wold=0
1877          wold=w
1878
1879
1880 ! 06-nov-2005 rce - define nwmax; calc dwmin from nwmax
1881          nwmax = 200
1882 !        dwmin = min( dwmax, 0.01 )
1883          dwmin = (wmax - wmin)/(nwmax-1)
1884          dwmin = min( dwmax, dwmin )
1885          dwmin = max( 0.01,  dwmin )
1886
1887 !
1888 ! loop over updrafts, incrementing sums as you go
1889 ! the "200" is (arbitrary) upper limit for number of updrafts
1890 ! if integration finishes before this, OK; otherwise, ERROR
1891 !
1892          if (idiagaa.gt.0) then
1893              write(*,94700) ktau, grid_id, ii, jj, kk, nwmax
1894              write(*,94710) 'wbar,sigw,wdiab=', wbar, sigw, wdiab
1895              write(*,94710) 'wmin,wmax,dwmin,dwmax=', wmin, wmax, dwmin, dwmax
1896              write(*,94720) -1, w, wold, dw
1897          end if
1898 94700    format( / 'activate 47000 - ktau,id,ii,jj,kk,nwmax=', 6i5 )
1899 94710    format( 'activate 47000 - ', a, 6(1x,f11.5) )
1900 94720    format( 'activate 47000 - nw,w,wold,dw=', i5, 3(1x,f11.5) )
1901
1902          do 47000 nw = 1, nwmax
1903 41000       wnuc=w+wdiab
1904
1905             if (idiagaa.gt.0) write(*,94720) nw, w, wold, dw
1906
1907 !           write(6,*)'wnuc=',wnuc
1908             alw=alpha*wnuc
1909             sqrtalw=sqrt(alw)
1910             zeta=2.*sqrtalw*aten/(3.*sqrtg)
1911             etafactor1=2.*alw*sqrtalw
1912             if (isectional .gt. 0) then
1913 !              sectional model.
1914 !              use bulk properties
1915
1916               do n=1,ntype_aer
1917                  if(totn(n).gt.1.e-10)then
1918                     eta(1,n)=etafactor1/(totn(n)*beta*sqrtg)
1919                  else
1920                     eta(1,n)=1.e10
1921                  endif
1922               enddo
1923               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
1924                    maxd_asize,(/1/),gmsm,gmlnsig,f1,smax)
1925               lnsmax=log(smax)
1926               x=2*(log(gmsm(1))-lnsmax)/(3*sq2*gmlnsig(1))
1927               fnew=0.5*(1.-ERF_ALT(x))
1928
1929             else
1930
1931               do n=1,ntype_aer
1932               do m=1,nsize_aer(n)
1933                  eta(m,n)=etafactor1*etafactor2(m,n)
1934               enddo
1935               enddo
1936
1937               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
1938                    maxd_asize,nsize_aer,sm,alnsign,f1,smax)
1939 !             write(6,*)'w,smax=',w,smax
1940
1941               lnsmax=log(smax)
1942
1943               x=2*(lnsm(nsize_aer(1),1)-lnsmax)/(3*sq2*alnsign(nsize_aer(1),1))
1944               fnew=0.5*(1.-ERF_ALT(x))
1945
1946             endif
1947
1948             dwnew = dw
1949 ! 06-nov-2005 rce - "n" here should be "nw" (?)
1950 !           if(fnew-fold.gt.dfmax.and.n.gt.1)then
1951             if(fnew-fold.gt.dfmax.and.nw.gt.1)then
1952 !              reduce updraft increment for greater accuracy in integration
1953                if (dw .gt. 1.01*dwmin) then
1954                   dw=0.7*dw
1955                   dw=max(dw,dwmin)
1956                   w=wold+dw
1957                   go to 41000
1958                else
1959                   dwnew = dwmin
1960                endif
1961             endif
1962
1963             if(fnew-fold.lt.dfmin)then
1964 !              increase updraft increment to accelerate integration
1965                dwnew=min(1.5*dw,dwmax)
1966             endif
1967             fold=fnew
1968
1969             z=(w-wbar)/(sigw*sq2)
1970             gaus=exp(-z*z)
1971             fnmin=1.
1972             xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
1973 !           write(6,*)'xmincoeff=',xmincoeff
1974
1975
1976             do 44002 n=1,ntype_aer
1977             do 44000 m=1,nsize_aer(n)
1978                if ( bin_is_empty(m,n) ) then
1979                    fn(m,n)=0.
1980                    fs(m,n)=0.
1981                    fm(m,n)=0.
1982                else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
1983 !                 sectional
1984 !                  within-section dn/dx = a + b*x
1985                   xcut=xmincoeff-third*lnhygro(m,n)
1986 !                 ycut=(exp(xcut)/rhi(m,n))**3
1987 ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
1988 ! if (ycut > yhi), then actual value of ycut is unimportant,
1989 ! so do the following to avoid overflow
1990                   lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
1991                   lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
1992                   ycut=exp(lnycut)
1993 !                 write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
1994 !                   if(lnsmax.lt.lnsmn(m,n))then
1995                   if(ycut.gt.yhi(m,n))then
1996                      fn(m,n)=0.
1997                      fs(m,n)=0.
1998                      fm(m,n)=0.
1999                   elseif(ycut.lt.ylo(m,n))then
2000                      fn(m,n)=1.
2001                      fs(m,n)=1.
2002                      fm(m,n)=1.
2003                   elseif ( bin_is_narrow(m,n) ) then
2004 ! 04-nov-2005 rce - for extremely narrow bins,
2005 ! do zero activation if xcut>xmean, 100% activation otherwise
2006                      if (ycut.gt.ymean(m,n)) then
2007                         fn(m,n)=0.
2008                         fs(m,n)=0.
2009                         fm(m,n)=0.
2010                      else
2011                         fn(m,n)=1.
2012                         fs(m,n)=1.
2013                         fm(m,n)=1.
2014                      endif
2015                   else
2016                      phiyy=ycut/yhi(m,n)
2017                      fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
2018                      if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
2019                       if (idiag_fnsm_prob .gt. 0) then
2020                         print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
2021                         print *,'na,volc       =', na(m,n), volc(m,n)
2022                         print *,'asub,bsub     =', asub(m,n), bsub(m,n)
2023                         print *,'yhi,ycut      =', yhi(m,n), ycut
2024                       endif
2025                      endif
2026
2027                      if (fn(m,n) .le. zero) then
2028 ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
2029                         fn(m,n)=zero
2030                         fs(m,n)=zero
2031                         fm(m,n)=zero
2032                      else if (fn(m,n) .ge. one) then
2033 ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
2034                         fn(m,n)=one
2035                         fs(m,n)=one
2036                         fm(m,n)=one
2037                      else
2038 ! 10-nov-2005 rce - otherwise, calc fm and check it
2039                         fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) +   &
2040                                   third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
2041                         if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
2042                          if (idiag_fnsm_prob .gt. 0) then
2043                            print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
2044                            print *,'na,volc,fn    =', na(m,n), volc(m,n), fn(m,n)
2045                            print *,'asub,bsub     =', asub(m,n), bsub(m,n)
2046                            print *,'yhi,ycut     =', yhi(m,n), ycut
2047                          endif
2048                         endif
2049                         if (fm(m,n) .le. fn(m,n)) then
2050 ! 10-nov-2005 rce - if fm=fn, then fs must =fn
2051                            fm(m,n)=fn(m,n)
2052                            fs(m,n)=fn(m,n)
2053                         else if (fm(m,n) .ge. one) then
2054 ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
2055                            fm(m,n)=one
2056                            fs(m,n)=one
2057                            fn(m,n)=one
2058                         else
2059 ! 10-nov-2005 rce - these two checks assure that the mean size
2060 ! of the activated & interstitial particles will be between rlo & rhi
2061                            dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n))
2062                            fm(m,n) = min( fm(m,n), dumaa )
2063                            dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n))
2064                            fm(m,n) = min( fm(m,n), dumaa )
2065 ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
2066                            betayy = ylo(m,n)/yhi(m,n)
2067                            dumaa = phiyy**twothird
2068                            dumbb = betayy**twothird
2069                            fs(m,n) =   &
2070                               (asub(m,n)*(1.0-phiyy*dumaa) +   &
2071                                   0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) /   &
2072                               (asub(m,n)*(1.0-betayy*dumbb) +   &
2073                                   0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
2074                            fs(m,n)=max(fs(m,n),fn(m,n))
2075                            fs(m,n)=min(fs(m,n),fm(m,n))
2076                         endif
2077                      endif
2078                   endif
2079
2080                else
2081 !                 modal
2082                   x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
2083                   fn(m,n)=0.5*(1.-ERF_ALT(x))
2084                   arg=x-sq2*alnsign(m,n)
2085                   fs(m,n)=0.5*(1.-ERF_ALT(arg))
2086                   arg=x-1.5*sq2*alnsign(m,n)
2087                   fm(m,n)=0.5*(1.-ERF_ALT(arg))
2088 !                 print *,'w,x,fn,fs,fm=',w,x,fn(m,n),fs(m,n),fm(m,n)
2089                endif
2090
2091 !                     fn(m,n)=1.  !test
2092 !                     fs(m,n)=1.
2093 !                     fm(m,n)=1.
2094                fnmin=min(fn(m,n),fnmin)
2095 !               integration is second order accurate
2096 !               assumes linear variation of f*gaus with w
2097                wb=(w+wold)
2098                fnbar=(fn(m,n)*gaus+fnold(m,n)*gold)
2099                fsbar=(fs(m,n)*gaus+fsold(m,n)*gold)
2100                fmbar=(fm(m,n)*gaus+fmold(m,n)*gold)
2101                if((top.and.w.lt.0.).or.(.not.top.and.w.gt.0.))then
2102                   sumflxn(m,n)=sumflxn(m,n)+sixth*(wb*fnbar           &
2103                       +(fn(m,n)*gaus*w+fnold(m,n)*gold*wold))*dw
2104                   sumflxs(m,n)=sumflxs(m,n)+sixth*(wb*fsbar           &
2105                       +(fs(m,n)*gaus*w+fsold(m,n)*gold*wold))*dw
2106                   sumflxm(m,n)=sumflxm(m,n)+sixth*(wb*fmbar           &
2107                       +(fm(m,n)*gaus*w+fmold(m,n)*gold*wold))*dw
2108                endif
2109                sumfn(m,n)=sumfn(m,n)+0.5*fnbar*dw
2110 !              write(6,'(a,9g10.2)')'lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw=', &
2111 !                lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw
2112                fnold(m,n)=fn(m,n)
2113                sumfs(m,n)=sumfs(m,n)+0.5*fsbar*dw
2114                fsold(m,n)=fs(m,n)
2115                sumfm(m,n)=sumfm(m,n)+0.5*fmbar*dw
2116                fmold(m,n)=fm(m,n)
2117
2118 44000       continue      ! m=1,nsize_aer(n)
2119 44002       continue      ! n=1,ntype_aer
2120
2121 !           same form as sumflxm(m,n) but replace the fm/fmold(m,n) with 1.0
2122             sumflx_fullact = sumflx_fullact &
2123                            + sixth*(wb*(gaus+gold) + (gaus*w + gold*wold))*dw
2124 !            sumg=sumg+0.5*(gaus+gold)*dw
2125             gold=gaus
2126             wold=w
2127             dw=dwnew
2128
2129             if(nw.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 48000
2130             w=w+dw
2131
2132 47000    continue      ! nw = 1, nwmax
2133
2134
2135          print *,'do loop is too short in activate'
2136          print *,'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
2137          print *,'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
2138          print *,'wnuc=',wnuc
2139          do n=1,ntype_aer
2140             print *,'ntype=',n
2141             print *,'na=',(na(m,n),m=1,nsize_aer(n))
2142             print *,'fn=',(fn(m,n),m=1,nsize_aer(n))
2143          end do
2144 !   dump all subr parameters to allow testing with standalone code
2145 !   (build a driver that will read input and call activate)
2146          print *,'top,wbar,sigw,wdiab,tair,rhoair,ntype_aer='
2147          print *, top,wbar,sigw,wdiab,tair,rhoair,ntype_aer
2148          print *,'na='
2149          print *, na
2150          print *,'volc='
2151          print *, volc
2152          print *,'sigman='
2153          print *, sigman
2154          print *,'hygro='
2155          print *, hygro
2156
2157          print *,'subr activate error 31 - i,j,k =', ii, jj, kk
2158 ! 06-nov-2005 rce - if integration fails, repeat it once with additional diagnostics
2159          if (ipass_nwloop .eq. 1) then
2160              ipass_nwloop = 2
2161              idiagaa = 2
2162              goto 40000
2163          end if
2164          call wrf_error_fatal("STOP: activate before 48000")
2165
2166 48000    continue
2167
2168
2169 !         ndist(n)=ndist(n)+1
2170          if(.not.top.and.w.lt.wmaxf)then
2171
2172 !            contribution from all updrafts stronger than wmax
2173 !            assuming constant f (close to fmax)
2174             wnuc=w+wdiab
2175
2176             z1=(w-wbar)/(sigw*sq2)
2177             z2=(wmaxf-wbar)/(sigw*sq2)
2178             integ=sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(z1)-ERFC_NUM_RECIPES(z2))
2179 !            consider only upward flow into cloud base when estimating flux
2180             wf1=max(w,zero)
2181             zf1=(wf1-wbar)/(sigw*sq2)
2182             gf1=exp(-zf1*zf1)
2183             wf2=max(wmaxf,zero)
2184             zf2=(wf2-wbar)/(sigw*sq2)
2185             gf2=exp(-zf2*zf2)
2186             gf=(gf1-gf2)
2187             integf=wbar*sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(zf1)-ERFC_NUM_RECIPES(zf2))+sigw*sigw*gf
2188
2189             do n=1,ntype_aer
2190             do m=1,nsize_aer(n)
2191                sumflxn(m,n)=sumflxn(m,n)+integf*fn(m,n)
2192                sumfn(m,n)=sumfn(m,n)+fn(m,n)*integ
2193                sumflxs(m,n)=sumflxs(m,n)+integf*fs(m,n)
2194                sumfs(m,n)=sumfs(m,n)+fs(m,n)*integ
2195                sumflxm(m,n)=sumflxm(m,n)+integf*fm(m,n)
2196                sumfm(m,n)=sumfm(m,n)+fm(m,n)*integ
2197             end do
2198             end do
2199 !           same form as sumflxm(m,n) but replace the fm(m,n) with 1.0
2200             sumflx_fullact = sumflx_fullact + integf
2201 !            sumg=sumg+integ
2202          endif
2203
2204
2205          do n=1,ntype_aer
2206          do m=1,nsize_aer(n)
2207
2208 !           fn(m,n)=sumfn(m,n)/(sumg)
2209             fn(m,n)=sumfn(m,n)/(sq2*sqpi*sigw)
2210             fluxn(m,n)=sumflxn(m,n)/(sq2*sqpi*sigw)
2211             if(fn(m,n).gt.1.01)then
2212              if (idiag_fnsm_prob .gt. 0) then
2213                print *,'fn=',fn(m,n),' > 1 - activate err41'
2214                print *,'w,m,n,na,am=',w,m,n,na(m,n),am(m,n)
2215                print *,'integ,sumfn,sigw=',integ,sumfn(m,n),sigw
2216                print *,'subr activate error - i,j,k =', ii, jj, kk
2217              endif
2218              fluxn(m,n) = fluxn(m,n)/fn(m,n)
2219             endif
2220
2221             fs(m,n)=sumfs(m,n)/(sq2*sqpi*sigw)
2222             fluxs(m,n)=sumflxs(m,n)/(sq2*sqpi*sigw)
2223             if(fs(m,n).gt.1.01)then
2224              if (idiag_fnsm_prob .gt. 0) then
2225                print *,'fs=',fs(m,n),' > 1 - activate err42'
2226                print *,'m,n,isectional=',m,n,isectional
2227                print *,'alnsign(m,n)=',alnsign(m,n)
2228                print *,'rcut,rlo(m,n),rhi(m,n)',exp(xcut),rlo(m,n),rhi(m,n)
2229                print *,'w,m,na,am=',w,m,na(m,n),am(m,n)
2230                print *,'integ,sumfs,sigw=',integ,sumfs(m,n),sigw
2231              endif
2232              fluxs(m,n) = fluxs(m,n)/fs(m,n)
2233             endif
2234
2235 !           fm(m,n)=sumfm(m,n)/(sumg)
2236             fm(m,n)=sumfm(m,n)/(sq2*sqpi*sigw)
2237             fluxm(m,n)=sumflxm(m,n)/(sq2*sqpi*sigw)
2238             if(fm(m,n).gt.1.01)then
2239              if (idiag_fnsm_prob .gt. 0) then
2240                print *,'fm(',m,n,')=',fm(m,n),' > 1 - activate err43'
2241              endif
2242              fluxm(m,n) = fluxm(m,n)/fm(m,n)
2243             endif
2244
2245          end do
2246          end do
2247 !        same form as fluxm(m,n)
2248          flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
2249
2250       goto 60000
2251 !.......................................................................
2252 !
2253 !   end calc. of activation fractions/fluxes
2254 !   for spectrum of updrafts (end of section 2)
2255 !
2256 !.......................................................................
2257
2258 !.......................................................................
2259 !
2260 !   start calc. of activation fractions/fluxes
2261 !   for (single) uniform updraft (start of section 3)
2262 !
2263 !.......................................................................
2264 50000 continue
2265
2266          wnuc=wbar+wdiab
2267 !         write(6,*)'uniform updraft =',wnuc
2268
2269 ! 04-nov-2005 rce - moved the code for "wnuc.le.0" code to here
2270          if(wnuc.le.0.)then
2271             do n=1,ntype_aer
2272             do m=1,nsize_aer(n)
2273                fn(m,n)=0
2274                fluxn(m,n)=0
2275                fs(m,n)=0
2276                fluxs(m,n)=0
2277                fm(m,n)=0
2278                fluxm(m,n)=0
2279             end do
2280             end do
2281             flux_fullact=0.
2282             return
2283          endif
2284
2285             w=wbar
2286             alw=alpha*wnuc
2287             sqrtalw=sqrt(alw)
2288             zeta=2.*sqrtalw*aten/(3.*sqrtg)
2289
2290             if (isectional .gt. 0) then
2291 !              sectional model.
2292 !              use bulk properties
2293               do n=1,ntype_aer
2294               if(totn(n).gt.1.e-10)then
2295                  eta(1,n)=2*alw*sqrtalw/(totn(n)*beta*sqrtg)
2296               else
2297                  eta(1,n)=1.e10
2298               endif
2299               end do
2300                call maxsat(zeta,eta,maxd_atype,ntype_aer, &
2301                     maxd_asize,(/1/),gmsm,gmlnsig,f1,smax)
2302
2303             else
2304
2305               do n=1,ntype_aer
2306               do m=1,nsize_aer(n)
2307                  if(na(m,n).gt.1.e-10)then
2308                     eta(m,n)=2*alw*sqrtalw/(na(m,n)*beta*sqrtg)
2309                  else
2310                     eta(m,n)=1.e10
2311                  endif
2312               end do
2313               end do
2314
2315               call maxsat(zeta,eta,maxd_atype,ntype_aer, &
2316                    maxd_asize,nsize_aer,sm,alnsign,f1,smax)
2317
2318             endif
2319
2320             lnsmax=log(smax)
2321             xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
2322
2323             do 55002 n=1,ntype_aer
2324             do 55000 m=1,nsize_aer(n)
2325
2326 ! 04-nov-2005 rce - check for bin_is_empty here too, just like earlier
2327                if ( bin_is_empty(m,n) ) then
2328                    fn(m,n)=0.
2329                    fs(m,n)=0.
2330                    fm(m,n)=0.
2331
2332                else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
2333 !                 sectional
2334 !                  within-section dn/dx = a + b*x
2335                   xcut=xmincoeff-third*lnhygro(m,n)
2336 !                 ycut=(exp(xcut)/rhi(m,n))**3
2337 ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
2338 ! if (ycut > yhi), then actual value of ycut is unimportant,
2339 ! so do the following to avoid overflow
2340                   lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
2341                   lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
2342                   ycut=exp(lnycut)
2343 !                 write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
2344 !                   if(lnsmax.lt.lnsmn(m,n))then
2345                   if(ycut.gt.yhi(m,n))then
2346                      fn(m,n)=0.
2347                      fs(m,n)=0.
2348                      fm(m,n)=0.
2349 !                   elseif(lnsmax.gt.lnsmx(m,n))then
2350                   elseif(ycut.lt.ylo(m,n))then
2351                      fn(m,n)=1.
2352                      fs(m,n)=1.
2353                      fm(m,n)=1.
2354                   elseif ( bin_is_narrow(m,n) ) then
2355 ! 04-nov-2005 rce - for extremely narrow bins,
2356 ! do zero activation if xcut>xmean, 100% activation otherwise
2357                      if (ycut.gt.ymean(m,n)) then
2358                         fn(m,n)=0.
2359                         fs(m,n)=0.
2360                         fm(m,n)=0.
2361                      else
2362                         fn(m,n)=1.
2363                         fs(m,n)=1.
2364                         fm(m,n)=1.
2365                      endif
2366                   else
2367                      phiyy=ycut/yhi(m,n)
2368                      fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
2369                      if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
2370                       if (idiag_fnsm_prob .gt. 0) then
2371                         print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
2372                         print *,'na,volc       =', na(m,n), volc(m,n)
2373                         print *,'asub,bsub     =', asub(m,n), bsub(m,n)
2374                         print *,'yhi,ycut      =', yhi(m,n), ycut
2375                       endif
2376                      endif
2377
2378                      if (fn(m,n) .le. zero) then
2379 ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
2380                         fn(m,n)=zero
2381                         fs(m,n)=zero
2382                         fm(m,n)=zero
2383                      else if (fn(m,n) .ge. one) then
2384 ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
2385                         fn(m,n)=one
2386                         fs(m,n)=one
2387                         fm(m,n)=one
2388                      else
2389 ! 10-nov-2005 rce - otherwise, calc fm and check it
2390                         fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) +   &
2391                                   third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
2392                         if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
2393                          if (idiag_fnsm_prob .gt. 0) then
2394                            print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
2395                            print *,'na,volc,fn    =', na(m,n), volc(m,n), fn(m,n)
2396                            print *,'asub,bsub     =', asub(m,n), bsub(m,n)
2397                            print *,'yhi,ycut      =', yhi(m,n), ycut
2398                          endif
2399                         endif
2400                         if (fm(m,n) .le. fn(m,n)) then
2401 ! 10-nov-2005 rce - if fm=fn, then fs must =fn
2402                            fm(m,n)=fn(m,n)
2403                            fs(m,n)=fn(m,n)
2404                         else if (fm(m,n) .ge. one) then
2405 ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
2406                            fm(m,n)=one
2407                            fs(m,n)=one
2408                            fn(m,n)=one
2409                         else
2410 ! 10-nov-2005 rce - these two checks assure that the mean size
2411 ! of the activated & interstitial particles will be between rlo & rhi
2412                            dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n))
2413                            fm(m,n) = min( fm(m,n), dumaa )
2414                            dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n))
2415                            fm(m,n) = min( fm(m,n), dumaa )
2416 ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
2417                            betayy = ylo(m,n)/yhi(m,n)
2418                            dumaa = phiyy**twothird
2419                            dumbb = betayy**twothird
2420                            fs(m,n) =   &
2421                               (asub(m,n)*(1.0-phiyy*dumaa) +   &
2422                                   0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) /   &
2423                               (asub(m,n)*(1.0-betayy*dumbb) +   &
2424                                   0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
2425                            fs(m,n)=max(fs(m,n),fn(m,n))
2426                            fs(m,n)=min(fs(m,n),fm(m,n))
2427                         endif
2428                      endif
2429
2430                   endif
2431
2432                else
2433 !                 modal
2434                   x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
2435                   fn(m,n)=0.5*(1.-ERF_ALT(x))
2436                   arg=x-sq2*alnsign(m,n)
2437                   fs(m,n)=0.5*(1.-ERF_ALT(arg))
2438                   arg=x-1.5*sq2*alnsign(m,n)
2439                   fm(m,n)=0.5*(1.-ERF_ALT(arg))
2440                endif
2441
2442 !                     fn(m,n)=1. ! test
2443 !                     fs(m,n)=1.
2444 !                     fm(m,n)=1.
2445                 if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then
2446                    fluxn(m,n)=fn(m,n)*w
2447                    fluxs(m,n)=fs(m,n)*w
2448                    fluxm(m,n)=fm(m,n)*w
2449                 else
2450                    fluxn(m,n)=0
2451                    fluxs(m,n)=0
2452                    fluxm(m,n)=0
2453                endif
2454
2455 55000       continue      ! m=1,nsize_aer(n)
2456 55002       continue      ! n=1,ntype_aer
2457
2458             if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then
2459                flux_fullact = w
2460             else
2461                flux_fullact = 0.0
2462             endif
2463
2464 ! 04-nov-2005 rce - moved the code for "wnuc.le.0" from here
2465 ! to near the start the uniform undraft section
2466
2467 !.......................................................................
2468 !
2469 !   end calc. of activation fractions/fluxes
2470 !   for (single) uniform updraft (end of section 3)
2471 !
2472 !.......................................................................
2473
2474
2475
2476 60000 continue
2477
2478
2479 !            do n=1,ntype_aer
2480 !            do m=1,nsize_aer(n)
2481 !                write(6,'(a,2i3,5e10.1)')'n,m,na,wbar,sigw,fn,fm=',n,m,na(m,n),wbar,sigw,fn(m,n),fm(m,n)
2482 !            end do
2483 !            end do
2484
2485
2486       return
2487       end subroutine activate
2488
2489
2490
2491 !----------------------------------------------------------------------
2492 !----------------------------------------------------------------------
2493       subroutine maxsat(zeta,eta, &
2494                         maxd_atype,ntype_aer,maxd_asize,nsize_aer, &
2495                         sm,alnsign,f1,smax)
2496
2497 !      Calculates maximum supersaturation for multiple competing aerosol
2498 !      modes. Note that maxsat_init must be called before calling this
2499 !      subroutine.
2500
2501 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
2502 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
2503
2504       implicit none
2505
2506       integer, intent(in) :: maxd_atype
2507       integer, intent(in) :: ntype_aer
2508       integer, intent(in) :: maxd_asize
2509       integer, intent(in) :: nsize_aer(maxd_atype) ! number of size bins
2510       real, intent(in) :: sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
2511       real, intent(in) :: zeta, eta(maxd_asize,maxd_atype)
2512       real, intent(in) :: alnsign(maxd_asize,maxd_atype) ! ln(sigma)
2513       real, intent(in) :: f1(maxd_asize,maxd_atype)
2514       real, intent(out) :: smax ! maximum supersaturation
2515
2516       real :: g1, g2
2517       real thesum
2518       integer m ! size index
2519       integer n ! type index
2520
2521       do n=1,ntype_aer
2522       do m=1,nsize_aer(n)
2523          if(zeta.gt.1.e5*eta(m,n) .or. &
2524               sm(m,n)*sm(m,n).gt.1.e5*eta(m,n))then
2525 !           weak forcing. essentially none activated
2526             smax=1.e-20
2527          else
2528 !           significant activation of this mode. calc activation all modes.
2529             go to 1
2530          endif
2531       end do
2532       end do
2533
2534       return
2535
2536   1   continue
2537
2538       thesum=0
2539       do n=1,ntype_aer
2540       do m=1,nsize_aer(n)
2541          if(eta(m,n).gt.1.e-20)then
2542             g1=sqrt(zeta/eta(m,n))
2543             g1=g1*g1*g1
2544             g2=sm(m,n)/sqrt(eta(m,n)+3*zeta)
2545             g2=sqrt(g2)
2546             g2=g2*g2*g2
2547             thesum=thesum + &
2548                  (f1(m,n)*g1+(1.+0.25*alnsign(m,n))*g2)/(sm(m,n)*sm(m,n))
2549          else
2550             thesum=1.e20
2551          endif
2552       end do
2553       end do
2554
2555       smax=1./sqrt(thesum)
2556
2557       return
2558       end subroutine maxsat
2559
2560
2561
2562 !----------------------------------------------------------------------
2563 !----------------------------------------------------------------------
2564       subroutine maxsat_init(maxd_atype, ntype_aer, &
2565            maxd_asize, nsize_aer, alnsign, f1)
2566
2567 !     Calculates the f1 paramter needed by maxsat.
2568
2569 !     Abdul-Razzak and Ghan, A parameterization of aerosol activation.
2570 !     2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
2571
2572       implicit none
2573
2574       integer, intent(in)  :: maxd_atype
2575       integer, intent(in)  :: ntype_aer ! number of aerosol types
2576       integer, intent(in)  :: maxd_asize
2577       integer, intent(in)  :: nsize_aer(maxd_atype) ! number of size bins
2578       real,    intent(in)  :: alnsign(maxd_asize,maxd_atype) ! ln(sigma)
2579       real,    intent(out) :: f1(maxd_asize,maxd_atype)
2580
2581       integer m ! size index
2582       integer n ! type index
2583
2584 !  calculate and save f1(sigma), assumes sigma is invariant
2585 !  between calls to this init routine
2586
2587       do n=1,ntype_aer
2588          do m=1,nsize_aer(n)
2589             f1(m,n)=0.5*exp(2.5*alnsign(m,n)*alnsign(m,n))
2590          end do
2591       end do
2592
2593       end subroutine maxsat_init
2594
2595
2596
2597 !----------------------------------------------------------------------
2598 !----------------------------------------------------------------------
2599 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3);
2600 !     grid_id, ktau, i, j, isize, itype added to arg list to assist debugging
2601        subroutine loadaer(chem,k,kmn,kmx,num_chem,cs,npv, &
2602                           dlo_sect,dhi_sect,maxd_acomp, ncomp,                &
2603                           grid_id, ktau, i, j, isize, itype,   &
2604                           numptr_aer, numptrcw_aer, dens_aer,   &
2605                           massptr_aer, massptrcw_aer,   &
2606                           maerosol, maerosolcw,                 &
2607                           maerosol_tot, maerosol_totcw,         &
2608                           naerosol, naerosolcw,                 &
2609                           vaerosol, vaerosolcw)
2610
2611       implicit none
2612
2613 !      load aerosol number, surface, mass concentrations
2614
2615 !      input
2616
2617        integer, intent(in) ::  num_chem ! maximum number of consituents
2618        integer, intent(in) ::  k,kmn,kmx
2619        real,    intent(in) ::  chem(kmn:kmx,num_chem) ! aerosol mass, number mixing ratios
2620        real,    intent(in) ::  cs  ! air density (kg/m3)
2621        real,    intent(in) ::  npv ! number per volume concentration (/m3)
2622        integer, intent(in) ::  maxd_acomp,ncomp
2623        integer, intent(in) ::  numptr_aer,numptrcw_aer
2624        integer, intent(in) ::  massptr_aer(maxd_acomp), massptrcw_aer(maxd_acomp)
2625        real,    intent(in) ::  dens_aer(maxd_acomp) ! aerosol material density (g/cm3)
2626        real,    intent(in) ::  dlo_sect,dhi_sect ! minimum, maximum diameter of section (cm)
2627        integer, intent(in) ::  grid_id, ktau, i, j, isize, itype
2628
2629 !      output
2630
2631        real, intent(out) ::  naerosol                ! interstitial number conc (/m3)
2632        real, intent(out) ::  naerosolcw              ! activated    number conc (/m3)
2633        real, intent(out) ::  maerosol(maxd_acomp)   ! interstitial mass conc (kg/m3)
2634        real, intent(out) ::  maerosolcw(maxd_acomp) ! activated    mass conc (kg/m3)
2635        real, intent(out) ::  maerosol_tot   ! total-over-species interstitial mass conc (kg/m3)
2636        real, intent(out) ::  maerosol_totcw ! total-over-species activated    mass conc (kg/m3)
2637        real, intent(out) ::  vaerosol       ! interstitial volume conc (m3/m3)
2638        real, intent(out) ::  vaerosolcw     ! activated volume conc (m3/m3)
2639
2640 !      internal
2641
2642        integer lnum,lnumcw,l,ltype,lmass,lmasscw,lsfc,lsfccw
2643        real num_at_dhi, num_at_dlo
2644        real npv_at_dhi, npv_at_dlo
2645        real, parameter :: pi = 3.1415926526
2646        real specvol ! inverse aerosol material density (m3/kg)
2647
2648           lnum=numptr_aer
2649           lnumcw=numptrcw_aer
2650           maerosol_tot=0.
2651           maerosol_totcw=0.
2652           vaerosol=0.
2653           vaerosolcw=0.
2654           do l=1,ncomp
2655              lmass=massptr_aer(l)
2656              lmasscw=massptrcw_aer(l)
2657              maerosol(l)=chem(k,lmass)*cs
2658              maerosol(l)=max(maerosol(l),0.)
2659              maerosolcw(l)=chem(k,lmasscw)*cs
2660              maerosolcw(l)=max(maerosolcw(l),0.)
2661              maerosol_tot=maerosol_tot+maerosol(l)
2662              maerosol_totcw=maerosol_totcw+maerosolcw(l)
2663 ! [ 1.e-3 factor because dens_aer is (g/cm3), specvol is (m3/kg) ]
2664              specvol=1.0e-3/dens_aer(l)
2665              vaerosol=vaerosol+maerosol(l)*specvol
2666              vaerosolcw=vaerosolcw+maerosolcw(l)*specvol
2667 !            write(6,'(a,3e12.2)')'maerosol,dens_aer,vaerosol=',maerosol(l),dens_aer(l),vaerosol
2668           enddo
2669
2670           if(lnum.gt.0)then
2671 !            aerosol number predicted
2672 ! [ 1.0e6 factor because because dhi_ & dlo_sect are (cm), vaerosol is (m3) ]
2673              npv_at_dhi = 6.0e6/(pi*dhi_sect*dhi_sect*dhi_sect)
2674              npv_at_dlo = 6.0e6/(pi*dlo_sect*dlo_sect*dlo_sect)
2675
2676              naerosol=chem(k,lnum)*cs
2677              naerosolcw=chem(k,lnumcw)*cs
2678              num_at_dhi = vaerosol*npv_at_dhi
2679              num_at_dlo = vaerosol*npv_at_dlo
2680              naerosol = max( num_at_dhi, min( num_at_dlo, naerosol ) )
2681
2682 !            write(6,'(a,5e10.1)')'naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect', &
2683 !                          naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect
2684              num_at_dhi = vaerosolcw*npv_at_dhi
2685              num_at_dlo = vaerosolcw*npv_at_dlo
2686              naerosolcw = max( num_at_dhi, min( num_at_dlo, naerosolcw ) )
2687           else
2688 !            aerosol number diagnosed from mass and prescribed size
2689              naerosol=vaerosol*npv
2690              naerosol=max(naerosol,0.)
2691              naerosolcw=vaerosolcw*npv
2692              naerosolcw=max(naerosolcw,0.)
2693           endif
2694
2695
2696        return
2697        end subroutine loadaer
2698
2699
2700
2701 !-----------------------------------------------------------------------
2702         real function erfc_num_recipes( x )
2703 !
2704 !   from press et al, numerical recipes, 1990, page 164
2705 !
2706         implicit none
2707         real x
2708         double precision erfc_dbl, dum, t, zz
2709
2710         zz = abs(x)
2711         t = 1.0/(1.0 + 0.5*zz)
2712
2713 !       erfc_num_recipes =
2714 !     &   t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
2715 !     &   t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
2716 !     &                                    t*(-1.13520398 +
2717 !     &   t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
2718
2719         dum =  ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
2720           t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
2721                                            t*(-1.13520398 +   &
2722           t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
2723
2724         erfc_dbl = t * exp(dum)
2725         if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
2726
2727         erfc_num_recipes = erfc_dbl
2728
2729         return
2730         end function erfc_num_recipes
2731
2732 !-----------------------------------------------------------------------
2733     real function erf_alt( x )
2734
2735     implicit none
2736
2737     real,intent(in) :: x
2738
2739     erf_alt = 1. - erfc_num_recipes(x)
2740
2741     end function erf_alt
2742
2743 END MODULE module_mixactivate