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