wrfv2_fire/phys/module_mixactivate.F

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