phys/module_cam_mp_microp_aero.F

   1 #define WRF_PORT
   2 #if ( WRF_CHEM == 1 )
   3 #       include "../chem/MODAL_AERO_CPP_DEFINES.h"
   4 #else
   5 #       define MODAL_AERO
   6 #       define MODAL_AERO_3MODE
   7 #endif
   8 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   9 module microp_aero
  10
  11 !---------------------------------------------------------------------------------
  12 ! Purpose:
  13 !   CAM Interface for aerosol activation
  14 !
  15 ! Author: Andrew Gettelman
  16 ! Based on code from: Hugh Morrison, Xiaohong Liu and Steve Ghan
  17 ! May 2010
  18 ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
  19 !                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
  20 ! for questions contact Andrew Gettelman  (andrew@ucar.edu)
  21 ! Modifications: A. Gettelman Nov 2010  - changed to support separation of
  22 !                microphysics and macrophysics and concentrate aerosol information here
  23 !
  24 !---------------------------------------------------------------------------------
  25
  26   use shr_kind_mod,  only: r8=>shr_kind_r8
  27 #ifndef WRF_PORT
  28   use spmd_utils,    only: masterproc
  29   use ppgrid,        only: pcols, pver, pverp
  30 #else
  31   use module_cam_support, only: masterproc, pcols, pver, pverp
  32 #endif
  33   use physconst,     only: gravit, rair, tmelt, cpair, rh2o, r_universal, mwh2o, rhoh2o
  34   use physconst,     only: latvap, latice
  35 #ifndef WRF_PORT
  36   use abortutils,    only: endrun
  37 #else
  38    use module_cam_support, only: endrun
  39 #endif
  40   use error_function, only: erf,erfc
  41   use wv_saturation,  only: estblf, hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice, &
  42                             vqsatd2, vqsatd2_single,polysvp
  43 #ifndef WRF_PORT
  44   use cam_history,    only: addfld, add_default, phys_decomp, outfld
  45   use cam_logfile,    only: iulog
  46   use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_props
  47   use phys_control,   only: phys_getopts
  48   use cldwat2m_macro, only: rhmini, rhmaxi
  49 #else
  50   use module_cam_support, only: addfld, add_default, phys_decomp, outfld, iulog
  51 #endif
  52
  53
  54  implicit none
  55  private
  56  save
  57
  58  public :: ini_microp_aero, microp_aero_ts
  59
  60 ! Private module data
  61
  62    logical             :: wsubTKE  ! If .true. (.false.), use UW's TKE (kkvh) for computing wsub
  63
  64
  65       real(r8), private:: rn_dst1, rn_dst2, rn_dst3, rn_dst4 !dust number mean radius for contact freezing
  66       real(r8), private:: pi    ! pi
  67       real(r8), private:: tt0   ! Freezing temperature
  68       real(r8), private:: qsmall !min mixing ratio
  69
  70   real(r8), private::  r              !Dry air Gas constant
  71
  72
  73 ! activate parameters
  74
  75       integer, private:: psat
  76       parameter (psat=6) ! number of supersaturations to calc ccn concentration
  77       real(r8), private:: aten
  78
  79       real(r8), allocatable, private:: alogsig(:) ! natl log of geometric standard dev of aerosol
  80       real(r8), allocatable, private:: exp45logsig(:)
  81       real(r8), allocatable, private:: argfactor(:)
  82       real(r8), allocatable, private:: amcube(:) ! cube of dry mode radius (m)
  83       real(r8), allocatable, private:: smcrit(:) ! critical supersatuation for activation
  84       real(r8), allocatable, private:: lnsm(:) ! ln(smcrit)
  85       real(r8), private:: amcubesulfate(pcols) ! cube of dry mode radius (m) of sulfate
  86       real(r8), private:: smcritsulfate(pcols) ! critical supersatuation for activation of sulfate
  87       real(r8), allocatable, private:: amcubefactor(:) ! factors for calculating mode radius
  88       real(r8), allocatable, private:: smcritfactor(:) ! factors for calculating critical supersaturation
  89       real(r8), private:: super(psat)
  90       real(r8), private:: alogten,alog2,alog3,alogaten
  91       real(r8), private, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
  92                (/0.02,0.05,0.1,0.2,0.5,1.0/)
  93       real(r8), allocatable, private:: ccnfact(:,:)
  94
  95       real(r8), allocatable, private:: f1(:),f2(:) ! abdul-razzak functions of width
  96       real(r8), private:: third, sixth,zero
  97       real(r8), private:: sq2, sqpi
  98
  99       integer :: naer_all    ! number of aerosols affecting climate
 100       integer :: idxsul = -1 ! index in aerosol list for sulfate
 101       integer :: idxdst1 = -1 ! index in aerosol list for dust1
 102       integer :: idxdst2 = -1 ! index in aerosol list for dust2
 103       integer :: idxdst3 = -1 ! index in aerosol list for dust3
 104       integer :: idxdst4 = -1 ! index in aerosol list for dust4
 105
 106       integer :: idxbcphi = -1 ! index in aerosol list for Soot (BCPHIL)
 107
 108       ! aerosol properties
 109       character(len=20), allocatable :: aername(:)
 110       real(r8), allocatable :: dryrad_aer(:)
 111       real(r8), allocatable :: density_aer(:)
 112       real(r8), allocatable :: hygro_aer(:)
 113       real(r8), allocatable :: dispersion_aer(:)
 114       real(r8), allocatable :: num_to_mass_aer(:)
 115
 116 contains
 117
 118 !===============================================================================
 119
 120 subroutine ini_microp_aero
 121
 122 !-----------------------------------------------------------------------
 123 !
 124 ! Purpose:
 125 ! initialize constants for aerosols needed by microphysics
 126 ! called from stratiform.F90
 127 !
 128 ! Author: Andrew Gettelman May 2010
 129 !
 130 !-----------------------------------------------------------------------
 131 #ifdef MODAL_AERO
 132     use ndrop,           only: activate_init
 133     use constituents,    only: cnst_name
 134 #ifndef WRF_PORT
 135     use cam_history,     only: fieldname_len
 136     use spmd_utils,      only: masterproc
 137     use modal_aero_data, only: cnst_name_cw, &
 138                                lmassptr_amode, lmassptrcw_amode, &
 139                                nspec_amode, ntot_amode, numptr_amode, numptrcw_amode, ntot_amode
 140 #else
 141     use module_cam_support, only: fieldname_len, masterproc
 142     use modal_aero_data,    only: cnst_name_cw => cnst_name_cw_mp, &
 143          lmassptr_amode => lmassptr_amode_mp, lmassptrcw_amode => lmassptrcw_amode_mp , &
 144          nspec_amode    => nspec_amode_mp   , ntot_amode       => ntot_amode_mp,        &
 145          numptr_amode   => numptr_amode_mp  , numptrcw_amode   => numptrcw_amode_mp
 146 #endif
 147
 148     integer                        :: lphase, lspec
 149     character(len=fieldname_len)   :: tmpname
 150     character(len=fieldname_len+3) :: fieldname
 151     character(128)                 :: long_name
 152     character(8)                   :: unit
 153     logical                        :: history_aerosol      ! Output the MAM aerosol tendencies
 154
 155 #endif
 156
 157    integer k
 158
 159    integer l,m, iaer
 160    real(r8) surften       ! surface tension of water w/respect to air (N/m)
 161    real(r8) arg
 162    real(r8) derf
 163
 164    character(len=16) :: eddy_scheme = ' '
 165    logical           :: history_microphysics     ! output variables for microphysics diagnostics package
 166 #ifndef WRF_PORT
 167    ! Query the PBL eddy scheme
 168    call phys_getopts(eddy_scheme_out              = eddy_scheme,           &
 169                      history_microphysics_out     = history_microphysics   )
 170    wsubTKE = .false.
 171    if (trim(eddy_scheme) .eq. 'diag_TKE') wsubTKE = .true.
 172 #else
 173    history_microphysics =.FALSE.
 174    wsubTKE              =.TRUE.
 175 #endif
 176
 177    ! Access the physical properties of the aerosols that are affecting the climate
 178    ! by using routines from the rad_constituents module.
 179 #ifndef WRF_PORT
 180    call rad_cnst_get_info(0, naero=naer_all)
 181    allocate( &
 182       aername(naer_all),        &
 183       dryrad_aer(naer_all),     &
 184       density_aer(naer_all),    &
 185       hygro_aer(naer_all),      &
 186       dispersion_aer(naer_all), &
 187       num_to_mass_aer(naer_all) )
 188
 189    do iaer = 1, naer_all
 190       call rad_cnst_get_aer_props(0, iaer, &
 191          aername         = aername(iaer), &
 192          dryrad_aer      = dryrad_aer(iaer), &
 193          density_aer     = density_aer(iaer), &
 194          hygro_aer       = hygro_aer(iaer), &
 195          dispersion_aer  = dispersion_aer(iaer), &
 196          num_to_mass_aer = num_to_mass_aer(iaer) )
 197
 198
 199       ! Look for sulfate aerosol in this list (Bulk aerosol only)
 200       if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
 201       if (trim(aername(iaer)) == 'DUST1') idxdst1 = iaer
 202       if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
 203       if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
 204       if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
 205       if (trim(aername(iaer)) == 'BCPHIL') idxbcphi = iaer
 206
 207
 208 !addfield calls for outputting aerosol number concentration
 209
 210       call addfld(aername(iaer)//'_m3', 'm-3', pver, 'A', 'aerosol number concentration', phys_decomp)
 211
 212    end do
 213
 214    if (masterproc) then
 215       write(iulog,*) 'ini_micro: iaer, name, dryrad, density, hygro, dispersion, num_to_mass'
 216 #ifdef WRF_PORT
 217         call wrf_message(iulog)
 218 #endif
 219       do iaer = 1, naer_all
 220          write(iulog,*) iaer, aername(iaer), dryrad_aer(iaer), density_aer(iaer), hygro_aer(iaer), &
 221             dispersion_aer(iaer), num_to_mass_aer(iaer)
 222 #ifdef WRF_PORT
 223          call wrf_message(iulog)
 224 #endif
 225       end do
 226       if (idxsul < 1) then
 227          write(iulog,*) 'ini_micro: SULFATE aerosol properties NOT FOUND'
 228 #ifdef WRF_PORT
 229          call wrf_message(iulog)
 230 #endif
 231       else
 232          write(iulog,*) 'ini_micro: SULFATE aerosol properties FOUND at index ',idxsul
 233 #ifdef WRF_PORT
 234          call wrf_message(iulog)
 235 #endif
 236       end if
 237    end if
 238 #endif
 239    call addfld ('CCN1    ','#/cm3   ',pver, 'A','CCN concentration at S=0.02%',phys_decomp)
 240    call addfld ('CCN2    ','#/cm3   ',pver, 'A','CCN concentration at S=0.05%',phys_decomp)
 241    call addfld ('CCN3    ','#/cm3   ',pver, 'A','CCN concentration at S=0.1%',phys_decomp)
 242    call addfld ('CCN4    ','#/cm3   ',pver, 'A','CCN concentration at S=0.2%',phys_decomp)
 243    call addfld ('CCN5    ','#/cm3   ',pver, 'A','CCN concentration at S=0.5%',phys_decomp)
 244    call addfld ('CCN6    ','#/cm3   ',pver, 'A','CCN concentration at S=1.0%',phys_decomp)
 245
 246    call add_default('CCN3',  1, ' ' )
 247
 248    call addfld ('NIHF    ','#/m3   ',pver, 'A','Activated Ice Number Concentation due to homogenous freezing',phys_decomp)
 249    call addfld ('NIDEP    ','#/m3   ',pver, 'A','Activated Ice Number Concentation due to deposition nucleation',phys_decomp)
 250    call addfld ('NIIMM    ','#/m3   ',pver, 'A','Activated Ice Number Concentation due to immersion freezing',phys_decomp)
 251    call addfld ('NIMEY    ','#/m3   ',pver, 'A','Activated Ice Number Concentation due to meyers deposition',phys_decomp)
 252
 253 ! contact freezing due to dust
 254 ! dust number mean radius (m), Zender et al JGR 2003 assuming number mode radius of 0.6 micron, sigma=2
 255
 256    rn_dst1=0.258e-6_r8
 257    rn_dst2=0.717e-6_r8
 258    rn_dst3=1.576e-6_r8
 259    rn_dst4=3.026e-6_r8
 260
 261 ! smallest mixing ratio considered in microphysics
 262
 263    qsmall = 1.e-18_r8
 264
 265 ! set parameters for droplet activation, following abdul-razzak and ghan 2000, JGR
 266
 267 !      mathematical constants
 268
 269       zero=0._r8
 270       third=1./3._r8
 271       sixth=1./6._r8
 272       sq2=sqrt(2._r8)
 273       pi=4._r8*atan(1.0_r8)
 274       sqpi=sqrt(pi)
 275
 276       r= rair        !Dry air Gas constant: note units(phys_constants are in J/K/kmol)
 277
 278 ! freezing temperature
 279       tt0=273.15_r8  !should be tmelt
 280
 281       surften=0.076_r8
 282       aten=2.*mwh2o*surften/(r_universal*tt0*rhoh2o)
 283       alogaten=log(aten)
 284       alog2=log(2._r8)
 285       alog3=log(3._r8)
 286       super(:)=0.01*supersat(:)
 287 #ifndef WRF_PORT
 288       allocate(alogsig(naer_all),      &
 289                exp45logsig(naer_all),  &
 290                argfactor(naer_all),    &
 291                amcube(naer_all),       &
 292                smcrit(naer_all),       &
 293                lnsm(naer_all),         &
 294                amcubefactor(naer_all), &
 295                smcritfactor(naer_all), &
 296                ccnfact(psat,naer_all), &
 297                f1(naer_all),           &
 298                f2(naer_all)            )
 299
 300       do m=1,naer_all
 301           alogsig(m)=log(dispersion_aer(m))
 302           exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
 303           argfactor(m)=2./(3.*sqrt(2.)*alogsig(m))
 304           f1(m)=0.5*exp(2.5*alogsig(m)*alogsig(m))
 305           f2(m)=1.+0.25*alogsig(m)
 306           amcubefactor(m)=3._r8/(4._r8*pi*exp45logsig(m)*density_aer(m))
 307           smcritfactor(m)=2._r8*aten*sqrt(aten/(27._r8*max(1.e-10_r8,hygro_aer(m))))
 308           amcube(m)=amcubefactor(m)/(num_to_mass_aer(m))
 309           if(hygro_aer(m).gt.1.e-10) then
 310              smcrit(m)=smcritfactor(m)/sqrt(amcube(m))
 311           else
 312              smcrit(m)=100.
 313           endif
 314           lnsm(m)=log(smcrit(m))
 315
 316           do l=1,psat
 317              arg=argfactor(m)*log(smcrit(m)/super(l))
 318              if(arg<2)then
 319                 if(arg<-2)then
 320                    ccnfact(l,m)=1.e-6
 321                 else
 322                    ccnfact(l,m)=1.e-6_r8*0.5_r8*erfc(arg)
 323                 endif
 324              else
 325                 ccnfact(l,m) = 0.
 326              endif
 327           enddo
 328       end do
 329
 330 #ifdef MODAL_AERO
 331  call phys_getopts( history_aerosol_out        = history_aerosol)
 332  call activate_init
 333
 334 ! Add dropmixnuc tendencies for all modal aerosol species
 335     do m = 1, ntot_amode
 336     do lphase = 1, 2
 337     do lspec = 0, nspec_amode(m)+1   ! loop over number + chem constituents + water
 338        unit = 'kg/m2/s'
 339        if (lspec == 0) then   ! number
 340           unit = '#/m2/s'
 341           if (lphase == 1) then
 342              l = numptr_amode(m)
 343           else
 344              l = numptrcw_amode(m)
 345           endif
 346        else if (lspec <= nspec_amode(m)) then   ! non-water mass
 347           if (lphase == 1) then
 348              l = lmassptr_amode(lspec,m)
 349           else
 350              l = lmassptrcw_amode(lspec,m)
 351           endif
 352        else   ! water mass
 353           cycle
 354        end if
 355        if (lphase == 1) then
 356           tmpname = cnst_name(l)
 357        else
 358           tmpname = cnst_name_cw(l)
 359        end if
 360
 361        fieldname = trim(tmpname) // '_mixnuc1'
 362        long_name = trim(tmpname) // ' dropmixnuc mixnuc column tendency'
 363        call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
 364        if ( history_aerosol ) then
 365           call add_default( fieldname, 1, ' ' )
 366           if ( masterproc ) write(*,'(2a)') 'microp_driver_init addfld - ', fieldname
 367        endif
 368
 369     end do   ! lspec
 370     end do   ! lphase
 371     end do   ! m
 372 #endif
 373 #endif
 374  return
 375  end subroutine ini_microp_aero
 376
 377
 378 !===============================================================================
 379 !activation routine for each timestep goes here...
 380
 381       subroutine microp_aero_ts (   &
 382    lchnk, ncol, deltatin, t, ttend,        &
 383    qn, qc, qi,               &
 384    nc, ni, p, pdel, cldn,                   &
 385    liqcldf, icecldf,                        &
 386    cldo, pint, rpdel, zm, omega,            &
 387 #ifdef MODAL_AERO
 388    qaer, cflx, qaertend, dgnumwet,dgnum, &
 389 #else
 390    aer_mmr,                                 &
 391 #endif
 392    kkvh, tke, turbtype, smaw, wsub, wsubi, &
 393    naai, naai_hom, npccn, rndst, nacon     &
 394 #ifdef WRF_PORT
 395    ,qqcw                                   &
 396 #endif
 397                                            )
 398
 399    use wv_saturation, only: vqsatd, vqsatd_water
 400 #ifndef WRF_PORT
 401    use constituents,  only: pcnst
 402 #else
 403    use module_cam_support, only: pcnst => pcnst_mp
 404 #endif
 405 #ifdef MODAL_AERO
 406    use ndrop,         only: dropmixnuc
 407 #ifndef WRF_PORT
 408    use modal_aero_data, only:numptr_amode, modeptr_accum, modeptr_coarse, modeptr_aitken, &
 409                               lptr_dust_a_amode,lptr_nacl_a_amode,ntot_amode
 410 #else
 411    use modal_aero_data, only:numptr_amode => numptr_amode_mp, modeptr_accum => modeptr_accum_mp, &
 412         modeptr_coarse    => modeptr_coarse_mp   , modeptr_aitken    => modeptr_aitken_mp, &
 413         lptr_dust_a_amode => lptr_dust_a_amode_mp, lptr_nacl_a_amode => lptr_nacl_a_amode_mp, &
 414         ntot_amode        => ntot_amode_mp, modeptr_coardust => modeptr_coardust_mp  !BSINGH - Added modeptr_coarsedust for MAM7 compliance
 415
 416 #endif
 417 #endif
 418
 419    ! input arguments
 420    integer,  intent(in) :: lchnk
 421    integer,  intent(in) :: ncol
 422    real(r8), intent(in) :: deltatin             ! time step (s)
 423    real(r8), intent(in) :: t(pcols,pver)       ! input temperature (K)
 424    real(r8), intent(in) :: ttend(pcols,pver)    ! non-microphysical temperature tendency (K/s)
 425    real(r8), intent(in) :: qn(pcols,pver)       ! input h20 vapor mixing ratio (kg/kg)
 426    ! note: all input cloud variables are grid-averaged
 427    real(r8), intent(in) :: qc(pcols,pver)    ! cloud water mixing ratio (kg/kg)
 428    real(r8), intent(in) :: qi(pcols,pver)    ! cloud ice mixing ratio (kg/kg)
 429    real(r8), intent(in) :: nc(pcols,pver)    ! cloud water number conc (1/kg)
 430    real(r8), intent(in) :: ni(pcols,pver)    ! cloud ice number conc (1/kg)
 431    real(r8), intent(in) :: p(pcols,pver)        ! air pressure (pa)
 432    real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across level (pa)
 433    real(r8), intent(in) :: cldn(pcols,pver)     ! cloud fraction
 434    real(r8), intent(in) :: icecldf(pcols,pver)  ! ice cloud fraction
 435    real(r8), intent(in) :: liqcldf(pcols,pver)  ! liquid cloud fraction
 436    real(r8), intent(in) :: cldo(pcols,pver)  ! old cloud fraction
 437    real(r8), intent(in) :: pint(pcols,pverp)    ! air pressure layer interfaces (pa)
 438    real(r8), intent(in) :: rpdel(pcols,pver)    ! inverse pressure difference across level (pa)
 439    real(r8), intent(in) :: zm(pcols,pver)       ! geopotential height of model levels (m)
 440    real(r8), intent(in) :: omega(pcols,pver)    ! vertical velocity (Pa/s)
 441 #ifdef MODAL_AERO
 442    real(r8), intent(in) :: qaer(pcols,pver,pcnst) ! aerosol number and mass mixing ratios
 443    real(r8), intent(in) :: cflx(pcols,pcnst)      ! constituent flux from surface
 444    real(r8), intent(inout) :: qaertend(pcols,pver,pcnst) ! qaer tendency (1/s)
 445 #ifdef WRF_PORT
 446    real(r8), intent(inout) :: qqcw(pcols,pver,pcnst) ! cloud-borne aerosol
 447 #endif
 448    real(r8), intent(in) :: dgnumwet(pcols,pver,ntot_amode) ! aerosol mode diameter
 449    real(r8), intent(in) :: dgnum(pcols,pver,ntot_amode) ! aerosol mode dry diameter
 450 #else
 451    real(r8), intent(in) :: aer_mmr(:,:,:)       ! aerosol mass mixing ratio
 452 #endif
 453    real(r8), intent(in) :: kkvh(pcols,pver+1) ! vertical eddy diff coef (m2 s-1)
 454    real(r8), intent(in) :: tke(pcols,pver+1)    ! TKE from the UW PBL scheme (m2 s-2)
 455    real(r8), intent(in) :: turbtype(pcols,pver+1) ! Turbulence type at each interface
 456    real(r8), intent(in) :: smaw(pcols,pver+1)     ! Normalized instability function of momentum. 0 <=  <= 4.964 and 1 at neutral condition.
 457
 458    ! output arguments
 459
 460    real(r8), intent(out) :: wsub(pcols,pver)    ! diagnosed sub-grid vertical velocity st. dev. (m/s)
 461    real(r8), intent(out) :: wsubi(pcols,pver)   ! diagnosed sub-grid vertical velocity ice (m/s)
 462    real(r8), intent(out) :: naai(pcols,pver)    ! number of activated aerosol for ice nucleation
 463    real(r8), intent(out) :: naai_hom(pcols,pver)! number of activated aerosol for ice nucleation (homogeneous freezing only)
 464    real(r8), intent(out) :: npccn(pcols,pver)   ! number of CCN (liquid activated)
 465    real(r8), intent(out) :: rndst(pcols,pver,4) ! radius of 4 dust bins for contact freezing (from microp_aero_ts)
 466    real(r8), intent(out) :: nacon(pcols,pver,4) ! number in 4 dust bins for contact freezing  (from microp_aero_ts)
 467
 468
 469 ! local workspace
 470 ! all units mks unless otherwise stated
 471
 472    real(r8) :: tkem(pcols,pver)     ! Layer-mean TKE
 473    real(r8) :: smm(pcols,pver)      ! Layer-mean instability function
 474    real(r8) :: relhum(pcols,pver) ! relative humidity
 475    real(r8) :: cldm(pcols,pver)   ! cloud fraction
 476    real(r8) :: icldm(pcols,pver)   ! ice cloud fraction
 477    real(r8) :: lcldm(pcols,pver)   ! liq cloud fraction
 478    real(r8) :: nfice(pcols,pver) !fice variable
 479    real(r8) :: dumfice
 480    real(r8) :: qcld          ! total cloud water
 481         real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
 482         real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
 483         real(r8) :: nctend_mixnuc(pcols,pver)
 484         real(r8) :: deltat        ! sub-time step (s)
 485         real(r8) :: nihf2,niimm2,nidep2,nimey2,dum2 !dummy variables for activated ice nuclei by diffferent processes
 486         real(r8) arg  !variable for CCN number (BAM)
 487
 488          integer ftrue  ! integer used to determine cloud depth
 489          real(r8) :: dum        ! temporary dummy variable
 490
 491         real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg)
 492 !       real(r8) :: t(pcols,pver) ! temperature (K)
 493         real(r8) :: ncloc(pcols,pver) ! local variable for nc
 494         real(r8) :: rho(pcols,pver) ! air density (kg m-3)
 495         real(r8) :: mincld  ! minimum allowed cloud fraction
 496
 497 ! used in contact freezing via dust particles
 498         real(r8)  tcnt, viscosity, mfp
 499         real(r8)  slip1, slip2, slip3, slip4
 500         real(r8)  dfaer1, dfaer2, dfaer3, dfaer4
 501         real(r8)  nacon1,nacon2,nacon3,nacon4
 502
 503         real(r8) dmc,ssmc,dstrn  ! variables for modal scheme.
 504
 505 ! returns from function/subroutine calls
 506         real(r8) :: tsp(pcols,pver)      ! saturation temp (K)
 507         real(r8) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
 508         real(r8) :: qsphy(pcols,pver)      ! saturation mixing ratio (kg/kg): hybrid rh
 509         real(r8) :: qs(pcols)            ! liquid-ice weighted sat mixing rat (kg/kg)
 510         real(r8) :: es(pcols)            ! liquid-ice weighted sat vapor press (pa)
 511         real(r8) :: esl(pcols,pver)      ! liquid sat vapor pressure (pa)
 512         real(r8) :: esi(pcols,pver)      ! ice sat vapor pressure (pa)
 513         real(r8) :: gammas(pcols)        ! parameter for cond/evap of cloud water
 514
 515 ! new aerosol variables
 516         real(r8), allocatable :: naermod(:) ! aerosol number concentration (/m3)
 517         real(r8), allocatable :: naer2(:,:,:)   ! new aerosol number concentration (/m3)
 518
 519         real(r8), allocatable :: maerosol(:,:)   ! aerosol mass conc (kg/m3)
 520         real(r8) naer(pcols)
 521         real(r8) ccn(pcols,pver,psat)        ! number conc of aerosols activated at supersat
 522         character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
 523
 524 !output for ice nucleation
 525           real(r8) :: nimey(pcols,pver) !output number conc of ice nuclei due to meyers deposition (1/m3)
 526           real(r8) :: nihf(pcols,pver)  !output number conc of ice nuclei due to heterogenous freezing (1/m3)
 527           real(r8) :: nidep(pcols,pver) !output number conc of ice nuclei due to deoposion nucleation (hetero nuc) (1/m3)
 528           real(r8) :: niimm(pcols,pver) !output number conc of ice nuclei due to immersion freezing (hetero nuc) (1/m3)
 529
 530 ! loop array variables
 531          integer i,k,nstep,n, l
 532          integer ii,kk, m
 533          integer, allocatable :: ntype(:)
 534
 535 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 536
 537 ! initialize  output fields for ice nucleation
 538     nimey(1:ncol,1:pver)=0._r8
 539     nihf(1:ncol,1:pver)=0._r8
 540     nidep(1:ncol,1:pver)=0._r8
 541     niimm(1:ncol,1:pver)=0._r8
 542
 543 ! initialize output
 544     naai(1:ncol,1:pver)=0._r8
 545     naai_hom(1:ncol,1:pver)=0._r8
 546     npccn(1:ncol,1:pver)=0._r8
 547     nacon(1:ncol,1:pver,:)=0._r8
 548
 549 ! set default or fixed dust bins for contact freezing
 550     rndst(1:ncol,1:pver,1)=rn_dst1
 551     rndst(1:ncol,1:pver,2)=rn_dst2
 552     rndst(1:ncol,1:pver,3)=rn_dst3
 553     rndst(1:ncol,1:pver,4)=rn_dst4
 554
 555 ! TKE initialization
 556     tkem(1:ncol,1:pver)=0._r8
 557     smm(1:ncol,1:pver)=0._r8
 558
 559     allocate(naermod(naer_all), &
 560       naer2(pcols,pver,naer_all), &
 561       maerosol(1,naer_all), &
 562       ntype(naer_all))
 563
 564 ! assign variable deltat for sub-stepping...
 565         deltat=deltatin
 566
 567 ! initialize multi-level fields
 568         q(1:ncol,1:pver)=qn(1:ncol,1:pver)
 569         ncloc(1:ncol,1:pver)=nc(1:ncol,1:pver)
 570
 571 ! parameters for scheme
 572         mincld=0.0001_r8
 573
 574 ! initialize time-varying parameters
 575
 576         do k=1,pver
 577            do i=1,ncol
 578               rho(i,k)=p(i,k)/(r*t(i,k))
 579            end do
 580         end do
 581
 582 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 583 ! More refined computation of sub-grid vertical velocity
 584 ! Set to be zero at the surface by initialization.
 585
 586         if( wsubTKE ) then
 587
 588             do i = 1, ncol
 589             do k = 1, pver
 590                tkem(i,k) = 0.5_r8 * ( tke(i,k)  +  tke(i,k+1) )
 591                smm(i,k)  = 0.5_r8 * ( smaw(i,k) + smaw(i,k+1) )
 592                if( turbtype(i,k) .eq. 3._r8 ) then       ! Bottom entrainment interface
 593                    tkem(i,k) = 0.5_r8 *  tke(i,k+1)
 594                    smm(i,k)  = 0.5_r8 * smaw(i,k+1)
 595                elseif( turbtype(i,k+1) .eq. 4._r8 ) then ! Top entrainment interface
 596                    tkem(i,k) = 0.5_r8 *  tke(i,k)
 597                    smm(i,k)  = 0.5_r8 * smaw(i,k)
 598                endif
 599                smm(i,k) = 0.259_r8*smm(i,k)
 600                smm(i,k) = max(smm(i,k), 0.4743_r8)
 601             end do
 602             end do
 603
 604         endif
 605
 606            do i=1,ncol
 607               ftrue=0
 608               do k=1,pver
 609 ! get sub-grid vertical velocity from diff coef.
 610 ! following morrison et al. 2005, JAS
 611 ! assume mixing length of 30 m
 612
 613               dum=(kkvh(i,k)+kkvh(i,k+1))/2._r8/30._r8
 614 ! use maximum sub-grid vertical vel of 10 m/s
 615               dum=min(dum,10._r8)
 616 ! set wsub to value at current vertical level
 617               wsub(i,k)=dum
 618             ! new computation
 619               if( wsubTKE ) then
 620                   wsub(i,k) = sqrt(0.5_r8*(tke(i,k)+tke(i,k+1))*(2._r8/3._r8))
 621                   wsub(i,k) = min(wsub(i,k),10._r8)
 622               endif
 623               wsubi(i,k) = max(0.001_r8,wsub(i,k))
 624               wsubi(i,k) = min(wsubi(i,k),0.2_r8)
 625               wsub(i,k)  = max(0.20_r8,wsub(i,k))
 626            end do
 627         end do
 628
 629 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 630 !Get humidity and saturation vapor pressures
 631
 632         do k=1,pver
 633
 634 ! find wet bulk temperature and saturation value for provisional t and q without
 635 ! condensation
 636
 637         call vqsatd_water(t(1,k),p(1,k),es,qs,gammas,ncol) ! use rhw
 638
 639         do i=1,ncol
 640
 641         esl(i,k)=polysvp(t(i,k),0)
 642         esi(i,k)=polysvp(t(i,k),1)
 643
 644 ! hm fix, make sure when above freezing that esi=esl, not active yet
 645         if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
 646
 647         relhum(i,k)=q(i,k)/qs(i)
 648
 649 ! get cloud fraction, check for minimum
 650
 651            cldm(i,k)=max(cldn(i,k),mincld)
 652
 653            icldm(i,k)=max(icecldf(i,k),mincld)
 654            lcldm(i,k)=max(liqcldf(i,k),mincld)
 655
 656
 657 ! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)
 658
 659         nfice(i,k)=0._r8
 660         dumfice=qc(i,k)+qi(i,k)
 661         if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then
 662            nfice(i,k)=qi(i,k)/dumfice
 663         endif
 664
 665 #ifndef MODAL_AERO
 666
 667         do m=1,naer_all
 668            ntype(m)=1
 669            maerosol(1,m)=aer_mmr(i,k,m)*rho(i,k)
 670
 671 ! For bulk aerosols only:
 672 !   set number nucleated for sulfate based on Lohmann et al 2000 (JGR) eq 2
 673 !           Na=340.*(massSO4)^0.58  where Na=cm-3 and massSO4=ug/m3
 674 !   convert units to Na [m-3] and SO4 [kgm-3]
 675 !   Na(m-3)= 1.e6 cm3 m-3 Na(cm-3)=340. *(massSO4[kg/m3]*1.e9ug/kg)^0.58
 676 !   or Na(m-3)= 1.e6* 340.*(1.e9ug/kg)^0.58 * (massSO4[kg/m3])^0.58
 677
 678            if(m .eq. idxsul) then
 679 !              naer2(i,k,m)= 5.64259e13_r8 * maerosol(1,m)**0.58
 680               naer2(i,k,m)=maerosol(1,m)*num_to_mass_aer(m)*2._r8
 681            else
 682               naer2(i,k,m)=maerosol(1,m)*num_to_mass_aer(m)
 683            endif
 684         enddo
 685 #endif
 686
 687 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 688 !ICE Nucleation
 689
 690         if (t(i,k).lt.tmelt - 5._r8) then
 691
 692
 693 ! get ice nucleation rate
 694
 695            call nucleati(wsubi(i,k),t(i,k),relhum(i,k),icldm(i,k),qc(i,k),nfice(i,k),rho(i,k),  &
 696 #ifdef MODAL_AERO
 697                          qaer(i,k,:)*rho(i,k),dgnum(i,k,:),1,naer_all,dum2,nihf2,niimm2,nidep2,nimey2)
 698 #else
 699                          naer2(i,k,:)/25._r8,1,naer_all,dum2,nihf2,niimm2,nidep2,nimey2)
 700 #endif
 701
 702            naai(i,k)=dum2
 703            naai_hom(i,k)=nihf2
 704            nihf(i,k)=nihf2
 705            niimm(i,k)=niimm2
 706            nidep(i,k)=nidep2
 707            nimey(i,k)=nimey2
 708         end if
 709
 710
 711 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 712 !droplet activation for bulk aerosol
 713
 714 #ifndef MODAL_AERO
 715         if (qc(i,k).ge.qsmall) then
 716
 717 ! get droplet activation rate
 718
 719         call activate(wsub(i,k),t(i,k),rho(i,k), &
 720                  naer2(i,k,:), naer_all,naer_all, maerosol,  &
 721                  dispersion_aer,hygro_aer, density_aer, dum2)
 722         dum = dum2
 723
 724         else
 725         dum = 0._r8
 726         end if
 727
 728  !note: deltat = deltat/2. to account for sub step in microphysics
 729         npccn(i,k) = (dum - nc(i,k)/lcldm(i,k))/(deltat/2._r8)*lcldm(i,k)
 730 #endif
 731
 732         end do ! i loop
 733         end do ! k loop
 734
 735 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 736 !droplet activation for modal aerosol
 737
 738 #ifdef MODAL_AERO
 739 !     partition cloud fraction into liquid water part
 740
 741       do k=1,pver
 742       do i=1,ncol
 743          qcld=qc(i,k)+qi(i,k)
 744          if(qcld.gt.qsmall)then
 745             lcldn(i,k)=cldn(i,k)*qc(i,k)/qcld
 746             lcldo(i,k)=cldo(i,k)*qc(i,k)/qcld
 747          else
 748             lcldn(i,k)=0._r8
 749             lcldo(i,k)=0._r8
 750          endif
 751       enddo
 752       enddo
 753
 754       call dropmixnuc(lchnk, ncol, ncloc, nctend_mixnuc, t, omega,  &
 755                     p, pint, pdel, rpdel, zm, kkvh, wsub, lcldn, lcldo,     &
 756                     qaer, cflx, qaertend, deltat                    &
 757 #ifdef WRF_PORT
 758                     , qqcw                                          &
 759 #endif
 760                                                                     )
 761
 762       npccn(:ncol,:)= nctend_mixnuc(:ncol,:)
 763 #endif
 764
 765
 766 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 767 ! Contact freezing  (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
 768 ! estimate rndst and nanco for 4 dust bins here to pass to MG microphysics
 769
 770
 771       do k=1,pver
 772       do i=1,ncol
 773
 774          if (t(i,k).lt.269.15_r8) then
 775
 776 #ifdef MODAL_AERO
 777 !BSINGH - Following #if added for MAM7 compliance
 778 #ifdef MODAL_AERO_7MODE
 779              nacon(i,k,3)=qaer(i,k,numptr_amode(modeptr_coardust))*rho(i,k)
 780              rndst(i,k,3)=0.5_r8*dgnumwet(i,k,modeptr_coardust)
 781 #elif ( defined MODAL_AERO_3MODE )
 782 ! For modal aerosols:
 783 !  use size '3' for dust coarse mode...
 784 !  scale by dust fraction in coarse mode
 785                dmc= qaer(i,k,lptr_dust_a_amode(modeptr_coarse))
 786                ssmc=qaer(i,k,lptr_nacl_a_amode(modeptr_coarse))
 787                if (dmc.gt.0.0_r8) then
 788                   nacon(i,k,3)=dmc/(ssmc+dmc) * qaer(i,k,numptr_amode(modeptr_coarse))*rho(i,k) !*tcnt
 789                else
 790                   nacon(i,k,3)=0._r8
 791                endif
 792
 793                !also redefine parameters based on size...
 794
 795                rndst(i,k,3)=0.5_r8*dgnumwet(i,k,modeptr_coarse)
 796
 797 !BSINGH - Following #endif added for MAM7 compliance
 798 # endif
 799            if (rndst(i,k,3).le.0._r8) then
 800               rndst(i,k,3)=rn_dst3
 801            endif
 802
 803 #else
 804
 805 !For Bulk Aerosols: set equal to aerosol number for dust for bins 2-4 (bin 1=0)
 806
 807            if (idxdst2.gt.0) then
 808               nacon(i,k,2)=naer2(i,k,idxdst2) !*tcnt ! 1/m3
 809            endif
 810            if (idxdst3.gt.0) then
 811               nacon(i,k,3)=naer2(i,k,idxdst3) !*tcnt
 812            endif
 813            if (idxdst4.gt.0) then
 814               nacon(i,k,4)=naer2(i,k,idxdst4) ! *tcnt
 815            endif
 816 #endif
 817         endif
 818       enddo
 819       enddo
 820
 821 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 822 !bulk aerosol ccn concentration (modal does it in ndrop, from dropmixnuc)
 823
 824 !        ccn concentration as diagnostic
 825 #ifndef MODAL_AERO
 826          do k=1,pver
 827             ccn(:ncol,k,:) = 0.
 828             do m=1,naer_all
 829                if(m.eq.idxsul)then
 830                ! Lohmann treatment for sulfate has variable size distribution
 831                   do i=1,ncol
 832                      if (naer2(i,k,m).gt.0._r8) then
 833                         amcubesulfate(i)=amcubefactor(m)*aer_mmr(i,k,m)*rho(i,k)/(naer2(i,k,m))
 834                         smcritsulfate(i)=smcritfactor(m)/sqrt(amcubesulfate(i))
 835                      else
 836                         smcritsulfate(i)=1._r8
 837                      endif
 838                   end do
 839                end if
 840                do l=1,psat
 841                   if(m.eq.idxsul)then
 842                      do i=1,ncol
 843                         arg=argfactor(m)*log(smcritsulfate(i)/super(l))
 844                         if(arg<2)then
 845                            if(arg<-2)then
 846                               ccnfact(l,m)=1.0e-6_r8
 847                            else
 848                               ccnfact(l,m)=0.5e-6_r8*erfc(arg)
 849                            endif
 850                         else
 851                            ccnfact(l,m)=0.0_r8
 852                         endif
 853                         ccn(i,k,l)=ccn(i,k,l)+naer2(i,k,m)*ccnfact(l,m)
 854                      end do
 855                   else
 856                   ccn(:ncol,k,l)=ccn(:ncol,k,l)+naer2(:ncol,k,m)*ccnfact(l,m)
 857                endif
 858             enddo
 859          enddo
 860       enddo
 861
 862       do l=1,psat
 863          call outfld(ccn_name(l), ccn(1,1,l)    , pcols, lchnk   )
 864       enddo
 865 #endif
 866
 867       do l=1,naer_all
 868          call outfld(aername(l)//'_m3', naer2(1,1,l), pcols, lchnk)
 869       enddo
 870
 871 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 872 ! output activated liquid and ice (convert from #/kg -> #/m3)
 873   do i = 1,ncol
 874      do k=1,pver
 875         nihf(i,k)=nihf(i,k)*rho(i,k)
 876         niimm(i,k)=niimm(i,k)*rho(i,k)
 877         nidep(i,k)=nidep(i,k)*rho(i,k)
 878         nimey(i,k)=nimey(i,k)*rho(i,k)
 879      end do
 880   end do
 881   call outfld('NIHF',nihf,    pcols,lchnk)
 882   call outfld('NIIMM',niimm,    pcols,lchnk)
 883   call outfld('NIDEP',nidep,    pcols,lchnk)
 884   call outfld('NIMEY',nimey,    pcols,lchnk)
 885
 886       deallocate( &
 887          naermod,  &
 888          naer2,    &
 889          maerosol, &
 890          ntype     )
 891
 892 return
 893 end subroutine microp_aero_ts
 894
 895 !##############################################################################
 896
 897       subroutine activate(wbar, tair, rhoair,  &
 898                           na, pmode, nmode, ma, sigman, hygro, rhodry, nact)
 899
 900 !      calculates number fraction of aerosols activated as CCN
 901 !      assumes an internal mixture within each of up to pmode multiple aerosol modes
 902 !      a gaussiam spectrum of updrafts can be treated.
 903
 904 !      mks units
 905
 906 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
 907 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
 908
 909       use physconst, only: rair, epsilo, cpair, rh2o, latvap, gravit,   &
 910                                  rhoh2o, mwh2o, r_universal
 911       use wv_saturation, only: estblf, epsqs
 912
 913       implicit none
 914
 915
 916 !      input
 917
 918       integer pmode,ptype ! dimension of modes, types in modes
 919       real(r8) wbar          ! grid cell mean vertical velocity (m/s)
 920       real(r8) tair          ! air temperature (K)
 921       real(r8) rhoair        ! air density (kg/m3)
 922       real(r8) na(pmode)           ! aerosol number concentration (/m3)
 923       integer nmode      ! number of aerosol modes
 924       real(r8) ma(pmode)     ! aerosol mass concentration (kg/m3)
 925       real(r8) rhodry(pmode) ! density of aerosol material
 926
 927       real(r8) sigman(pmode)  ! geometric standard deviation of aerosol size distribution
 928       real(r8) hygro(pmode)  ! hygroscopicity of aerosol mode
 929
 930
 931 !      output
 932
 933       real(r8) nact      ! number fraction of aerosols activated
 934
 935 !      local
 936
 937 !      external erf,erfc
 938 !      real(r8) erf,erfc
 939 #if (defined AIX)
 940 #define ERF erf
 941 #define ERFC erfc
 942 #else
 943 #define ERF derf
 944 #define ERFC derfc
 945       real(r8) derf,derfc
 946 #endif
 947
 948       integer, parameter:: nx=200
 949       integer :: maxmodes
 950       real(r8) surften       ! surface tension of water w/respect to air (N/m)
 951       data surften/0.076/
 952       save surften
 953       real(r8) p0     ! reference pressure (Pa)
 954       data p0/1013.25e2/
 955       save p0
 956
 957       real(r8), allocatable :: volc(:) ! total aerosol volume  concentration (m3/m3)
 958       real(r8) tmass ! total aerosol mass concentration (g/cm3)
 959       real(r8) rm ! number mode radius of aerosol at max supersat (cm)
 960       real(r8) pres ! pressure (Pa)
 961       real(r8) path ! mean free path (m)
 962       real(r8) diff ! diffusivity (m2/s)
 963       real(r8) conduct ! thermal conductivity (Joule/m/sec/deg)
 964       real(r8) diff0,conduct0
 965       real(r8) qs ! water vapor saturation mixing ratio
 966       real(r8) dqsdt ! change in qs with temperature
 967       real(r8) dqsdp ! change in qs with pressure
 968       real(r8) gloc ! thermodynamic function (m2/s)
 969       real(r8) zeta
 970       real(r8), allocatable :: eta(:)
 971       real(r8), allocatable :: smc(:)
 972       real(r8) lnsmax ! ln(smax)
 973       real(r8) alpha
 974       real(r8) gammaloc
 975       real(r8) beta
 976       real(r8) sqrtg
 977       real(r8) alogam
 978       real(r8) rlo,rhi,xint1,xint2,xint3,xint4
 979       real(r8) w,wnuc,wb
 980
 981       real(r8) alw,sqrtalw
 982       real(r8) smax
 983       real(r8) x,arg
 984       real(r8) xmincoeff,xcut,volcut,surfcut
 985       real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
 986       real(r8) :: etafactor1,etafactor2max
 987       real(r8),allocatable :: etafactor2(:)
 988       real(r8) es
 989       integer m,n
 990
 991       real(r8),allocatable :: amcubeloc(:)
 992       real(r8),allocatable :: lnsmloc(:)
 993
 994       maxmodes = naer_all
 995       allocate( &
 996          volc(maxmodes),       &
 997          eta(maxmodes),        &
 998          smc(maxmodes),        &
 999          etafactor2(maxmodes), &
1000          amcubeloc(maxmodes),  &
1001          lnsmloc(maxmodes)     )
1002
1003       if(maxmodes<pmode)then
1004          write(iulog,*)'maxmodes,pmode in activate =',maxmodes,pmode
1005 #ifdef WRF_PORT
1006          call wrf_message(iulog)
1007 #endif
1008          call endrun('activate')
1009       endif
1010
1011       nact=0._r8
1012
1013       if(nmode.eq.1.and.na(1).lt.1.e-20)return
1014
1015       if(wbar.le.0.)return
1016
1017       pres=rair*rhoair*tair
1018       diff0=0.211e-4*(p0/pres)*(tair/tt0)**1.94
1019       conduct0=(5.69+0.017*(tair-tt0))*4.186e2*1.e-5 ! convert to J/m/s/deg
1020       es = estblf(tair)
1021       qs = epsilo*es/(pres-(1.0_r8 - epsqs)*es)
1022       dqsdt=latvap/(rh2o*tair*tair)*qs
1023       alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1./(rair*tair))
1024       gammaloc=(1+latvap/cpair*dqsdt)/(rhoair*qs)
1025 !     growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
1026 !     should depend on mean radius of mode to account for gas kinetic effects
1027       gloc=1./(rhoh2o/(diff0*rhoair*qs)                                    &
1028           +latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1.))
1029       sqrtg=sqrt(gloc)
1030       beta=4.*pi*rhoh2o*gloc*gammaloc
1031       etafactor2max=1.e10/(alpha*wbar)**1.5 ! this should make eta big if na is very small.
1032
1033       do m=1,nmode
1034 !         internal mixture of aerosols
1035           volc(m)=ma(m)/(rhodry(m)) ! only if variable size dist
1036 ! pjj/cray - use TINY intrinsic
1037 !        if(volc(m).gt.1.e-39.and.na(m).gt.1.e-39)then
1038          if(volc(m).gt. TINY(volc) .and.na(m).gt. TINY(na))then
1039             etafactor2(m)=1./(na(m)*beta*sqrtg)  !fixed or variable size dist
1040 !            number mode radius (m)
1041             amcubeloc(m)=(3.*volc(m)/(4.*pi*exp45logsig(m)*na(m)))  ! only if variable size dist
1042             smc(m)=smcrit(m) ! only for prescribed size dist
1043
1044             if(hygro(m).gt.1.e-10)then   ! loop only if variable size dist
1045                smc(m)=2.*aten*sqrt(aten/(27.*hygro(m)*amcubeloc(m)))
1046             else
1047               smc(m)=100.
1048             endif
1049          else
1050             smc(m)=1.
1051             etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
1052          endif
1053          lnsmloc(m)=log(smc(m)) ! only if variable size dist
1054       enddo
1055
1056 !         single  updraft
1057          wnuc=wbar
1058
1059             w=wbar
1060             alw=alpha*wnuc
1061             sqrtalw=sqrt(alw)
1062             zeta=2.*sqrtalw*aten/(3.*sqrtg)
1063             etafactor1=2.*alw*sqrtalw
1064
1065             do m=1,nmode
1066                eta(m)=etafactor1*etafactor2(m)
1067             enddo
1068
1069             call maxsat(zeta,eta,nmode,smc,smax)
1070
1071             lnsmax=log(smax)
1072             xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
1073
1074             nact=0._r8
1075             do m=1,nmode
1076                x=2*(lnsmloc(m)-lnsmax)/(3*sq2*alogsig(m))
1077                nact=nact+0.5*(1.-erf(x))*na(m)
1078             enddo
1079             nact=nact/rhoair ! convert from #/m3 to #/kg
1080
1081       deallocate( &
1082          volc,       &
1083          eta,        &
1084          smc,        &
1085          etafactor2, &
1086          amcubeloc,  &
1087          lnsmloc     )
1088
1089       return
1090       end subroutine activate
1091
1092       subroutine maxsat(zeta,eta,nmode,smc,smax)
1093
1094 !      calculates maximum supersaturation for multiple
1095 !      competing aerosol modes.
1096
1097 !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1098 !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1099
1100       implicit none
1101
1102     ! integer pmode
1103     ! parameter (pmode=naer_all)
1104       integer nmode ! number of modes
1105       real(r8) :: smc(:) ! critical supersaturation for number mode radius
1106       real(r8) zeta
1107       real(r8) :: eta(:)
1108       real(r8) smax ! maximum supersaturation
1109       integer m  ! mode index
1110       real(r8) sum, g1, g2
1111
1112       do m=1,nmode
1113          if(zeta.gt.1.e5*eta(m).or.smc(m)*smc(m).gt.1.e5*eta(m))then
1114 !            weak forcing. essentially none activated
1115             smax=1.e-20
1116          else
1117 !            significant activation of this mode. calc activation all modes.
1118             go to 1
1119          endif
1120       enddo
1121
1122       return
1123
1124   1   continue
1125
1126       sum=0
1127       do m=1,nmode
1128          if(eta(m).gt.1.e-20)then
1129             g1=sqrt(zeta/eta(m))
1130             g1=g1*g1*g1
1131             g2=smc(m)/sqrt(eta(m)+3*zeta)
1132             g2=sqrt(g2)
1133             g2=g2*g2*g2
1134             sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
1135          else
1136             sum=1.e20
1137          endif
1138       enddo
1139
1140       smax=1./sqrt(sum)
1141
1142       return
1143
1144       end subroutine maxsat
1145
1146
1147 subroutine nucleati(wbar, tair, relhum, cldn, qc, nfice, rhoair, &
1148 #ifdef MODAL_AERO
1149        qaerpt, dgnum, ptype, naer_all, nuci, onihf, oniimm, onidep, onimey)
1150 #else
1151        na, ptype, naer_all, nuci, onihf, oniimm, onidep, onimey)
1152 #endif
1153
1154 !---------------------------------------------------------------
1155 ! Purpose:
1156 !  The parameterization of ice nucleation.
1157 !
1158 ! Method: The current method is based on Liu & Penner (2005)
1159 !  It related the ice nucleation with the aerosol number, temperature and the
1160 !  updraft velocity. It includes homogeneous freezing of sulfate, immersion
1161 !  freezing of soot, and Meyers et al. (1992) deposition nucleation
1162 !
1163 ! Authors: Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010
1164 !----------------------------------------------------------------
1165 #ifdef MODAL_AERO
1166 #ifndef WRF_PORT
1167       use modal_aero_data, only:numptr_amode, modeptr_accum, modeptr_coarse, modeptr_aitken, &
1168                                  ntot_amode, sigmag_amode, &
1169                                  lptr_dust_a_amode,lptr_nacl_a_amode
1170       use constituents, only: pcnst
1171 #else
1172       use modal_aero_data, only:numptr_amode => numptr_amode_mp, modeptr_accum => modeptr_accum_mp, &
1173            modeptr_coarse    => modeptr_coarse_mp   , modeptr_aitken    => modeptr_aitken_mp, &
1174            ntot_amode        => ntot_amode_mp       , sigmag_amode      => sigmag_amode_mp,   &
1175            lptr_dust_a_amode => lptr_dust_a_amode_mp, lptr_nacl_a_amode => lptr_nacl_a_amode_mp, &
1176            modeptr_coardust => modeptr_coardust_mp  !BSINGH - Added modeptr_coarsedust for MAM7 compliance
1177       use module_cam_support, only: pcnst =>pcnst_mp
1178 #endif
1179 #endif
1180
1181 !-----------------------------------------------------
1182 ! Input Arguments
1183 !
1184   integer ptype, naer_all
1185   real(r8), intent(in) :: wbar                ! grid cell mean vertical velocity (m/s)
1186   real(r8), intent(in) :: tair                ! temperature (K)
1187   real(r8), intent(in) :: relhum              ! relative humidity with respective to liquid
1188   real(r8), intent(in) :: cldn                ! new value of cloud fraction    (fraction)
1189   real(r8), intent(in) :: qc                  ! liquid water mixing ratio (kg/kg)
1190   real(r8), intent(in) :: nfice               ! ice mass fraction
1191   real(r8), intent(in) :: rhoair              ! air density (kg/m3)
1192 #ifdef MODAL_AERO
1193   real(r8), intent(in) :: qaerpt(pcnst) ! aerosol number and mass mixing ratios
1194   real(r8), intent(in) :: dgnum(ntot_amode)   ! aerosol mode dry diameter (m)
1195 #else
1196   real(r8), intent(in) :: na(naer_all)        ! aerosol number concentration (/m3)
1197 #endif
1198
1199
1200 !
1201 ! Output Arguments
1202 !
1203   real(r8), intent(out) :: nuci               ! ice number nucleated (#/kg)
1204   real(r8), intent(out) :: onihf              ! nucleated number from homogeneous freezing of so4
1205   real(r8), intent(out) :: oniimm             ! nucleated number from immersion freezing
1206   real(r8), intent(out) :: onidep             ! nucleated number from deposition nucleation
1207   real(r8), intent(out) :: onimey             ! nucleated number from deposition nucleation  (meyers: mixed phase)
1208
1209 !
1210 ! Local workspace
1211 !
1212   real(r8)  so4_num                                      ! so4 aerosol number (#/cm^3)
1213   real(r8)  soot_num                                     ! soot (hydrophilic) aerosol number (#/cm^3)
1214   real(r8)  dst1_num,dst2_num,dst3_num,dst4_num          ! dust aerosol number (#/cm^3)
1215   real(r8)  dst_num                                      ! total dust aerosol number (#/cm^3)
1216   real(r8)  nihf                                         ! nucleated number from homogeneous freezing of so4
1217   real(r8)  niimm                                        ! nucleated number from immersion freezing
1218   real(r8)  nidep                                        ! nucleated number from deposition nucleation
1219   real(r8)  nimey                                        ! nucleated number from deposition nucleation (meyers)
1220   real(r8)  n1,ni                                        ! nucleated number
1221   real(r8)  tc,A,B,C,regm,RHw                            ! work variable
1222   real(r8)  esl,esi,deles                                ! work variable
1223   real(r8)  dst_scale
1224   real(r8)  subgrid
1225   real(r8) dmc,ssmc         ! variables for modal scheme.
1226
1227     so4_num=0.0_r8
1228     soot_num=0.0_r8
1229     dst_num=0.0_r8
1230     dst1_num = 0.0_r8
1231     dst2_num = 0.0_r8
1232     dst3_num = 0.0_r8
1233     dst4_num = 0.0_r8
1234
1235 !For modal aerosols, assume for the upper troposphere:
1236 ! soot = accumulation mode
1237 ! sulfate = aiken mode
1238 ! dust = coarse mode
1239 ! since modal has internal mixtures.
1240
1241 #ifdef MODAL_AERO
1242       soot_num = qaerpt(numptr_amode(modeptr_accum))*1.0e-6_r8 !#/cm^3
1243 !BSINGH - Following #if added for MAM7 compliance
1244 #ifdef MODAL_AERO_7MODE
1245       dst_num=qaerpt(numptr_amode(modeptr_coardust))*1.0e-6_r8 !#/cm^3
1246 #elif ( defined MODAL_AERO_3MODE )
1247       dmc= qaerpt(lptr_dust_a_amode(modeptr_coarse))
1248       ssmc=qaerpt(lptr_nacl_a_amode(modeptr_coarse))
1249       if (dmc.gt.0._r8) then
1250          dst_num=dmc/(ssmc+dmc) * qaerpt(numptr_amode(modeptr_coarse))*1.0e-6_r8 !#/cm^3
1251       else
1252          dst_num=0.0_r8
1253       endif
1254 !BSINGH - Following #endif added for MAM7 compliance
1255 #endif
1256       if (dgnum(modeptr_aitken) .gt. 0._r8  ) then
1257          so4_num  = qaerpt(numptr_amode(modeptr_aitken))*1.0e-6_r8 & !#/cm^3, only allow so4 with D>0.1 um in ice nucleation
1258                * (0.5_r8 - 0.5_r8*erf(log(0.1e-6_r8/dgnum(modeptr_aitken))/  &
1259                  (2._r8**0.5_r8*log(sigmag_amode(modeptr_aitken)))))
1260       else
1261          so4_num = 0.0_r8
1262       endif
1263       so4_num = max(0.0_r8,so4_num)
1264 #else
1265     if(idxsul .gt. 0) then
1266        so4_num=na(idxsul)*1.0e-6_r8 ! #/cm^3
1267     end if
1268
1269 !continue above philosophy here....
1270
1271     if(idxbcphi .gt. 0) then
1272       soot_num=na(idxbcphi)*1.0e-6_r8 !#/cm^3
1273     end if
1274
1275     if(idxdst1 .gt. 0) then
1276        dst1_num=na(idxdst1)  *1.0e-6_r8 !#/cm^3
1277     end if
1278
1279     if(idxdst2 .gt. 0) then
1280        dst2_num=na(idxdst2)*1.0e-6_r8 !#/cm^3
1281     end if
1282
1283     if(idxdst3 .gt. 0) then
1284        dst3_num=na(idxdst3)*1.0e-6_r8 !#/cm^3
1285     end if
1286
1287     if(idxdst4 .gt. 0) then
1288        dst4_num=na(idxdst4)*1.0e-6_r8 !#/cm^3
1289     end if
1290
1291     dst_num =dst1_num+dst2_num+dst3_num+dst4_num
1292
1293 #endif
1294
1295 ! no soot nucleation
1296     soot_num=0.0_r8
1297
1298     ni=0._r8
1299     tc=tair-273.15_r8
1300
1301     ! initialize
1302     niimm=0._r8
1303     nidep=0._r8
1304     nihf=0._r8
1305
1306     if(so4_num.ge.1.0e-10_r8 .and. (soot_num+dst_num).ge.1.0e-10_r8 .and. cldn.gt.0._r8) then
1307
1308 !-----------------------------
1309 ! RHw parameterization for heterogeneous immersion nucleation
1310     A = 0.0073_r8
1311     B = 1.477_r8
1312     C = 131.74_r8
1313     RHw=(A*tc*tc+B*tc+C)*0.01_r8   ! RHi ~ 120-130%
1314
1315     subgrid = 1.2_r8
1316
1317     if((tc.le.-35.0_r8) .and. ((relhum*polysvp(tair,0)/polysvp(tair,1)*subgrid).ge.1.2_r8)) then ! use higher RHi threshold
1318
1319        A = -1.4938_r8 * log(soot_num+dst_num) + 12.884_r8
1320        B = -10.41_r8  * log(soot_num+dst_num) - 67.69_r8
1321        regm = A * log(wbar) + B
1322
1323        if(tc.gt.regm) then    ! heterogeneous nucleation only
1324          if(tc.lt.-40._r8 .and. wbar.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
1325            call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
1326            niimm=0._r8
1327            nidep=0._r8
1328            n1=nihf
1329          else
1330            call hetero(tc,wbar,soot_num+dst_num,niimm,nidep)
1331            nihf=0._r8
1332            n1=niimm+nidep
1333          endif
1334        elseif (tc.lt.regm-5._r8) then ! homogeneous nucleation only
1335          call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
1336          niimm=0._r8
1337          nidep=0._r8
1338          n1=nihf
1339        else        ! transition between homogeneous and heterogeneous: interpolate in-between
1340          if(tc.lt.-40._r8 .and. wbar.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
1341            call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
1342            niimm=0._r8
1343            nidep=0._r8
1344            n1=nihf
1345          else
1346
1347            call hf(regm-5._r8,wbar,relhum,subgrid,so4_num,nihf)
1348            call hetero(regm,wbar,soot_num+dst_num,niimm,nidep)
1349
1350            if(nihf.le.(niimm+nidep)) then
1351              n1=nihf
1352            else
1353              n1=(niimm+nidep)*((niimm+nidep)/nihf)**((tc-regm)/5._r8)
1354            endif
1355          endif
1356        endif
1357
1358        ni=n1
1359
1360     endif
1361     endif
1362 1100  continue
1363
1364 ! deposition/condensation nucleation in mixed clouds (-40<T<0C) (Meyers, 1992)
1365     if(tc.lt.0._r8 .and. tc.gt.-37._r8 .and. qc.gt.1.e-12_r8) then
1366       esl = polysvp(tair,0)     ! over water in mixed clouds
1367       esi = polysvp(tair,1)     ! over ice
1368       deles = (esl - esi)
1369       nimey=1.e-3_r8*exp(12.96_r8*deles/esi - 0.639_r8)
1370     else
1371       nimey=0._r8
1372     endif
1373
1374     nuci=ni+nimey
1375     if(nuci.gt.9999._r8.or.nuci.lt.0._r8) then
1376        write(iulog, *) 'Warning: incorrect ice nucleation number (nuci reset =0)'
1377 #ifdef WRF_PORT
1378        call wrf_message(iulog)
1379 #endif
1380        write(iulog, *) ni, tair, relhum, wbar, nihf, niimm, nidep,deles,esi,dst2_num,dst3_num,dst4_num
1381 #ifdef WRF_PORT
1382        call wrf_message(iulog)
1383 #endif
1384        nuci=0._r8
1385     endif
1386
1387     nuci=nuci*1.e+6_r8/rhoair    ! change unit from #/cm3 to #/kg
1388     onimey=nimey*1.e+6_r8/rhoair
1389     onidep=nidep*1.e+6_r8/rhoair
1390     oniimm=niimm*1.e+6_r8/rhoair
1391     onihf=nihf*1.e+6_r8/rhoair
1392
1393   return
1394   end subroutine nucleati
1395
1396   subroutine hetero(T,ww,Ns,Nis,Nid)
1397
1398     real(r8), intent(in)  :: T, ww, Ns
1399     real(r8), intent(out) :: Nis, Nid
1400
1401     real(r8) A11,A12,A21,A22,B11,B12,B21,B22
1402     real(r8) A,B,C
1403
1404 !---------------------------------------------------------------------
1405 ! parameters
1406
1407       A11 = 0.0263_r8
1408       A12 = -0.0185_r8
1409       A21 = 2.758_r8
1410       A22 = 1.3221_r8
1411       B11 = -0.008_r8
1412       B12 = -0.0468_r8
1413       B21 = -0.2667_r8
1414       B22 = -1.4588_r8
1415
1416 !     ice from immersion nucleation (cm^-3)
1417
1418       B = (A11+B11*log(Ns)) * log(ww) + (A12+B12*log(Ns))
1419       C =  A21+B21*log(Ns)
1420
1421       Nis = exp(A22) * Ns**B22 * exp(B*T) * ww**C
1422       Nis = min(Nis,Ns)
1423
1424       Nid = 0.0_r8    ! don't include deposition nucleation for cirrus clouds when T<-37C
1425
1426       return
1427   end subroutine hetero
1428
1429
1430   subroutine hf(T,ww,RH,subgrid,Na,Ni)
1431
1432       real(r8), intent(in)  :: T, ww, RH, subgrid, Na
1433       real(r8), intent(out) :: Ni
1434
1435       real(r8)    A1_fast,A21_fast,A22_fast,B1_fast,B21_fast,B22_fast
1436       real(r8)    A2_fast,B2_fast
1437       real(r8)    C1_fast,C2_fast,k1_fast,k2_fast
1438       real(r8)    A1_slow,A2_slow,B1_slow,B2_slow,B3_slow
1439       real(r8)    C1_slow,C2_slow,k1_slow,k2_slow
1440       real(r8)    regm
1441       real(r8)    A,B,C
1442       real(r8)    RHw
1443
1444 !---------------------------------------------------------------------
1445 ! parameters
1446
1447       A1_fast  =0.0231_r8
1448       A21_fast =-1.6387_r8  !(T>-64 deg)
1449       A22_fast =-6.045_r8   !(T<=-64 deg)
1450       B1_fast  =-0.008_r8
1451       B21_fast =-0.042_r8   !(T>-64 deg)
1452       B22_fast =-0.112_r8   !(T<=-64 deg)
1453       C1_fast  =0.0739_r8
1454       C2_fast  =1.2372_r8
1455
1456       A1_slow  =-0.3949_r8
1457       A2_slow  =1.282_r8
1458       B1_slow  =-0.0156_r8
1459       B2_slow  =0.0111_r8
1460       B3_slow  =0.0217_r8
1461       C1_slow  =0.120_r8
1462       C2_slow  =2.312_r8
1463
1464       Ni = 0.0_r8
1465
1466 !----------------------------
1467 !RHw parameters
1468       A = 6.0e-4_r8*log(ww)+6.6e-3_r8
1469       B = 6.0e-2_r8*log(ww)+1.052_r8
1470       C = 1.68_r8  *log(ww)+129.35_r8
1471       RHw=(A*T*T+B*T+C)*0.01_r8
1472
1473       if((T.le.-37.0_r8) .and. ((RH*subgrid).ge.RHw)) then
1474
1475         regm = 6.07_r8*log(ww)-55.0_r8
1476
1477         if(T.ge.regm) then    ! fast-growth regime
1478
1479           if(T.gt.-64.0_r8) then
1480             A2_fast=A21_fast
1481             B2_fast=B21_fast
1482           else
1483             A2_fast=A22_fast
1484             B2_fast=B22_fast
1485           endif
1486
1487           k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww))
1488           k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww)
1489
1490           Ni = k1_fast*Na**(k2_fast)
1491           Ni = min(Ni,Na)
1492
1493         else       ! slow-growth regime
1494
1495           k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww))
1496           k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww)
1497
1498           Ni = k1_slow*Na**(k2_slow)
1499           Ni = min(Ni,Na)
1500
1501         endif
1502
1503       end if
1504
1505       return
1506   end subroutine hf
1507
1508 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1509 ! error function in single precision
1510 !
1511 !    Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
1512 !    You may use, copy, modify this code for any purpose and
1513 !    without fee. You may distribute this ORIGINAL package.
1514
1515       real(r8) function derf(x)
1516       implicit real (a - h, o - z)
1517 ! pjj/cray function type
1518 !     real(r8) a,b,x
1519       real(r8) a,b,x
1520       dimension a(0 : 64), b(0 : 64)
1521       integer i,k
1522       data (a(i), i = 0, 12) / &
1523          0.00000000005958930743d0, -0.00000000113739022964d0, &
1524          0.00000001466005199839d0, -0.00000016350354461960d0, &
1525          0.00000164610044809620d0, -0.00001492559551950604d0, &
1526          0.00012055331122299265d0, -0.00085483269811296660d0, &
1527          0.00522397762482322257d0, -0.02686617064507733420d0, &
1528          0.11283791670954881569d0, -0.37612638903183748117d0, &
1529          1.12837916709551257377d0 /
1530       data (a(i), i = 13, 25) / &
1531          0.00000000002372510631d0, -0.00000000045493253732d0, &
1532          0.00000000590362766598d0, -0.00000006642090827576d0, &
1533          0.00000067595634268133d0, -0.00000621188515924000d0, &
1534          0.00005103883009709690d0, -0.00037015410692956173d0, &
1535          0.00233307631218880978d0, -0.01254988477182192210d0, &
1536          0.05657061146827041994d0, -0.21379664776456006580d0, &
1537          0.84270079294971486929d0 /
1538       data (a(i), i = 26, 38) / &
1539          0.00000000000949905026d0, -0.00000000018310229805d0, &
1540          0.00000000239463074000d0, -0.00000002721444369609d0, &
1541          0.00000028045522331686d0, -0.00000261830022482897d0, &
1542          0.00002195455056768781d0, -0.00016358986921372656d0, &
1543          0.00107052153564110318d0, -0.00608284718113590151d0, &
1544          0.02986978465246258244d0, -0.13055593046562267625d0, &
1545          0.67493323603965504676d0 /
1546       data (a(i), i = 39, 51) / &
1547          0.00000000000382722073d0, -0.00000000007421598602d0, &
1548          0.00000000097930574080d0, -0.00000001126008898854d0, &
1549          0.00000011775134830784d0, -0.00000111992758382650d0, &
1550          0.00000962023443095201d0, -0.00007404402135070773d0, &
1551          0.00050689993654144881d0, -0.00307553051439272889d0, &
1552          0.01668977892553165586d0, -0.08548534594781312114d0, &
1553          0.56909076642393639985d0 /
1554       data (a(i), i = 52, 64) / &
1555          0.00000000000155296588d0, -0.00000000003032205868d0, &
1556          0.00000000040424830707d0, -0.00000000471135111493d0, &
1557          0.00000005011915876293d0, -0.00000048722516178974d0, &
1558          0.00000430683284629395d0, -0.00003445026145385764d0, &
1559          0.00024879276133931664d0, -0.00162940941748079288d0, &
1560          0.00988786373932350462d0, -0.05962426839442303805d0, &
1561          0.49766113250947636708d0 /
1562       data (b(i), i = 0, 12) / &
1563          -0.00000000029734388465d0, 0.00000000269776334046d0, &
1564          -0.00000000640788827665d0, -0.00000001667820132100d0, &
1565          -0.00000021854388148686d0, 0.00000266246030457984d0, &
1566          0.00001612722157047886d0, -0.00025616361025506629d0, &
1567          0.00015380842432375365d0, 0.00815533022524927908d0, &
1568          -0.01402283663896319337d0, -0.19746892495383021487d0,&
1569          0.71511720328842845913d0 /
1570       data (b(i), i = 13, 25) / &
1571          -0.00000000001951073787d0, -0.00000000032302692214d0, &
1572          0.00000000522461866919d0, 0.00000000342940918551d0, &
1573          -0.00000035772874310272d0, 0.00000019999935792654d0, &
1574          0.00002687044575042908d0, -0.00011843240273775776d0, &
1575          -0.00080991728956032271d0, 0.00661062970502241174d0, &
1576          0.00909530922354827295d0, -0.20160072778491013140d0, &
1577          0.51169696718727644908d0 /
1578       data (b(i), i = 26, 38) / &
1579          0.00000000003147682272d0, -0.00000000048465972408d0, &
1580          0.00000000063675740242d0, 0.00000003377623323271d0, &
1581          -0.00000015451139637086d0, -0.00000203340624738438d0,&
1582          0.00001947204525295057d0, 0.00002854147231653228d0, &
1583          -0.00101565063152200272d0, 0.00271187003520095655d0, &
1584          0.02328095035422810727d0, -0.16725021123116877197d0, &
1585          0.32490054966649436974d0 /
1586       data (b(i), i = 39, 51) / &
1587          0.00000000002319363370d0, -0.00000000006303206648d0, &
1588          -0.00000000264888267434d0, 0.00000002050708040581d0, &
1589          0.00000011371857327578d0, -0.00000211211337219663d0, &
1590          0.00000368797328322935d0, 0.00009823686253424796d0, &
1591          -0.00065860243990455368d0, -0.00075285814895230877d0,&
1592          0.02585434424202960464d0, -0.11637092784486193258d0, &
1593          0.18267336775296612024d0 /
1594       data (b(i), i = 52, 64) / &
1595          -0.00000000000367789363d0, 0.00000000020876046746d0, &
1596          -0.00000000193319027226d0, -0.00000000435953392472d0, &
1597          0.00000018006992266137d0, -0.00000078441223763969d0, &
1598          -0.00000675407647949153d0, 0.00008428418334440096d0, &
1599          -0.00017604388937031815d0, -0.00239729611435071610d0, &
1600          0.02064129023876022970d0, -0.06905562880005864105d0, &
1601          0.09084526782065478489d0 /
1602       w = abs(x)
1603       if (w .lt. 2.2d0) then
1604           t = w * w
1605           k = int(t)
1606           t = t - k
1607           k = k * 13
1608           y = ((((((((((((a(k) * t + a(k + 1)) * t + &
1609              a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + &
1610              a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + &
1611              a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + &
1612              a(k + 11)) * t + a(k + 12)) * w
1613       else if (w .lt. 6.9d0) then
1614           k = int(w)
1615           t = w - k
1616           k = 13 * (k - 2)
1617           y = (((((((((((b(k) * t + b(k + 1)) * t + &
1618              b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + &
1619              b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + &
1620              b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + &
1621              b(k + 11)) * t + b(k + 12)
1622           y = y * y
1623           y = y * y
1624           y = y * y
1625           y = 1 - y * y
1626       else
1627           y = 1
1628       end if
1629       if (x .lt. 0) y = -y
1630       derf = y
1631       end function derf
1632 !
1633
1634 end module microp_aero
1635