chem/module_cam_mam_coag.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! module_cam_mam_coag.F
   4 ! adapted from cam3 modal_aero_coag.F90 by r.c.easter, august 2010
   5 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   6 ! 2010-08-04 rce notes
   7 !     > change pcnstxx to a subr argument because gas_pcnst is
   8 !       not a constant/parameter in wrf-chem
   9 !
  10 !--------------------------------------------------------------
  11 #include "MODAL_AERO_CPP_DEFINES.h"
  12 ! modal_aero_coag.F90
  13
  14
  15 !----------------------------------------------------------------------
  16 !BOP
  17 !
  18 ! !MODULE: modal_aero_coag --- modal aerosol coagulation
  19 !
  20 ! !INTERFACE:
  21    module modal_aero_coag
  22
  23 ! !USES:
  24    use shr_kind_mod,    only:  r8 => shr_kind_r8
  25    use shr_kind_mod,    only:  r4 => shr_kind_r4
  26 #ifndef WRF_PORT
  27    use chem_mods,       only:  gas_pcnst
  28 #else
  29    use module_cam_support, only:  gas_pcnst => gas_pcnst_modal_aero
  30 #endif
  31    use modal_aero_data, only:  maxd_aspectype
  32
  33   implicit none
  34   private
  35   save
  36
  37 ! !PUBLIC MEMBER FUNCTIONS:
  38   public modal_aero_coag_sub, modal_aero_coag_init
  39
  40 ! !PUBLIC DATA MEMBERS:
  41 #ifndef WRF_PORT
  42   integer, parameter :: pcnstxx = gas_pcnst
  43 #endif
  44
  45 #if ( defined MODAL_AERO_7MODE )
  46   integer, parameter, public :: pair_option_acoag = 3
  47 #elif ( defined MODAL_AERO_3MODE )
  48   integer, parameter, public :: pair_option_acoag = 1
  49 #endif
  50 ! specifies pairs of modes for which coagulation is calculated
  51 !   1 -- [aitken-->accum]
  52 !   2 -- [aitken-->accum], and [pcarbon-->accum]
  53 !   3 -- [aitken-->accum], [pcarbon-->accum],
  54 !        and [aitken-->pcarbon--(aging)-->accum]
  55 ! other -- do no coag
  56
  57   integer, parameter, public :: maxpair_acoag = 10
  58   integer, parameter, public :: maxspec_acoag = maxd_aspectype
  59
  60   integer, public :: npair_acoag
  61   integer, public :: modefrm_acoag(maxpair_acoag)
  62   integer, public :: modetoo_acoag(maxpair_acoag)
  63   integer, public :: modetooeff_acoag(maxpair_acoag)
  64   integer, public :: nspecfrm_acoag(maxpair_acoag)
  65   integer, public :: lspecfrm_acoag(maxspec_acoag,maxpair_acoag)
  66   integer, public :: lspectoo_acoag(maxspec_acoag,maxpair_acoag)
  67
  68 ! !DESCRIPTION: This module implements ...
  69 !
  70 ! !REVISION HISTORY:
  71 !
  72 !   RCE 07.04.13:  Adapted from MIRAGE2 code
  73 !
  74 !EOP
  75 !----------------------------------------------------------------------
  76 !BOC
  77
  78 ! list private module data here
  79
  80 !EOC
  81 !----------------------------------------------------------------------
  82
  83
  84   contains
  85
  86 !----------------------------------------------------------------------
  87 !BOP
  88 ! !ROUTINE:  modal_aero_coag_sub --- ...
  89 !
  90 ! !INTERFACE:
  91    subroutine modal_aero_coag_sub(                               &
  92                         lchnk,    ncol,     nstep,               &
  93                         loffset,  deltat_main,                   &
  94                         latndx,   lonndx,                        &
  95                         t,        pmid,     pdel,                &
  96                         q,                                       &
  97                         dgncur_a,           dgncur_awet,         &
  98                         wetdens_a                                &
  99 #ifdef WRF_PORT
 100                         , pcnstxx                                &
 101 #endif
 102                                                                  )
 103
 104
 105 !----------------------------------------------------------------------
 106 !   Authors: R. Easter
 107 !----------------------------------------------------------------------
 108
 109 ! !USES:
 110 #ifndef WRF_PORT
 111    use mo_constants,     only: pi
 112 #endif
 113    use modal_aero_data
 114    use modal_aero_gasaerexch, only:  n_so4_monolayers_pcage, &
 115                                      soa_equivso4_factor
 116 #ifndef WRF_PORT
 117    use abortutils,       only: endrun
 118    use cam_history,      only: outfld, fieldname_len
 119    use chem_mods,        only: adv_mass
 120    use constituents,     only: pcnst, cnst_name
 121 #else
 122    use module_cam_support, only: endrun, outfld, fieldname_len, pcnst =>pcnst_runtime
 123    use constituents,     only: cnst_name
 124    use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
 125 #endif
 126    use physconst,        only: gravit, mwdry, r_universal
 127 #ifndef WRF_PORT
 128    use ppgrid,           only: pcols, pver
 129    use spmd_utils,       only: iam, masterproc
 130 #else
 131    use module_cam_support, only: pcols, pver, iam, masterproc
 132 #endif
 133
 134    implicit none
 135
 136 ! !PARAMETERS:
 137 #ifdef WRF_PORT
 138    integer,  intent(in)    :: pcnstxx              ! number of species in "incoming q array"
 139 #endif
 140    integer, intent(in)  :: lchnk            ! chunk identifier
 141    integer, intent(in)  :: ncol             ! number of columns in chunk
 142    integer, intent(in)  :: nstep            ! model step
 143    integer, intent(in)  :: loffset          ! offset applied to modal aero "pointers"
 144    integer, intent(in)  :: latndx(pcols), lonndx(pcols)
 145
 146    real(r8), intent(in) :: deltat_main      ! model timestep (s)
 147
 148    real(r8), intent(in) :: t(pcols,pver)    ! temperature (K)
 149    real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
 150    real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
 151
 152    real(r8), intent(inout) :: q(ncol,pver,pcnstxx)
 153                                             ! tracer mixing ratio (TMR) array
 154                                             ! *** MUST BE mol/mol-air or #/mol-air
 155                                             ! *** NOTE ncol & pcnstxx dimensions
 156    real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode)
 157                                  ! dry geo. mean dia. (m) of number distrib.
 158    real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
 159                                  ! wet geo. mean dia. (m) of number distrib.
 160    real(r8), intent(in) :: wetdens_a(pcols,pver,ntot_amode)
 161                                  ! density of wet aerosol (kg/m3)
 162
 163 ! !DESCRIPTION:
 164 !   computes changes due to coagulation involving
 165 !       aitken mode (modeptr_aitken) with accumulation mode (modeptr_accum)
 166 !   this version will
 167 !       compute changes to mass and number, but not to surface area
 168 !       calculates coagulation rate coefficients using either
 169 !           new CMAQ V4.6 fast method
 170 !           older cmaq slow method (direct gauss-hermite quadrature)
 171 !
 172 ! !REVISION HISTORY:
 173 !   RCE 07.04.15:  Adapted from MIRAGE2 code and CMAQ V4.6 code
 174 !
 175 !EOP
 176 !----------------------------------------------------------------------
 177 !BOC
 178
 179 ! local variables
 180         integer :: i, iok, ipair, ip_aitacc, ip_aitpca, ip_pcaacc, iq
 181         integer :: idomode(ntot_amode), iselfcoagdone(ntot_amode)
 182         integer :: jfreqcoag
 183         integer :: k
 184         integer :: l, l1, l2, la, lmz, lsfrm, lstoo, lunout
 185         integer :: modefrm, modetoo, mait, macc, mpca
 186         integer ::  n, nfreqcoag
 187
 188
 189         integer, save :: nerr = 0       ! number of errors for entire run
 190         integer, save :: nerrmax = 9999 ! maximum number of errors before abort
 191         integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1
 192
 193         logical, parameter :: fastcoag_flag = .true. ! selects coag rate-coef method
 194
 195         real(r8) :: aircon
 196         real(r8) :: deltat, deltatinv_main
 197         real(r8) :: dr_so4_monolayers_pcage
 198         real(r8) :: dryvol_a(pcols,pver,ntot_amode)
 199         real(r8) :: dumexp, dumloss, dumprod
 200         real(r8) :: dumsfc_frm_old, dumsfc_frm_new
 201         real(r8) :: dum_m2v
 202         real(r8) :: fac_m2v_aitage(maxd_aspectype), fac_m2v_pcarbon(maxd_aspectype)
 203         real(r8) :: fac_volsfc_pcarbon
 204         real(r8) :: lnsg_frm, lnsg_too
 205         real(r8) :: sg_frm, sg_too
 206         real(r8) :: tmpa, tmpb, tmpc, tmpf, tmpg, tmph, tmpn
 207         real(r8) :: tmp1, tmp2
 208         real(r8) :: tmp_qold
 209         real(r8) :: v2ncur_a_tmp
 210         real(r8) :: vol_core, vol_shell
 211         real(r8) :: wetdens_frm, wetdens_too, wetdgnum_frm, wetdgnum_too
 212         real(r8) :: xbetaij0, xbetaij2i, xbetaij2j, xbetaij3, &
 213                     xbetaii0, xbetaii2,  xbetajj0,  xbetajj2
 214         real(r8) :: xferamt, xferfracvol, xferfrac_pcage, xferfrac_max
 215         real(r8) :: xnumbconc(ntot_amode)
 216         real(r8) :: xnumbconcavg(ntot_amode), xnumbconcnew(ntot_amode)
 217         real(r8) :: ybetaij0(maxpair_acoag), ybetaij3(maxpair_acoag)
 218         real(r8) :: ybetaii0(maxpair_acoag), ybetajj0(maxpair_acoag)
 219
 220         real(r8) :: dqdt(ncol,pver,pcnstxx)  ! TMR "dq/dt" array - NOTE dims
 221         logical  :: dotend(pcnst)            ! identifies the species that
 222                                              ! tendencies are computed for
 223         real(r8) :: qsrflx(pcols)
 224
 225         character(len=fieldname_len)   :: tmpname
 226         character(len=fieldname_len+3) :: fieldname
 227
 228 ! begin
 229 !   check if any coagulation pairs exist
 230         if (npair_acoag <= 0) return
 231
 232 !--------------------------------------------------------------------------------
 233    if (ldiag1 > 0) then
 234    if (nstep <= 3) then
 235    do i = 1, ncol
 236    if (lonndx(i) /= 37) cycle
 237    if (latndx(i) /= 23) cycle
 238    if (nstep > 3)       cycle
 239    write( *, '(/a,i7,i5,2(2x,2i5))' )   &
 240          '*** modal_aero_coag_sub -- nstep, iam, lat, lon, pcols, ncol =',   &
 241          nstep, iam, latndx(i), lonndx(i), pcols, ncol
 242    end do
 243    end if
 244 !  if (ncol /= -999888777) return
 245    if (nstep > 3) call endrun( 'modal_aero_coag_sub -- nstep>3 testing halt' )
 246    end if   ! (ldiag1 > 0)
 247 !--------------------------------------------------------------------------------
 248
 249         dotend(:) = .false.
 250         dqdt(1:ncol,:,:) = 0.0
 251
 252         lunout = 6
 253
 254
 255 !
 256 !   determine if coagulation will be done on this time-step
 257 !   currently coagulation is done every 3 hours
 258 !
 259 !       deltat = 3600.0*3.0
 260         deltat = deltat_main
 261         nfreqcoag = max( 1, nint( deltat/deltat_main ) )
 262         jfreqcoag = nfreqcoag/2
 263         xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8)   ! 1-eps
 264
 265         if (nfreqcoag .gt. 1) then
 266             if ( mod(nstep,nfreqcoag) .ne. jfreqcoag ) return
 267         end if
 268
 269 !
 270 !   set idomode
 271 !
 272         idomode(:) = 0
 273         do ipair = 1, npair_acoag
 274             idomode(modefrm_acoag(ipair)) = 1
 275             idomode(modetoo_acoag(ipair)) = 1
 276         end do
 277
 278 !
 279 !   other init
 280 !
 281         macc = modeptr_accum
 282         mait = modeptr_aitken
 283         mpca = modeptr_pcarbon
 284
 285         fac_m2v_aitage(:) = 0.0
 286         fac_m2v_pcarbon(:) = 0.0
 287         if (pair_option_acoag == 3) then
 288 !   following ipair definitions MUST BE CONSISTENT with
 289 !   the coding in modal_aero_coag_init for pair_option_acoag == 3
 290             ip_aitacc = 1
 291             ip_pcaacc = 2
 292             ip_aitpca = 3
 293
 294         ! use 1 mol (bi-)sulfate = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
 295             dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10
 296
 297             ipair = ip_aitpca
 298             do iq = 1, nspecfrm_acoag(ipair)
 299                 lsfrm = lspecfrm_acoag(iq,ipair)
 300                 if (lsfrm == lptr_so4_a_amode(mait)) then
 301                     fac_m2v_aitage(iq) = specmw_so4_amode / specdens_so4_amode
 302                 else if (lsfrm == lptr_nh4_a_amode(mait)) then
 303                     fac_m2v_aitage(iq) = specmw_nh4_amode / specdens_nh4_amode
 304                 else if (lsfrm == lptr_soa_a_amode(mait)) then
 305                     fac_m2v_aitage(iq) = soa_equivso4_factor*   &
 306                                         (specmw_soa_amode / specdens_soa_amode)
 307 !   for soa, the soa_equivso4_factor converts the soa volume into an
 308 !       so4(+nh4) volume that has same hygroscopicity contribution as soa
 309 !   this allows aging calculations to be done in terms of the amount
 310 !       of (equivalent) so4(+nh4) in the shell
 311 !   (see modal_aero_gasaerexch)
 312                 end if
 313             end do
 314
 315             do l = 1, nspec_amode(mpca)
 316                 l2 = lspectype_amode(l,mpca)
 317 !   fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
 318 !               [m3-AP/kmol-AP]    = [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
 319                 fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
 320             end do
 321
 322             fac_volsfc_pcarbon = exp( 2.5*(alnsg_amode(mpca)**2) )
 323         else
 324             ip_aitacc = -999888777
 325             ip_pcaacc = -999888777
 326             ip_aitpca = -999888777
 327         end if
 328
 329 !
 330 !   loop over levels and columns to calc the coagulation
 331 !
 332 !   integrate coagulation changes over deltat = nfreqcoag*deltat_main
 333 !   then compute tendencies as
 334 !      dqdt = (q(t+deltat) - q(t))/deltat_main
 335 !   because tendencies are applied (in physics_update) over deltat_main
 336 !
 337         deltat = nfreqcoag*deltat_main
 338         deltatinv_main = 1.0_r8/(deltat_main*(1.0_r8 + 1.0e-15_r8))
 339
 340 main_k: do k = 1, pver
 341 main_i: do i = 1, ncol
 342
 343 !   air molar density (kmol/m3)
 344         aircon = (pmid(i,k)/(r_universal*t(i,k)))
 345
 346 !   calculate number conc. (#/m3) for modes doing coagulation
 347         do n = 1, ntot_amode
 348             if (idomode(n) .gt. 0) then
 349                 xnumbconc(n) = q(i,k,numptr_amode(n)-loffset)*aircon
 350                 xnumbconc(n) = max( 0.0_r8, xnumbconc(n) )
 351             end if
 352             iselfcoagdone(n) = 0
 353         end do
 354
 355 !
 356 !   calculate coagulation rates for each pair
 357 !
 358 main_ipair1: do ipair = 1, npair_acoag
 359
 360         modefrm = modefrm_acoag(ipair)
 361         modetoo = modetoo_acoag(ipair)
 362
 363 !
 364 ! compute coagulation rates using cmaq "fast" method
 365 !    (based on E. Whitby's approximation approach)
 366 ! here subr. arguments are all in mks unit
 367 !
 368         call getcoags_wrapper_f(                                       &
 369           t(i,k), pmid(i,k),                                           &
 370           dgncur_awet(i,k,modefrm),     dgncur_awet(i,k,modetoo),      &
 371           sigmag_amode(modefrm),        sigmag_amode(modetoo),         &
 372           alnsg_amode(modefrm),         alnsg_amode(modetoo),          &
 373           wetdens_a(i,k,modefrm),       wetdens_a(i,k,modetoo),        &
 374           xbetaij0, xbetaij2i, xbetaij2j, xbetaij3,                    &
 375           xbetaii0, xbetaii2,  xbetajj0,  xbetajj2                     )
 376
 377
 378 !   test diagnostics begin --------------------------------------------
 379         if (ldiag2 > 0) then
 380         if (nstep <= 3) then
 381         if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
 382         if ((mod(k-1,5) == 0) .or. (k>=23)) then
 383
 384         wetdgnum_frm = dgncur_awet(i,k,modefrm)
 385         wetdgnum_too = dgncur_awet(i,k,modetoo)
 386         wetdens_frm  = wetdens_a(i,k,modefrm)
 387         wetdens_too  = wetdens_a(i,k,modetoo)
 388         sg_frm   = sigmag_amode(modefrm)
 389         sg_too   = sigmag_amode(modetoo)
 390         lnsg_frm = alnsg_amode(modefrm)
 391         lnsg_too = alnsg_amode(modetoo)
 392
 393         call getcoags_wrapper_f(                                       &
 394           t(i,k), pmid(i,k),                                           &
 395           wetdgnum_frm,                   wetdgnum_too,                &
 396           sg_frm,                         sg_too,                      &
 397           lnsg_frm,                       lnsg_too,                    &
 398           wetdens_frm,                    wetdens_too,                 &
 399           xbetaij0, xbetaij2i, xbetaij2j, xbetaij3,                    &
 400           xbetaii0, xbetaii2,  xbetajj0,  xbetajj2                     )
 401
 402
 403             write(lunout,9801)
 404             write(lunout,9810) 'nstep,lat,lon,k,ipair   ',   &
 405                 nstep, latndx(i), lonndx(i), k, ipair
 406             write(lunout,9820) 'tk, pmb, aircon, pdel   ',   &
 407                 t(i,k), pmid(i,k)*1.0e-2, aircon, pdel(i,k)*1.0e-2
 408             write(lunout,9820) 'wetdens-cgs, sg      f/t',   &
 409                 wetdens_frm*1.0e-3, wetdens_too*1.0e-3,   &
 410                 sg_frm, sg_too
 411             write(lunout,9820) 'dgnwet-um, dgndry-um f/t',   &
 412                 1.0e6*wetdgnum_frm, 1.0e6*wetdgnum_too,   &
 413                 1.0e6*dgncur_a(i,k,modefrm), 1.0e6*dgncur_a(i,k,modetoo)
 414             write(lunout,9820) 'xbeta ij0, ij3, ii0, jj0',   &
 415                 xbetaij0, xbetaij3, xbetaii0, xbetajj0
 416             write(lunout,9820) 'xbeta ij2i & j, ii2, jj2',   &
 417                 xbetaij2i, xbetaij2j, xbetaii2, xbetajj2
 418             write(lunout,9820) 'numbii, numbjj, deltat  ',   &
 419                 xnumbconc(modefrm), xnumbconc(modetoo), deltat
 420             write(lunout,9820) 'loss ij3, ii0, jj0      ',   &
 421                 (xbetaij3*xnumbconc(modetoo)*deltat),   &
 422                 (xbetaij0*xnumbconc(modetoo)*deltat+    &
 423                  xbetaii0*xnumbconc(modefrm)*deltat),   &
 424                 (xbetajj0*xnumbconc(modetoo)*deltat)
 425  9801   format( / 72x, 'ACOAG' )
 426  9810   format( 'ACOAG ', a, 2i8, 3i7, 3(1pe15.6) )
 427  9820   format( 'ACOAG ', a, 4(1pe15.6) )
 428  9830   format( 'ACOAG ', a, i1, a, 4(1pe15.6) )
 429         end if
 430         end if
 431         end if
 432         end if   ! (ldiag2 > 0)
 433 !   test diagnostics end ----------------------------------------------
 434
 435         ybetaij0(ipair) = xbetaij0
 436         ybetaij3(ipair) = xbetaij3
 437         ybetaii0(ipair) = xbetaii0
 438         ybetajj0(ipair) = xbetajj0
 439
 440         end do main_ipair1
 441
 442
 443
 444         if ( (pair_option_acoag == 1) .or.   &
 445              (pair_option_acoag == 2) ) then
 446 !
 447 !   calculate number and mass changes for pair_option_acoag == 1,2
 448 !
 449 main_ipair2: do ipair = 1, npair_acoag
 450
 451         modefrm = modefrm_acoag(ipair)
 452         modetoo = modetoo_acoag(ipair)
 453
 454 !   calculate number changes
 455 !   apply self-coagulation losses only once to a mode (when iselfcoagdone=0)
 456 !       first calc change to "too" mode
 457 !       next  calc change to "frm" mode, using average number conc of "too"
 458         if ( (mprognum_amode(modetoo) > 0) .and.   &
 459              (iselfcoagdone(modetoo) <= 0) ) then
 460             iselfcoagdone(modetoo) = 1
 461             tmpn = xnumbconc(modetoo)
 462             xnumbconcnew(modetoo) = tmpn/(1.0_r8 + deltat*ybetajj0(ipair)*tmpn)
 463             xnumbconcavg(modetoo) = 0.5_r8*(xnumbconcnew(modetoo) + tmpn)
 464             lstoo = numptr_amode(modetoo) - loffset
 465             q(i,k,lstoo) = xnumbconcnew(modetoo)/aircon
 466             dqdt(i,k,lstoo) = (xnumbconcnew(modetoo)-tmpn)*deltatinv_main/aircon
 467         end if
 468
 469         if ( (mprognum_amode(modefrm) > 0) .and.   &
 470              (iselfcoagdone(modefrm) <= 0) ) then
 471             iselfcoagdone(modefrm) = 1
 472             tmpn = xnumbconc(modefrm)
 473             tmpa = deltat*ybetaij0(ipair)*xnumbconcavg(modetoo)
 474             tmpb = deltat*ybetaii0(ipair)
 475             tmpc = tmpa + tmpb*tmpn
 476             if (abs(tmpc) < 0.01_r8) then
 477                 xnumbconcnew(modefrm) = tmpn*exp(-tmpc)
 478             else if (abs(tmpa) < 0.001_r8) then
 479                 xnumbconcnew(modefrm) =   &
 480                     exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
 481             else
 482                 tmpf = tmpb*tmpn/tmpc
 483                 tmpg = exp(-tmpa)
 484                 tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
 485                 xnumbconcnew(modefrm) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
 486             end if
 487             xnumbconcavg(modefrm) = 0.5_r8*(xnumbconcnew(modefrm) + tmpn)
 488             lsfrm = numptr_amode(modefrm) - loffset
 489             q(i,k,lsfrm) = xnumbconcnew(modefrm)/aircon
 490             dqdt(i,k,lsfrm) = (xnumbconcnew(modefrm)-tmpn)*deltatinv_main/aircon
 491         end if
 492
 493 !   calculate mass changes
 494 !     xbetaij3*xnumbconc(modetoo) = first order loss rate for modefrm volume
 495 !     xferfracvol = fraction of modefrm volume transferred to modetoo over deltat
 496         dumloss = ybetaij3(ipair)*xnumbconcavg(modetoo)
 497         xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
 498         xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
 499
 500         do iq = 1, nspecfrm_acoag(ipair)
 501             lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 502             lstoo = lspectoo_acoag(iq,ipair) - loffset
 503             if (lsfrm > 0) then
 504                 xferamt = q(i,k,lsfrm)*xferfracvol
 505                 dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
 506                 q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
 507                 if (lstoo > 0) then
 508                     dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
 509                     q(i,k,lstoo) = q(i,k,lstoo) + xferamt
 510                 end if
 511             end if
 512         end do
 513
 514         end do main_ipair2
 515
 516
 517         else if (pair_option_acoag == 3) then
 518 !
 519 !   calculate number and mass changes for pair_option_acoag == 3
 520 !
 521
 522 !   calculate number changes to accum mode
 523         if (mprognum_amode(macc) > 0) then
 524             tmpn = xnumbconc(macc)
 525             xnumbconcnew(macc) = tmpn/(1.0_r8 + deltat*ybetajj0(ip_aitacc)*tmpn)
 526             xnumbconcavg(macc) = 0.5_r8*(xnumbconcnew(macc) + tmpn)
 527             lstoo = numptr_amode(macc) - loffset
 528             q(i,k,lstoo) = xnumbconcnew(macc)/aircon
 529             dqdt(i,k,lstoo) = (xnumbconcnew(macc)-tmpn)*deltatinv_main/aircon
 530         end if
 531
 532 !   calculate number changes to primary carbon mode
 533         modefrm = modeptr_pcarbon
 534         if (mprognum_amode(mpca) > 0) then
 535             tmpn = xnumbconc(mpca)
 536             tmpa = deltat*ybetaij0(ip_pcaacc)*xnumbconcavg(macc)
 537             tmpb = deltat*ybetaii0(ip_pcaacc)
 538             tmpc = tmpa + tmpb*tmpn
 539             if (abs(tmpc) < 0.01_r8) then
 540                 xnumbconcnew(mpca) = tmpn*exp(-tmpc)
 541             else if (abs(tmpa) < 0.001_r8) then
 542                 xnumbconcnew(mpca) =   &
 543                     exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
 544             else
 545                 tmpf = tmpb*tmpn/tmpc
 546                 tmpg = exp(-tmpa)
 547                 tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
 548                 xnumbconcnew(mpca) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
 549             end if
 550             xnumbconcavg(mpca) = 0.5_r8*(xnumbconcnew(mpca) + tmpn)
 551             lsfrm = numptr_amode(mpca) - loffset
 552             q(i,k,lsfrm) = xnumbconcnew(mpca)/aircon
 553             dqdt(i,k,lsfrm) = (xnumbconcnew(mpca)-tmpn)*deltatinv_main/aircon
 554         end if
 555
 556 !   calculate number changes to aitken mode
 557         if (mprognum_amode(mait) > 0) then
 558             tmpn = xnumbconc(mait)
 559             tmpa = deltat*( ybetaij0(ip_aitacc)*xnumbconcavg(macc)   &
 560                           + ybetaij0(ip_aitpca)*xnumbconcavg(mpca) )
 561             tmpb = deltat*ybetaii0(ip_aitacc)
 562             tmpc = tmpa + tmpb*tmpn
 563             if (abs(tmpc) < 0.01_r8) then
 564                 xnumbconcnew(mait) = tmpn*exp(-tmpc)
 565             else if (abs(tmpa) < 0.001_r8) then
 566                 xnumbconcnew(mait) =   &
 567                     exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
 568             else
 569                 tmpf = tmpb*tmpn/tmpc
 570                 tmpg = exp(-tmpa)
 571                 tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
 572                 xnumbconcnew(mait) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
 573             end if
 574             xnumbconcavg(mait) = 0.5_r8*(xnumbconcnew(mait) + tmpn)
 575             lsfrm = numptr_amode(mait) - loffset
 576             q(i,k,lsfrm) = xnumbconcnew(mait)/aircon
 577             dqdt(i,k,lsfrm) = (xnumbconcnew(mait)-tmpn)*deltatinv_main/aircon
 578         end if
 579
 580
 581 !   calculate mass changes from aitken-->accum direct coagulation and
 582 !       aitken-->pcarbon-->accum coagulation/aging
 583 !   also calc volume of shell material (so4 & nh4 from aitken-->pcarbon)
 584         dumloss = ybetaij3(ip_aitacc)*xnumbconcavg(macc)   &
 585                 + ybetaij3(ip_aitpca)*xnumbconcavg(mpca)
 586         tmpa = ybetaij3(ip_aitpca)*xnumbconcavg(mpca)/max( dumloss, 1.0e-37_r8 )
 587         xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
 588         xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
 589         vol_shell = 0.0
 590
 591         ipair = ip_aitacc
 592         do iq = 1, nspecfrm_acoag(ipair)
 593             lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 594             lstoo = lspectoo_acoag(iq,ipair) - loffset
 595             if (lsfrm > 0) then
 596                 xferamt = q(i,k,lsfrm)*xferfracvol
 597                 dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
 598                 q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
 599                 if (lstoo > 0) then
 600                     dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
 601                     q(i,k,lstoo) = q(i,k,lstoo) + xferamt
 602                 end if
 603                 vol_shell = vol_shell + xferamt*tmpa*fac_m2v_aitage(iq)
 604             end if
 605         end do
 606
 607
 608 !   now calculate aging transfer fraction for pcarbon-->accum
 609 !   this duplicates the code in modal_aero_gasaerexch
 610         vol_core = 0.0
 611         do l = 1, nspec_amode(mpca)
 612             vol_core = vol_core + &
 613                 q(i,k,lmassptr_amode(l,mpca)-loffset)*fac_m2v_pcarbon(l)
 614         end do
 615         tmp1 = vol_shell*dgncur_a(i,k,mpca)*fac_volsfc_pcarbon
 616         tmp2 = 6.0_r8*dr_so4_monolayers_pcage*vol_core
 617         tmp2 = max( tmp2, 0.0_r8 )
 618         if (tmp1 >= tmp2) then
 619             xferfrac_pcage = xferfrac_max
 620         else
 621             xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
 622         end if
 623
 624
 625 !   calculate mass changes from pcarbon-->accum by direct coagulation
 626 !   and aging
 627         dumloss = ybetaij3(ip_pcaacc)*xnumbconcavg(macc)
 628         xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
 629         xferfracvol = xferfracvol + xferfrac_pcage
 630         xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
 631
 632         ipair = ip_pcaacc
 633         do iq = 1, nspecfrm_acoag(ipair)
 634             lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 635             lstoo = lspectoo_acoag(iq,ipair) - loffset
 636             if (lsfrm > 0) then
 637                 xferamt = q(i,k,lsfrm)*xferfracvol
 638                 dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
 639                 q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
 640                 if (lstoo > 0) then
 641                     dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
 642                     q(i,k,lstoo) = q(i,k,lstoo) + xferamt
 643                 end if
 644             end if
 645         end do
 646
 647         lsfrm = numptr_amode(mpca) - loffset
 648         lstoo = numptr_amode(macc) - loffset
 649         if (lsfrm > 0) then
 650             xferamt = q(i,k,lsfrm)*xferfrac_pcage
 651             dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
 652             q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
 653             if (lstoo > 0) then
 654                 dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
 655                 q(i,k,lstoo) = q(i,k,lstoo) + xferamt
 656             end if
 657         end if
 658
 659
 660
 661         else   ! (pair_option_acoag /= 1,2,3) then
 662
 663         write(lunout,*) '*** modal_aero_coag_sub error'
 664         write(lunout,*) '    cannot do _coag_sub error pair_option_acoag =', &
 665                 pair_option_acoag
 666         call endrun( 'modal_aero_coag_sub error' )
 667
 668
 669         end if   ! (pair_option_acoag == ...)
 670
 671
 672 !   test diagnostics begin --------------------------------------------
 673         if (ldiag3 > 0) then
 674         if (nstep <= 3) then
 675         if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
 676         if ((mod(k-1,5) == 0) .or. (k>=23)) then
 677            if (pair_option_acoag == 3) then
 678                 write(*,*)
 679                 write(lunout,9820) 'xnumbconcavg ait,acc,pca', &
 680                     xnumbconcavg(mait), xnumbconcavg(macc), xnumbconcavg(mpca)
 681                 write(lunout,9820) 'vshell, core            ', &
 682                     vol_shell, vol_core
 683                 write(lunout,9820) 'dr_mono, dgn            ', &
 684                     dr_so4_monolayers_pcage, dgncur_a(i,k,mpca)
 685                 write(lunout,9820) 'tmp1, tmp2              ', tmp1, tmp2
 686                 write(lunout,9820) 'xferfrac_age            ', xferfrac_pcage
 687            end if
 688
 689            do ipair = 1, npair_acoag
 690            modefrm = modefrm_acoag(ipair)
 691            modetoo = modetoo_acoag(ipair)
 692            if (npair_acoag > 1) then
 693                 write(lunout,*)
 694                 write(lunout,9810) 'ipair =   ', ipair
 695            end if
 696
 697            do iq = 1, nspecfrm_acoag(ipair)
 698            lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 699            lstoo = lspectoo_acoag(iq,ipair) - loffset
 700            if (lsfrm > 0) then
 701            tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
 702 !          write(lunout,9820) 'm1 frm dqdt/q0,dqdt,q0/1',   &
 703            write(lunout,9830) 'm', iq,   &
 704                                 ' frm dqdt/q0,dqdt,q0/1',   &
 705                 dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
 706            end if
 707            if (lstoo > 0) then
 708            tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
 709            write(lunout,9830) 'm', iq,   &
 710                                 ' too dqdt/q0,dqdt,q0/1',   &
 711                 dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
 712            end if
 713            end do   ! iq
 714
 715            lsfrm = numptr_amode(modefrm) - loffset
 716            lstoo = numptr_amode(modetoo) - loffset
 717            if (lsfrm > 0) then
 718            tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
 719            write(lunout,9820) 'n  frm dqdt/q0,dqdt,q0/1',   &
 720                 dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
 721            end if
 722            if (lstoo > 0) then
 723            tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
 724            write(lunout,9820) 'n  too dqdt/q0,dqdt,q0/1',   &
 725                 dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
 726            end if
 727
 728            end do   ! ipair
 729         end if
 730         end if
 731         end if
 732         end if   ! (ldiag3 > 0)
 733 !   test diagnostics end ----------------------------------------------
 734
 735
 736
 737         end do main_i
 738         end do main_k
 739
 740
 741 ! set dotend's
 742         do ipair = 1, npair_acoag
 743             modefrm = modefrm_acoag(ipair)
 744             modetoo = modetoo_acoag(ipair)
 745
 746             do iq = 1, nspecfrm_acoag(ipair)
 747                 lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 748                 lstoo = lspectoo_acoag(iq,ipair) - loffset
 749                 if (lsfrm > 0) dotend(lsfrm) = .true.
 750                 if (lstoo > 0) dotend(lstoo) = .true.
 751             end do
 752
 753             if (mprognum_amode(modefrm) > 0) then
 754                 lsfrm = numptr_amode(modefrm) - loffset
 755                 if (lsfrm > 0) dotend(lsfrm) = .true.
 756             end if
 757             if (mprognum_amode(modetoo) > 0) then
 758                 lstoo = numptr_amode(modetoo) - loffset
 759                 if (lstoo > 0) dotend(lstoo) = .true.
 760             end if
 761
 762         end do
 763
 764
 765 !   do history file column-tendency fields
 766         do l = loffset+1, pcnst
 767             lmz = l - loffset
 768             if ( .not. dotend(lmz) ) cycle
 769
 770             qsrflx(:) = 0.0_r8
 771             do k = 1, pver
 772             do i = 1, ncol
 773                 qsrflx(i) = qsrflx(i) + dqdt(i,k,lmz)*pdel(i,k)
 774             end do
 775             end do
 776             qsrflx(:) = qsrflx(:)*(adv_mass(lmz)/(gravit*mwdry))
 777             fieldname = trim(cnst_name(l)) // '_sfcoag1'
 778             call outfld( fieldname, qsrflx, pcols, lchnk )
 779 !           if (( masterproc ) .and. (nstep < 1)) &
 780 !               write(*,'(2(a,2x),1p,e11.3)') &
 781 !               'modal_aero_coag_sub outfld', fieldname, adv_mass(lmz)
 782         end do ! l = ...
 783
 784
 785         return
 786
 787
 788 !EOC
 789         end subroutine modal_aero_coag_sub
 790
 791
 792 !----------------------------------------------------------------------
 793 !----------------------------------------------------------------------
 794         subroutine modal_aero_coag_init
 795 !
 796 !   computes pointers for species transfer during coagulation
 797 !
 798         use modal_aero_data
 799         use modal_aero_gasaerexch, only:  &
 800                 modefrm_pcage, nspecfrm_pcage, lspecfrm_pcage, lspectoo_pcage
 801 #ifndef WRF_PORT
 802         use abortutils,      only: endrun
 803         use cam_history,     only: addfld, add_default, fieldname_len, phys_decomp
 804         use constituents,    only: pcnst, cnst_name
 805         use spmd_utils,      only: masterproc
 806         use phys_control,    only: phys_getopts
 807 #else
 808         use module_cam_support, only: endrun, addfld, add_default, fieldname_len, phys_decomp, &
 809              pcnst =>pcnst_runtime, masterproc
 810         use constituents,    only: cnst_name
 811 #endif
 812
 813         implicit none
 814
 815 !   local variables
 816         integer :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
 817         integer :: l, lsfrm, lstoo, lunout
 818         integer :: m, mfrm, mtoo, mtef
 819         integer :: nsamefrm, nsametoo, nspec
 820
 821         character(len=fieldname_len)   :: tmpname
 822         character(len=fieldname_len+3) :: fieldname
 823         character(128)                 :: long_name
 824         character(8)                   :: unit
 825
 826         logical :: dotend(pcnst)
 827         logical :: history_aerosol      ! Output the MAM aerosol tendencies
 828
 829         !-----------------------------------------------------------------------
 830 #ifndef WRF_PORT
 831         call phys_getopts( history_aerosol_out        = history_aerosol   )
 832 #else
 833         history_aerosol =.false.
 834 #endif
 835
 836         lunout = 6
 837 !
 838 !   define "from mode" and "to mode" for each coagulation pairing
 839 !       currently just a2-->a1 coagulation
 840 !
 841         if (pair_option_acoag == 1) then
 842             npair_acoag = 1
 843             modefrm_acoag(1) = modeptr_aitken
 844             modetoo_acoag(1) = modeptr_accum
 845             modetooeff_acoag(1) = modeptr_accum
 846         else if (pair_option_acoag == 2) then
 847             npair_acoag = 2
 848             modefrm_acoag(1) = modeptr_aitken
 849             modetoo_acoag(1) = modeptr_accum
 850             modetooeff_acoag(1) = modeptr_accum
 851             modefrm_acoag(2) = modeptr_pcarbon
 852             modetoo_acoag(2) = modeptr_accum
 853             modetooeff_acoag(2) = modeptr_accum
 854         else if (pair_option_acoag == 3) then
 855             npair_acoag = 3
 856             modefrm_acoag(1) = modeptr_aitken
 857             modetoo_acoag(1) = modeptr_accum
 858             modetooeff_acoag(1) = modeptr_accum
 859             modefrm_acoag(2) = modeptr_pcarbon
 860             modetoo_acoag(2) = modeptr_accum
 861             modetooeff_acoag(2) = modeptr_accum
 862             modefrm_acoag(3) = modeptr_aitken
 863             modetoo_acoag(3) = modeptr_pcarbon
 864             modetooeff_acoag(3) = modeptr_accum
 865             if (modefrm_pcage <= 0) then
 866                 write(*,*) '*** modal_aero_coag_init error'
 867                 write(*,*) '    pair_option_acoag, modefrm_pcage mismatch'
 868                 write(*,*) '    pair_option_acoag, modefrm_pcage =', &
 869                     pair_option_acoag, modefrm_pcage
 870                 call endrun( 'modal_aero_coag_init error' )
 871             end if
 872         else
 873             npair_acoag = 0
 874             return
 875         end if
 876
 877 !
 878 !   define species involved in each coagulation pairing
 879 !       (include aerosol water)
 880 !
 881 aa_ipair: do ipair = 1, npair_acoag
 882
 883         mfrm = modefrm_acoag(ipair)
 884         mtoo = modetoo_acoag(ipair)
 885         mtef = modetooeff_acoag(ipair)
 886         if ( (mfrm < 1) .or. (mfrm > ntot_amode) .or.   &
 887              (mtoo < 1) .or. (mtoo > ntot_amode) .or.   &
 888              (mtef < 1) .or. (mtef > ntot_amode) ) then
 889             write(*,*) '*** modal_aero_coag_init error'
 890             write(*,*) '    ipair, ntot_amode =', ipair, ntot_amode
 891             write(*,*) '    mfrm, mtoo, mtef  =', mfrm, mtoo, mtef
 892             call endrun( 'modal_aero_coag_init error' )
 893         end if
 894
 895
 896         mtoo = mtef   ! effective modetoo
 897         nspec = 0
 898 aa_iqfrm: do iqfrm = 1, nspec_amode(mfrm)
 899             lsfrm = lmassptr_amode(iqfrm,mfrm)
 900             if ((lsfrm .lt. 1) .or. (lsfrm .gt. pcnst)) cycle aa_iqfrm
 901
 902 ! find "too" species having same lspectype_amode as the "frm" species
 903 ! several species in a mode may have the same lspectype_amode, so also
 904 !    use the ordering as a criterion (e.g., 1st <--> 1st, 2nd <--> 2nd)
 905             iqfrm_aa = 1
 906             iqtoo_aa = 1
 907             if (iqfrm .gt. nspec_amode(mfrm)) then
 908                 iqfrm_aa = nspec_amode(mfrm) + 1
 909                 iqtoo_aa = nspec_amode(mtoo) + 1
 910             end if
 911             nsamefrm = 0
 912             do iq = iqfrm_aa, iqfrm
 913                 if ( lspectype_amode(iq   ,mfrm) .eq.   &
 914                      lspectype_amode(iqfrm,mfrm) ) then
 915                     nsamefrm = nsamefrm + 1
 916                 end if
 917             end do
 918             nsametoo = 0
 919             lstoo = 0
 920             do iqtoo = iqtoo_aa, nspec_amode(mtoo)
 921                 if ( lspectype_amode(iqtoo,mtoo) .eq.   &
 922                      lspectype_amode(iqfrm,mfrm) ) then
 923                     nsametoo = nsametoo + 1
 924                     if (nsametoo .eq. nsamefrm) then
 925                         lstoo = lmassptr_amode(iqtoo,mtoo)
 926                         exit
 927                     end if
 928                 end if
 929             end do
 930
 931             nspec = nspec + 1
 932             lspecfrm_acoag(nspec,ipair) = lsfrm
 933             lspectoo_acoag(nspec,ipair) = lstoo
 934         end do aa_iqfrm
 935
 936 !       lsfrm = lwaterptr_amode(mfrm)
 937 !       if ((lsfrm .ge. 1) .and. (lsfrm .le. pcnst)) then
 938 !           lstoo = lwaterptr_amode(mtoo)
 939 !           if ((lstoo .lt. 1) .or. (lstoo .gt. pcnst)) lstoo = 0
 940 !           nspec = nspec + 1
 941 !           lspecfrm_acoag(nspec,ipair) = lsfrm
 942 !           lspectoo_acoag(nspec,ipair) = lstoo
 943 !       end if
 944
 945         nspecfrm_acoag(ipair) = nspec
 946         end do aa_ipair
 947
 948 !
 949 !   output results
 950 !
 951         if ( masterproc ) then
 952
 953         write(lunout,9310)
 954
 955         do ipair = 1, npair_acoag
 956           mfrm = modefrm_acoag(ipair)
 957           mtoo = modetoo_acoag(ipair)
 958           mtef = modetooeff_acoag(ipair)
 959           write(lunout,9320) ipair, mfrm, mtoo, mtef
 960
 961           do iq = 1, nspecfrm_acoag(ipair)
 962             lsfrm = lspecfrm_acoag(iq,ipair)
 963             lstoo = lspectoo_acoag(iq,ipair)
 964             if (lstoo .gt. 0) then
 965                 write(lunout,9330) lsfrm, cnst_name(lsfrm),   &
 966                         lstoo, cnst_name(lstoo)
 967             else
 968                 write(lunout,9340) lsfrm, cnst_name(lsfrm)
 969             end if
 970           end do
 971
 972         end do ! ipair = ...
 973         write(lunout,*)
 974
 975         end if ! ( masterproc )
 976
 977 9310    format( / 'subr. modal_aero_coag_init' )
 978 9320    format( 'pair', i3, 5x, 'mode', i3, &
 979                 ' ---> mode', i3, '   eff', i3 )
 980 9330    format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
 981 9340    format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )
 982
 983
 984 !
 985 !   create history file column-tendency fields
 986 !
 987         dotend(:) = .false.
 988         do ipair = 1, npair_acoag
 989           do iq = 1, nspecfrm_acoag(ipair)
 990             l = lspecfrm_acoag(iq,ipair)
 991             if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
 992             l = lspectoo_acoag(iq,ipair)
 993             if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
 994           end do
 995
 996           m = modefrm_acoag(ipair)
 997           if ((m > 0) .and. (m <= ntot_amode)) then
 998             l = numptr_amode(m)
 999             if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
1000           end if
1001           m = modetoo_acoag(ipair)
1002           if ((m > 0) .and. (m <= ntot_amode)) then
1003             l = numptr_amode(m)
1004             if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
1005           end if
1006         end do ! ipair = ...
1007
1008         if (pair_option_acoag == 3) then
1009            do iq = 1, nspecfrm_pcage
1010               lsfrm = lspecfrm_pcage(iq)
1011               lstoo = lspectoo_pcage(iq)
1012               if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1013                  dotend(lsfrm) = .true.
1014                  if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1015                     dotend(lstoo) = .true.
1016                  end if
1017               end if
1018            end do
1019         end if
1020
1021         do l = 1, pcnst
1022             if ( .not. dotend(l) ) cycle
1023             tmpname = cnst_name(l)
1024             unit = 'kg/m2/s'
1025             do m = 1, ntot_amode
1026                 if (l == numptr_amode(m)) unit = '#/m2/s'
1027             end do
1028             fieldname = trim(tmpname) // '_sfcoag1'
1029             long_name = trim(tmpname) // ' modal_aero coagulation column tendency'
1030             call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1031             if ( history_aerosol ) then
1032                call add_default( fieldname, 1, ' ' )
1033             endif
1034             if ( masterproc ) write(*,'(3(a,2x))') &
1035                 'modal_aero_coag_init addfld', fieldname, unit
1036         end do ! l = ...
1037
1038
1039         return
1040         end subroutine modal_aero_coag_init
1041
1042
1043 !----------------------------------------------------------------------
1044 !----------------------------------------------------------------------
1045         subroutine calc_coag_coef4(   &
1046           dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj,   &
1047           wetddrydi, wetddrydj,   &
1048           temp, presscgs, lunerr, lunout, iok,   &
1049           betaij0, betaij2i, betaij2j, betaij3,   &
1050           betaii0, betaii2, betajj0, betajj2 )
1051 !
1052 !   computes the following coagulation rate "coefficients"
1053 !   for lognormal modes
1054 !
1055 !   self coagulation
1056 !       dNi/dt = - betaii0*Ni*Ni
1057 !       dNj/dt = - betajj0*Nj*Nj
1058 !       dSi/dt = - betaii2*Si*Ni
1059 !       dSj/dt = - betajj2*Sj*Nj
1060 !
1061 !   modei-modej coagulation
1062 !       dNi/dt = - betaij0*Ni*Nj
1063 !       dNj/dt = 0.
1064 !       dSi/dt = - betaij2i*Si*Nj
1065 !       dSj/dt = + betaij2j*Si*Nj
1066 !       dVi/dt  = - betaij3*Vi*Nj == - dVj/dt
1067 !
1068 !   Ni, Nj are number  concentrations (particles/cm3)
1069 !   Si, Sj are surface concentrations (cm2/cm3)
1070 !   Vi, Vj are volume  concentrations (cm3/cm3)
1071 !   rates are (N,S,V)/s
1072 !
1073 !   mode i is the smaller mode (e.g., Aitken)
1074 !   mode j is the larger  mode (e.g., accumulation)
1075 !
1076 !   input arguments
1077 !       dgni, dgnj = DRY median diameters for number distribution (cm)
1078 !       alnsgi, alnsgj = ln of geometric standard deviations (dimensionless)
1079 !       rhopi, rhopj = WET density (g/cm3)
1080 !       wetddrydi, wetddrydj = (WET median diameter)/(DRY median diameter)
1081 !       temp = air temperature (K)
1082 !       presscgs = air pressure (dynes/cm2)
1083 !       lunerr, lunout = logical unit for error or diagnostic output
1084 !
1085 !   output arguments
1086 !       iok = status flag  (+1 = success, 0/negative = failure)
1087 !       beta--- = see above
1088 !
1089
1090         implicit none
1091
1092 !   arguments
1093         integer, intent(in)  :: lunerr, lunout
1094         integer, intent(out) :: iok
1095         real(r4), intent(in) ::   &
1096           dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj,   &
1097           wetddrydi, wetddrydj,   &
1098           temp, presscgs
1099         real(r4), intent(out) ::   &
1100           betaij0, betaij2i, betaij2j, betaij3,   &
1101           betaii0, betaii2, betajj0, betajj2
1102
1103 !   local variables
1104         real(r4) airprs, airtemp
1105         real(r4) dgacc, dgatk, pdensac, pdensat
1106         real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij
1107         real(r4) dp2bar_mks_i, dp2bar_mks_j, dp3bar_mks_i
1108
1109 !   check for reasonable inputs
1110         iok = -1
1111         if ((dgni .lt. 1.e-8) .or. (dgni .gt. 1.)) then
1112             write(lunerr,9100) 'dgni', dgni
1113             return
1114         else if ((dgnj .lt. 1.e-8) .or. (dgnj .gt. 1.)) then
1115             write(lunerr,9100) 'dgnj', dgnj
1116             return
1117 !       else if (dgni .gt. dgnj) then
1118 !           write(lunerr,9110) dgni, dgnj
1119 !           return
1120         else if ((alnsgi .lt. 0.0) .or. (alnsgi .gt. 2.3)) then
1121             write(lunerr,9100) 'alnsgi', alnsgi
1122             return
1123         else if ((alnsgj .lt. 0.0) .or. (alnsgj .gt. 2.3)) then
1124             write(lunerr,9100) 'alnsgj', alnsgj
1125             return
1126         else if ((rhopi .lt. 0.01) .or. (rhopi .gt. 100.)) then
1127             write(lunerr,9100) 'rhopi', rhopi
1128             return
1129         else if ((rhopj .lt. 0.01) .or. (rhopj .gt. 100.)) then
1130             write(lunerr,9100) 'rhopj', rhopj
1131             return
1132         else if ((temp .lt. 10.) .or. (temp .gt. 370.)) then
1133             write(lunerr,9100) 'temp', temp
1134             return
1135         else if ((presscgs .lt. 1.e2) .or. (presscgs .gt. 1.5e6)) then
1136             write(lunerr,9100) 'presscgs', presscgs
1137             return
1138         end if
1139 9100    format( '*** subr. calc_coag_coef_4 - bad value for ', a,   &
1140                 ' = ', 1pe15.5 )
1141 !9110   format( '*** subr. calc_coag_coef_4 - dgni > dgnj -- ',
1142 !     +         2(1pe15.5) )
1143         iok = +1
1144
1145 !
1146 !   define mks variables
1147 !
1148         airtemp = temp
1149 !   dyne/cm2 to pa
1150         airprs = presscgs * 1.0e-1
1151
1152 !   cm to m
1153         dgatk = dgni * 1.0e-2
1154         dgacc = dgnj * 1.0e-2
1155
1156 !   g/cm3 to kg/m3
1157         pdensat = rhopi * 1.0e+3
1158         pdensac = rhopj * 1.0e+3
1159
1160 !
1161 !   call interace to binkowski routines
1162 !
1163         call bink_coag_rates(   &
1164           dgatk, dgacc, alnsgi, alnsgj, pdensat, pdensac,   &
1165           wetddrydi, wetddrydj,   &
1166           airtemp, airprs, lunerr, iok,   &
1167           batat, batac, bacat, bacac, c3ij )
1168
1169 !
1170 !   transfer the batat, batac, bacat, bacac values
1171 !   to the beta--- variables
1172 !
1173
1174 !   self coagulation, number
1175 !       dNi/dt = - betaii0*Ni*Ni
1176 !       dNXi/dt = - batat(1)*NXi*NXi
1177 !   where Ni is (particles/cm3), NXi is (particles/m3)
1178 !   the first-order loss rates are just (1/s) and thus are equal
1179 !       d[ln(Ni)]/dt = d[ln(NXi)]/dt
1180 !   so
1181 !       betaii0 = batat(1) * (NXi/Ni)
1182 !   conversion factors are
1183 !       1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
1184 !       (NXi/Ni) = 1.0e6
1185 !
1186         betaii0 = batat(1) * 1.0e-14
1187         betajj0 = bacac(1) * 1.0e-14
1188
1189 !   self coagulation, surface
1190 !       dSi/dt = - betaii2*Si*Ni
1191 !       dM2Xi/dt = - batat(1)*NXi*NXi
1192 !   where Si is surface in (cm2/cm3) and M2Xi is (surface/pi) in (m2/m3)
1193 !   the first-order loss rates are just (1/s) and thus are equal
1194 !       betaii2 = batat(2) * (NXi/M2Xi) * (NXi/Ni)
1195 !   conversion factors are
1196 !       1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
1197 !       (NXi/Ni) = 1.0e6
1198 !       (M2Xi/NXi) = dp2bar_mks_i/pi, where dp2bar_mks_i is the average
1199 !           Dp**2 in m**2
1200 !
1201         dp2bar_mks_j = (dgnj**2) * exp( 2.0*alnsgj*alnsgj ) * 1.0e-4
1202
1203         dp2bar_mks_i = (dgni**2) * exp( 2.0*alnsgi*alnsgi ) * 1.0e-4
1204         dp3bar_mks_i = (dgni**3) * exp( 4.5*alnsgi*alnsgi ) * 1.0e-6
1205
1206         betaii2 = (batat(2) / dp2bar_mks_i) * 1.0e-14
1207         betajj2 = (bacac(2) / dp2bar_mks_j) * 1.0e-14
1208
1209 !   modei-modej coagulation, number
1210 !       dNi/dt = - betaij0*Ni*Nj
1211 !       dNj/dt = 0.
1212 !   conversions as for self coagulation
1213 !
1214         betaij0  = batac(1) * 1.0e-14
1215
1216 !   modei-modej coagulation, surface
1217 !       dSi/dt = - betaij2i*Si*Nj
1218 !       dSj/dt = + betaij2j*Si*Nj
1219 !   conversions as for self coagulation
1220 !
1221         betaij2i = (batac(2) / dp2bar_mks_i) * 1.0e-14
1222         betaij2j = (bacat(2) / dp2bar_mks_i) * 1.0e-14
1223
1224 !   modei-modej coagulation, volume
1225 !       dVi/dt  = - betaij3*Vi*Nj
1226 !       dM3Xi/dt = - c3ij*NXi*NXj
1227 !   where Vi is (cm3/cm3) and M3Xi is (volume*6/pi) in (m3/m3)
1228 !   the first-order loss rates are just (1/s) and thus are equal
1229 !       betaij3 = c3ij * (NXi/M3Xi) * (NXj/Nj)
1230 !   conversion factors are
1231 !       1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
1232 !       (NXj/Nj) = 1.0e6
1233 !       (M3Xi/NXi) = dp3bar_mks_i*(6/pi), where dp3bar_mks_i is the average
1234 !           Dp**3 in m**3
1235 !
1236         betaij3  = ( c3ij / dp3bar_mks_i ) * 1.0e-14
1237
1238         return
1239         end subroutine calc_coag_coef4
1240
1241
1242
1243 !-----------------------------------------------------------------------
1244         subroutine bink_coag_rates(   &
1245           dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac,   &
1246           wetddrydat, wetddrydac,   &
1247           airtemp, airprs, lunerr, iok,   &
1248           batat, batac, bacat, bacac, c3ij )
1249
1250 !
1251 !   provides interface to F. Binkowski's intracoag and intercoag
1252 !
1253 !   this code was "cut" from F. Binkowski's aero_info_ae3.f & aero_subs3_ae3.f
1254 !
1255 !   computes the following coagulation rate "coefficients"
1256 !
1257 !   self coagulation
1258 !       dNXi/dt  = - (1.e-20*batat(1))*NXi*NXi
1259 !       dNXj/dt  = - (1.e-20*bacac(1))*NXj*NXj
1260 !       dM2Xi/dt = - (1.e-20*batat(2))*NXi*NXi
1261 !       dM2Xj/dt = - (1.e-20*bacac(2))*NXj*NXj
1262 !
1263 !   modei-modej coagulation
1264 !       dNXi/dt  = - (1.e-20*batac(1))*NXi*NXj
1265 !       dNXj/dt  = 0.
1266 !       dM2Xi/dt = - (1.e-20*batac(2))*NXi*NXj
1267 !       dM2Xj/dt = + (1.e-20*bacat(2))*NXi*NXj
1268 !       dM3Xi/dt =     - (1.e-20*c3ij)*NXi*NXj == - dM3Xj/dt
1269 !
1270 !   NXi,  NXj are number concentrations (particles/m3)
1271 !   M2Xi, M2Xj are 2nd moment (surface/pi)  concentrations (m2/m3)
1272 !   M3Xi, M3Xj are 2rd moment (volume*6/pi) concentrations (m3/m3)
1273 !
1274 !   in the above, mode i is aitken mode, mode j is accumulation mode
1275 !
1276 !   *** note that the batat, batac, bacat, bacac, c3ij
1277 !       all must be multiplied by 1.e-20 to undo the 1.e+20 scaling
1278 !       in the real*4 versions of intracoag and intercoag
1279 !
1280 !   input arguments
1281 !       dgatk, dgacc = DRY median diameters for number distribution (m)
1282 !           for aitken (atk or at) and accumulation (acc or ac) modes
1283 !       xxlsgat, xxlsgac = ln of geometric standard deviations (dimensionless)
1284 !       pdensat, pdensac = WET density (kg/m3)
1285 !       wetddrydat, wetddrydac = (WET median diameter)/(DRY median diameter)
1286 !       airtemp = air temperature (K)
1287 !       airprs = air pressure (Pa)
1288 !       lunerr  = logical unit for error output
1289 !
1290 !   output arguments
1291 !       iok = status flag  (+1 = success, 0/negative = failure)
1292 !       batat, batac, bacat, bacac, c3ij = see above
1293 !
1294         implicit none
1295
1296 !   arguments
1297         integer lunerr, iok
1298         real(r4) dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac,   &
1299           wetddrydat, wetddrydac,   &
1300           airtemp, airprs
1301         real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij
1302
1303 ! *** modal geometric mean diameters: [ m ]
1304 !     real dgatk   ! nuclei mode
1305 !     real dgacc   ! accumulation mode
1306 ! *** log of modal geometric standard deviation
1307 !     real xxlsgat  ! aitken mode
1308 !     real xxlsgac  ! accumulation mode
1309 ! *** average modal particle densities  [ kg/m**3 ]
1310 !     real pdensat         ! nuclei mode
1311 !     real pdensac         ! accumulation mode
1312
1313 !     real airtemp   ! air temperature [ k ]
1314 !     real airprs    ! air pressure in [ pa ]
1315
1316
1317 !   local variables
1318
1319       real(r4)   two3
1320       parameter( two3 = 2.0/3.0 )
1321
1322       real(r4)    avo      ! avogadro's constant [ 1/mol ]
1323       parameter ( avo = 6.0221367e23 )
1324
1325       real(r4)    rgasuniv ! universal gas constant [ j/mol-k ]
1326       parameter ( rgasuniv = 8.314510 )
1327
1328       real(r4)    boltz    ! boltzmann's constant [ j / k]
1329       parameter ( boltz = rgasuniv / avo )
1330
1331       real(r4)    p0       !  starting standard surface pressure [ pa ]
1332       parameter ( p0 = 101325.0 )
1333
1334       real(r4)    t0       !  starting standard surface temperature [ k ]
1335       parameter ( t0 = 288.15 )
1336
1337       real(r4) xlm    ! atmospheric mean free path [ m ]
1338       real(r4) amu    ! atmospheric dynamic viscosity [ kg/m s ]
1339
1340       real(r4) kfm, knc, lamda, sqrt_temp
1341
1342
1343 ! fsb calculate the square root of the ambient
1344 !      temperature for later use
1345
1346 ! *** calculate mean free path [ m ]:
1347 ! *** 6.6328e-8 is the sea level values given in table i.2.8
1348 ! *** on page 10 of u.s. standard atmosphere 1962
1349       xlm = 6.6328e-8 * p0 * airtemp / ( t0 * airprs )
1350
1351 ! ***   calculate dynamic viscosity [ kg m**-1 s**-1 ]:
1352 ! *** u.s. standard atmosphere 1962 page 14 expression
1353 !     for dynamic viscosity is:
1354 !     dynamic viscosity =  beta * t * sqrt(t) / ( t + s)
1355 !     where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
1356       sqrt_temp = sqrt( airtemp)
1357       amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )
1358
1359
1360
1361 ! *** coagulation
1362 ! *** set up coagulation rates
1363
1364 ! *** moment independent factors
1365             knc = two3 * boltz *  airtemp / amu
1366             lamda = xlm
1367
1368 ! ***  calculate the coagulation coefficients for use
1369 !      in the  aitken (nuclei) & accumulation modes
1370 ! ***  with gauss-hermite numerical quadrature
1371 !      using 10 abscissas.
1372
1373 ! *** aitken - aitken mode coagulation
1374       kfm = sqrt( 3.0 * boltz * airtemp / pdensat )
1375       call  intracoag_gh(lamda,   &
1376                          kfm, knc,   &
1377                          dgatk,   &
1378                          xxlsgat,   &
1379                          wetddrydat,   &
1380                          batat(2),   &
1381                          batat(1) )
1382
1383 ! *** accumulation - accumulation mode coagulation
1384       kfm = sqrt( 3.0 * boltz * airtemp / pdensac )
1385       call  intracoag_gh(lamda,   &
1386                          kfm, knc,   &
1387                          dgacc,   &
1388                          xxlsgac,   &
1389                          wetddrydac,   &
1390                          bacac(2),   &
1391                          bacac(1) )
1392
1393 ! *** aitken accumulation mode coagulation
1394       bacat(1) = 0.0 ! not used
1395       kfm  = sqrt( 6.0 * boltz * airtemp  /   &
1396                              ( pdensat + pdensac ) )
1397       call  intercoag_gh(lamda,   &
1398                          kfm, knc,   &
1399                          dgatk, dgacc ,   &
1400                          xxlsgat, xxlsgac ,   &
1401                          wetddrydat, wetddrydac,   &
1402                          batac(2),   &
1403                          bacat(2),   &
1404                          batac(1),   &
1405                          c3ij          )
1406
1407 !     c30atac = c3ij * cblk( vac0 ) * cblk( vat0 )
1408 !     loss =  c30atac / cblk( vat3 )
1409
1410       return
1411       end subroutine bink_coag_rates
1412
1413
1414
1415 ! ---------------------------------------------------------------------
1416 ! fsb subrs to do gauss-hermite numerical quadrature
1417
1418       subroutine intracoag_gh( lamda, kfm, knc,   &
1419            dg, xlnsig, wetddryd,   &
1420            quads11, quadn11)
1421
1422 !  fsb this version is for intramodal coagulation for number
1423 !       and second moment
1424
1425 ! fsb this version runs in real*4 arithmetic
1426
1427 ! *** this version calculates the coagulation coefficients
1428 !     using the harmonic mean approach for both fm and nc cases.
1429 ! *** does gauss-hermite quadrature for intra-modal
1430 !      coagulation integrals for 2nd moment
1431 !     for a lognormal distribution
1432 !      defined by dg,xlnsig,
1433 ! *** dg and xlnsig are the geometric mean diameters (meters)
1434 !      and the logarithms of the
1435 !     geometric standard deviations (dimensionless)
1436 !     whose meaning is defined below at the end of the routine
1437 !     ghxi, ghwi are the gauss-hermite weights and n is one-half the
1438 !     number of abscissas, since an even number of abscissas is used
1439 !
1440 !.......................................................................
1441 !   (following comments added by rc_easter)
1442 !
1443 !   computes the following coagulation rate "coefficients"
1444 !   for intramodal coagulation
1445 !       dNX/dt  = - (1.e-20*quadn11)*NX1*NX1
1446 !       dM2X/dt = - (1.e-20*quads11)*NX1*NX1
1447 !
1448 !   NX  is the mode's number concentration (particles/m3)
1449 !   M2X is the mode's 2nd moment (surface/pi) concentration (m2/m3)
1450 !
1451 !   input arguments
1452 !       lamda = mean free path (m)
1453 !       kfm, knc = constants used in free-molecular and near-continuum
1454 !           calculations (see subr bink_coag_rates)
1455 !       dg = DRY median diameter for number distribution (m)
1456 !       xlnsig = ln of geometric standard deviation (dimensionless)
1457 !       wetddryd = (WET median diameter)/(DRY median diameter)
1458 !
1459 !   output arguments
1460 !       quads11, quadn11 = coagulation rate coefficients (see above)
1461 !
1462 !   *** note that the quads11, quadn11
1463 !       are all scaled by 1.e+20 to avoid underflow, and they
1464 !       must be multiplied by 1.e-20 when they are applied
1465 !
1466 !   *** wetddryd1 was added because MIRAGE treats N/S/V of
1467 !       the DRY size distribution, so the quadrature should be done
1468 !       over the DRY size distribution.  However, the coagulation kernel
1469 !       must be computed using actual (WET) particle sizes
1470 !
1471 !.......................................................................
1472 !
1473       implicit none
1474
1475       integer i,j
1476
1477       real(r4) lamda ! mean free path
1478       real(r4) kfm, knc
1479       real(r4) dg, xlnsig, wetddryd
1480       real(r4) quads11, quadn11
1481       real(r4) pi
1482       parameter( pi = 3.14159265358979)
1483       real(r4) two3rds
1484       parameter( two3rds = 2.0d0 /  3.0d0 )
1485       real(r4) sqrt2
1486       parameter(sqrt2 = 1.41421356237309 )
1487       real(r4) sum1sfm, sum2sfm, sum1nfm, sum2nfm
1488       real(r4) sum1snc, sum2snc, sum1nnc, sum2nnc
1489       real(r4) xi, wxi, xf, dp1p,dp1m,dp1psq,dp1msq, dp1pwet, dp1mwet
1490       real(r4) v1p,v1m, a2p,a2m,v2p,v2m
1491       real(r4) yi,wyi,yf,dp2p,dp2m,dp2psq,dp2msq, dp2pwet, dp2mwet
1492       real(r4) dspp,dsmp,dspm, dsmm
1493       real(r4) bppfm,bmpfm,bpmfm,bmmfm
1494       real(r4) bppnc,bmpnc,bpmnc,bmmnc
1495       real(r4) xx1, xx2
1496       real(r4) xbsfm, xbsnc, xbnfm, xbnnc
1497       real(r4) betafm, betanc
1498
1499       real(r4)   a               ! approx cunningham corr. factor
1500       parameter( a = 1.246d0 )
1501
1502       real(r4)   twoa
1503       parameter( twoa = 2.0d0 * a )
1504
1505 ! *** has a fixed number of gauss-herimite abscissas (n)
1506       integer n    ! one-half the number of abscissas
1507       parameter ( n = 5 )
1508       real(r4), save :: ghxi(n) ! gauss-hermite abscissas
1509       real(r4), save :: ghwi(n) ! gauss-hermite weights
1510
1511 ! ** values from table 25.10 (page 924) of abramowitz and stegun,
1512 !    handbook of mathematical functions, national bureau of standards,
1513 !  december 1965.
1514
1515 !    breaks in number to facilitate comparison with printed table
1516
1517 ! *** tests show that 10 point is adquate.
1518
1519       data ghxi/0.342901327223705,   &
1520                 1.036610829789514,   &
1521                 1.756683649299882,   &
1522                 2.532731674232790,   &
1523                 3.436159118837738/
1524
1525       data ghwi/6.108626337353d-1,   &
1526                 2.401386110823d-1,   &
1527                 3.387439445548d-2,   &
1528                 1.343645746781d-3,   &
1529                 7.640432855233d-6/
1530
1531 ! *** the following expressions are from binkowski & shanker
1532 !     jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
1533 !     december 20, 1995
1534 ! ***  for free molecular eq. a5
1535         betafm(xx1, xx2) = kfm *   &
1536              sqrt(1.0 / xx1**3  + 1.0 / xx2**3 ) * (xx1 + xx2)**2
1537
1538 ! ***  for near continuum  eq. a6
1539         betanc(xx1, xx2) =  knc * (xx1 + xx2) *   &
1540                              ( 1.0d0 / xx1 + 1.0d0 / xx2  +   &
1541                            twoa * lamda * ( 1.0d0 / xx1 ** 2   &
1542                                           + 1.0d0 / xx2 **2 ) )
1543
1544
1545       sum1sfm = 0.0
1546       sum1snc = 0.0
1547
1548       sum1nfm = 0.0
1549       sum1nnc = 0.0
1550       do 1 i=1,n
1551
1552         sum2sfm = 0.0
1553         sum2snc = 0.0
1554         sum2nfm = 0.0
1555         sum2nnc = 0.0
1556
1557         xi = ghxi(i)
1558         wxi = ghwi(i)
1559         xf = exp( sqrt2 * xi *xlnsig)
1560         dp1p = dg*xf
1561         dp1m = dg/xf
1562         dp1psq = dp1p*dp1p
1563         dp1msq = dp1m*dp1m
1564         v1p = dp1p*dp1psq
1565         v1m = dp1m*dp1msq
1566
1567         dp1pwet = dp1p * wetddryd
1568         dp1mwet = dp1m * wetddryd
1569
1570       do 11 j=1,n
1571         yi = ghxi(j)
1572         wyi = ghwi(j)
1573         yf = exp( sqrt2 * yi * xlnsig)
1574         dp2p = dg*yf
1575         dp2m = dg/yf
1576         dp2psq = dp2p*dp2p
1577         dp2msq = dp2m*dp2m
1578         a2p = dp2psq
1579         a2m = dp2msq
1580         v2p =  dp2p*dp2psq
1581         v2m =dp2m*dp2msq
1582         dspp = 0.5*(v1p+v2p)**two3rds - a2p
1583         dsmp = 0.5*(v1m+v2p)**two3rds - a2p
1584         dspm = 0.5*(v1p+v2m)**two3rds - a2m
1585         dsmm = 0.5*(v1m+v2m)**two3rds - a2m
1586
1587         dp2pwet = dp2p * wetddryd
1588         dp2mwet = dp2m * wetddryd
1589
1590 !   scale by 1.0e+20 to avoid underflow
1591         bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
1592         bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
1593         bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
1594         bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20
1595
1596         bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
1597         bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
1598         bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
1599         bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20
1600
1601         sum2sfm = sum2sfm + wyi*(dspp * bppfm + dspm * bpmfm   &
1602                      +   dsmp * bmpfm + dsmm * bmmfm )
1603
1604         sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)
1605
1606         sum2snc = sum2snc + wyi*(dspp * bppnc + dspm * bpmnc   &
1607                      +   dsmp * bmpnc + dsmm * bmmnc )
1608         sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)
1609
1610    11 continue
1611       sum1sfm = sum1sfm + wxi * sum2sfm
1612       sum1nfm = sum1nfm + wxi * sum2nfm
1613
1614       sum1snc = sum1snc + wxi * sum2snc
1615       sum1nnc = sum1nnc + wxi * sum2nnc
1616
1617     1 continue
1618
1619       xbsfm   = -sum1sfm  / pi
1620       xbsnc   = -sum1snc  / pi
1621
1622 !     quads11 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
1623       quads11 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
1624                  * min(xbsfm,xbsnc)
1625
1626 ! *** quads11 is the intra-modal coagulation term for 2nd moment
1627
1628       xbnfm   = 0.5 * sum1nfm  / pi
1629       xbnnc   = 0.5 * sum1nnc  / pi
1630
1631
1632 !     quadn11 =  xbnfm * xbnnc / ( xbnfm + xbnnc )
1633       quadn11 =  ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) )   &
1634                  * min(xbnfm,xbnnc)
1635
1636 ! *** quadn11 is the intra-modal coagulation term for number
1637
1638
1639       return
1640       end subroutine intracoag_gh
1641
1642
1643 ! ---------------------------------------------------------------------
1644        subroutine intercoag_gh( lamda, kfm, knc, dg1, dg2,   &
1645            xlnsig1, xlnsig2,   &
1646            wetddryd1, wetddryd2,   &
1647            quads12, quads21, quadn12, quadv12 )
1648
1649 !  fsb this version is for intermodal coagulation for number,
1650 !       second, and third moments
1651 ! fsb this version runs in real*4 arithmetic
1652
1653 ! *** this version calculates the coagulation coefficients
1654 !     using the harmonic mean approach for both fm and nc cases.
1655 ! *** does gauss-hermite quadrature for inter-modal
1656 !      coagulation integrals for 2nd moment
1657 !     for two lognormal distributions
1658 !      defined by dg1,xlnsig1, dg2,xlnsig2
1659 ! *** dg and xlnsig are the geometric mean diameters (meters)
1660 !      and the logarithms of the
1661 !     geometric standard deviations (dimensionless)
1662 !     whose meaning is defined below at the end of the routine
1663 !     ghxi, ghwi are the gauss-hermite weights and n is one-half the
1664 !     number of abscissas, since an even number of abscissas is used
1665 !
1666 !.......................................................................
1667 !   (following comments added by rc_easter)
1668 !
1669 !   computes the following coagulation rate "coefficients"
1670 !   for mode1-mode2 coagulation
1671 !       dNX1/dt  = - (1.e-20*quadn12)*NX1*NX2
1672 !       dNX2/dt  = 0.
1673 !       dM2X1/dt = - (1.e-20*quads12)*NX1*NX2
1674 !       dM2X2/dt = + (1.e-20*quads21)*NX1*NX2
1675 !       dM3X1/dt = - (1.e-20*quads12)*NX1*NX2 == - dM3X2/dt
1676 !
1677 !   NX1,  NX2 are number concentrations (particles/m3)
1678 !   M2X1, M2X2 are 2nd moment (surface/pi)  concentrations (m2/m3)
1679 !   M3X1, M3X2 are 2rd moment (volume*6/pi) concentrations (m3/m3)
1680 !
1681 !   in the above, mode 1 is aitken mode, mode 2 is accumulation mode
1682 !
1683 !   input arguments
1684 !       lamda = mean free path (m)
1685 !       kfm, knc = constants used in free-molecular and near-continuum
1686 !           calculations (see subr bink_coag_rates)
1687 !       dg1, dg2 = DRY median diameters for number distribution (m)
1688 !           for aitken and accumulation modes
1689 !       xlnsig1, xlnsig2 = ln of geometric standard deviations (dimensionless)
1690 !       wetddryd1, wetddryd2 = (WET median diameter)/(DRY median diameter)
1691 !
1692 !   output arguments
1693 !       quads12, quads21, quadn12, quadv12 = coagulation rate
1694 !           coefficients (see above)
1695 !
1696 !   *** note that the quads12, quads21, quadn12, quadv12
1697 !       are all scaled by 1.e+20 to avoid underflow, and they
1698 !       must be multiplied by 1.e-20 when they are applied
1699 !
1700 !   *** wetddryd1/2 were added because MIRAGE treats N/S/V of
1701 !       the DRY size distribution, so the quadrature should be done
1702 !       over the DRY size distribution.  However, the coagulation kernel
1703 !       must be computed using actual (WET) particle sizes
1704 !
1705 !.......................................................................
1706
1707       implicit none
1708
1709       integer i,j
1710
1711       real(r4) lamda ! mean free path
1712       real(r4) kfm, knc
1713       real(r4) dg1, xlnsig1, dg2, xlnsig2
1714       real(r4) wetddryd1, wetddryd2
1715       real(r4) quads12, quads21, quadn12, quadv12
1716       real(r4) pi
1717       parameter( pi = 3.14159265358979)
1718       real(r4) two3rds
1719       parameter( two3rds = 2.0d0 /  3.0d0 )
1720       real(r4) sqrt2
1721       parameter(sqrt2 = 1.41421356237309 )
1722       real(r4) sum1s12fm, sum1s21fm, sum2s12fm, sum2s21fm
1723       real(r4) sum1nfm, sum2nfm
1724       real(r4) sum1vfm, sum2vfm
1725       real(r4) sum1s12nc, sum1s21nc, sum2s12nc, sum2s21nc
1726       real(r4) sum1nnc, sum2nnc
1727       real(r4) sum1vnc, sum2vnc
1728       real(r4) xi, wxi,xf, dp1p, dp1m, dp1psq, dp1msq, dp1pwet, dp1mwet
1729       real(r4) a1p, a1m, v1p, v1m
1730       real(r4) a2p, a2m, v2p, v2m
1731       real(r4) yi, wyi, yf, dp2p, dp2m, dp2psq, dp2msq, dp2pwet, dp2mwet
1732       real(r4) dspp, dsmp, dspm, dsmm
1733       real(r4) bppfm, bmpfm, bpmfm, bmmfm
1734       real(r4) bppnc, bmpnc, bpmnc, bmmnc
1735       real(r4) xx1, xx2
1736       real(r4) xbsfm, xbsnc, xbnfm, xbnnc, xbvfm, xbvnc
1737       real(r4) betafm, betanc
1738
1739       real(r4)   a               ! approx cunningham corr. factor
1740       parameter( a = 1.246d0 )
1741
1742       real(r4)   twoa
1743       parameter( twoa = 2.0d0 * a )
1744
1745 ! *** has a fixed number of gauss-herimite abscissas (n)
1746       integer n    ! one-half the number of abscissas
1747       parameter ( n = 5 )
1748       real(r4), save :: ghxi(n) ! gauss-hermite abscissas
1749       real(r4), save :: ghwi(n) ! gauss-hermite weights
1750
1751 ! ** values from table 25.10 (page 924) of abramowitz and stegun,
1752 !    handbook of mathematical functions, national bureau of standards,
1753 !  december 1965.
1754
1755 !    breaks in number to facilitate comparison with printed table
1756
1757 ! *** tests show that 10 point is adquate.
1758
1759       data ghxi/0.342901327223705,   &
1760                 1.036610829789514,   &
1761                 1.756683649299882,   &
1762                 2.532731674232790,   &
1763                 3.436159118837738/
1764
1765       data ghwi/6.108626337353d-1,   &
1766                 2.401386110823d-1,   &
1767                 3.387439445548d-2,   &
1768                 1.343645746781d-3,   &
1769                 7.640432855233d-6/
1770
1771
1772 ! *** the following expressions are from binkowski & shanker
1773 !     jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
1774 !     december 20, 1995
1775
1776 ! ***  for free molecular eq. a5
1777         betafm(xx1, xx2) = kfm *   &
1778              sqrt(1.0 / xx1**3  + 1.0 / xx2**3 ) * (xx1 + xx2)**2
1779
1780 ! ***  for near continuum  eq. a6
1781         betanc(xx1, xx2) =  knc * (xx1 + xx2) *   &
1782                              ( 1.0d0 / xx1 + 1.0d0 / xx2  +   &
1783                            twoa * lamda * ( 1.0d0 / xx1 ** 2   &
1784                                           + 1.0d0 / xx2 **2 ) )
1785
1786       sum1s12fm = 0.0
1787       sum1s12nc = 0.0
1788       sum1s21fm = 0.0
1789       sum1s21nc = 0.0
1790       sum1vnc = 0.0
1791       sum1vfm = 0.0
1792       sum1nfm = 0.0
1793       sum1nnc = 0.0
1794       do 1 i=1,n
1795
1796         sum2s12fm = 0.0
1797         sum2s12nc = 0.0
1798         sum2s21fm = 0.0
1799         sum2s21nc = 0.0
1800         sum2nfm = 0.0
1801         sum2nnc = 0.0
1802         sum2vnc = 0.0
1803         sum2vfm = 0.0
1804         xi = ghxi(i)
1805         wxi = ghwi(i)
1806         xf = exp( sqrt2 * xi *xlnsig1)
1807         dp1p = dg1*xf
1808         dp1m = dg1/xf
1809         dp1psq = dp1p*dp1p
1810         dp1msq = dp1m*dp1m
1811         a1p = dp1psq
1812         a1m = dp1msq
1813         v1p = dp1p*dp1psq
1814         v1m = dp1m*dp1msq
1815
1816         dp1pwet = dp1p * wetddryd1
1817         dp1mwet = dp1m * wetddryd1
1818
1819       do 11 j=1,n
1820         yi  = ghxi(j)
1821         wyi = ghwi(j)
1822         yf = exp( sqrt2 * yi * xlnsig2)
1823         dp2p = dg2*yf
1824         dp2m = dg2/yf
1825         dp2psq = dp2p*dp2p
1826         dp2msq = dp2m*dp2m
1827         a2p  = dp2psq
1828         a2m  = dp2msq
1829         v2p  =  dp2p*dp2psq
1830         v2m  = dp2m*dp2msq
1831         dspp = (v1p+v2p)**two3rds - a2p
1832         dsmp = (v1m+v2p)**two3rds - a2p
1833         dspm = (v1p+v2m)**two3rds - a2m
1834         dsmm = (v1m+v2m)**two3rds - a2m
1835
1836         dp2pwet = dp2p * wetddryd2
1837         dp2mwet = dp2m * wetddryd2
1838
1839 !   scale by 1.0e+20 to avoid underflow
1840         bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
1841         bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
1842         bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
1843         bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20
1844
1845         bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
1846         bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
1847         bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
1848         bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20
1849
1850
1851         sum2s12fm = sum2s12fm + wyi*(a1p * bppfm + a1p * bpmfm   &
1852                      +   a1m * bmpfm + a1m * bmmfm )
1853
1854         sum2s21fm = sum2s21fm + wyi*(dspp * bppfm + dspm * bpmfm   &
1855                      +   dsmp * bmpfm + dsmm * bmmfm )
1856
1857
1858         sum2s12nc = sum2s12nc + wyi*(a1p * bppnc + a1p * bpmnc   &
1859                      +   a1m * bmpnc + a1m * bmmnc )
1860
1861         sum2s21nc = sum2s21nc + wyi*(dspp * bppnc + dspm * bpmnc   &
1862                      +   dsmp * bmpnc + dsmm * bmmnc )
1863
1864         sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)
1865
1866         sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)
1867
1868         sum2vfm = sum2vfm + wyi*(v1p*(bppfm + bpmfm) +   &
1869                                  v1m*(bmpfm + bmmfm) )
1870
1871         sum2vnc = sum2vnc + wyi*(v1p*(bppnc + bpmnc) +   &
1872                                  v1m*(bmpnc + bmmnc) )
1873
1874    11 continue
1875
1876       sum1s12fm = sum1s12fm + wxi * sum2s12fm
1877       sum1s21fm = sum1s21fm + wxi * sum2s21fm
1878       sum1nfm   = sum1nfm + wxi * sum2nfm
1879       sum1vfm   = sum1vfm + wxi * sum2vfm
1880
1881       sum1s12nc = sum1s12nc + wxi * sum2s12nc
1882       sum1s21nc = sum1s21nc + wxi * sum2s21nc
1883       sum1nnc   = sum1nnc + wxi * sum2nnc
1884       sum1vnc   = sum1vnc + wxi * sum2vnc
1885
1886     1 continue
1887
1888 ! *** second moment intermodal coagulation coefficients
1889
1890 ! fsb note: the transfer of second moment is not symmetric.
1891 !     see equations a3 & a4 of binkowski & shankar (1995)
1892
1893 ! ***  to accumulation mode from aitken mode
1894
1895       xbsfm   = sum1s21fm  / pi
1896       xbsnc   = sum1s21nc  / pi
1897
1898 !     quads21 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
1899       quads21 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
1900                  * min(xbsfm,xbsnc)
1901
1902 ! *** from aitken mode to accumulation mode
1903
1904       xbsfm   = sum1s12fm  / pi
1905       xbsnc   = sum1s12nc  / pi
1906
1907 !     quads12 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
1908       quads12 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
1909                  * min(xbsfm,xbsnc)
1910
1911
1912
1913       xbnfm   = sum1nfm  / pi
1914       xbnnc   = sum1nnc  / pi
1915
1916 !     quadn12 =  xbnfm * xbnnc / ( xbnfm + xbnnc )
1917       quadn12 =  ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) )   &
1918                  * min(xbnfm,xbnnc)
1919
1920 ! *** quadn12 is the intermodal coagulation coefficient for number
1921
1922
1923        xbvfm = sum1vfm / pi
1924        xbvnc = sum1vnc / pi
1925
1926 !     quadv12 = xbvfm * xbvnc / ( xbvfm + xbvnc )
1927       quadv12 =  ( max(xbvfm,xbvnc) / ( xbvfm + xbvnc ) )   &
1928                  * min(xbvfm,xbvnc)
1929
1930 ! *** quadv12 is the intermodal coagulation coefficient for 3rd moment
1931
1932
1933       return
1934       end subroutine intercoag_gh
1935
1936
1937 !----------------------------------------------------------------------
1938 !----------------------------------------------------------------------
1939       subroutine getcoags_wrapper_f(              &
1940           airtemp, airprs,                        &
1941           dgatk, dgacc,                           &
1942           sgatk, sgacc,                           &
1943           xxlsgat, xxlsgac,                       &
1944           pdensat, pdensac,                       &
1945           betaij0, betaij2i, betaij2j, betaij3,   &
1946           betaii0, betaii2, betajj0, betajj2      )
1947 !
1948 ! interface to subr. getcoags
1949 !
1950 ! interface code adapted from subr. aeroproc of cmaq v4.6,
1951 !     with some of the parameter values from module aero_info_ae4
1952 !
1953       implicit none
1954
1955 ! *** arguments
1956
1957       real(r8), intent(in) :: airtemp  ! air temperature [ k ]
1958       real(r8), intent(in) :: airprs   ! air pressure in [ pa ]
1959
1960       real(r8), intent(in) :: dgatk    ! aitken mode geometric mean diameter [m]
1961       real(r8), intent(in) :: dgacc    ! accumulation mode geometric mean diam [m]
1962
1963       real(r8), intent(in) :: sgatk    ! aitken mode geometric standard deviation
1964       real(r8), intent(in) :: sgacc    ! accumulation mode geometric standard deviation
1965
1966       real(r8), intent(in) :: xxlsgat  ! natural log of geometric standard
1967       real(r8), intent(in) :: xxlsgac  !  deviations
1968
1969       real(r8), intent(in) :: pdensat  ! aitken mode particle density [ kg / m**3 ]
1970       real(r8), intent(in) :: pdensac  ! accumulation mode density [ kg / m**3 ]
1971
1972       real(r8), intent(out) :: betaij0, betaij2i, betaij2j, betaij3,   &
1973                                betaii0, betaii2,  betajj0,  betajj2
1974
1975
1976 ! *** local parameters
1977       real(r8), parameter :: p0 = 101325.0         !  standard surface pressure [ pa ]
1978       real(r8), parameter :: t0 = 288.15           !  standard surface temperature [ k ]
1979       real(r8), parameter :: avo = 6.0221367e23    ! avogadro's constant [ 1/mol ]
1980       real(r8), parameter :: rgasuniv = 8.314510   ! universal gas constant [ j/mol-k ]
1981       real(r8), parameter :: boltz = rgasuniv/avo  ! boltzmann's constant [ j / k]
1982       real(r8), parameter :: two3 = 2.0/3.0
1983
1984 ! *** local variables
1985       real(r8) amu            ! atmospheric dynamic viscosity [ kg/m s ]
1986       real(r8) sqrt_temp      ! square root of ambient temperature
1987       real(r8) lamda          ! mean free path [ m ]
1988
1989 ! *** intramodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
1990       real(r8)    batat( 2 )  ! aitken mode
1991       real(r8)    bacac( 2 )  ! accumulation mode
1992 ! *** intermodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
1993       real(r8)    batac( 2 )  ! aitken to accumulation
1994       real(r8)    bacat( 2 )  ! accumulation from aitken
1995 ! *** intermodal coagulation rate [ m**3/s ] ( 3rd moment )
1996       real(r8)    c3ij        ! aitken to accumulation
1997 ! *** 3rd moment intermodal transfer rate by coagulation
1998       real(r8)    c30atac     ! aitken to accumulation
1999
2000 ! *** near continnuum regime (independent of mode)
2001       real(r8)    knc         ! knc = two3 * boltz *  airtemp / amu
2002 ! *** free molecular regime (depends upon modal density)
2003       real(r8)    kfmat       ! kfmat = sqrt(3.0*boltz*airtemp/pdensat)
2004       real(r8)    kfmac       ! kfmac = sqrt(3.0*boltz*airtemp/pdensac)
2005       real(r8)    kfmatac     ! kfmatac = sqrt( 6.0 * boltz * airtemp /
2006                               !                ( pdensat + pdensac ) )
2007
2008       real(r8)    dumacc2, dumatk2, dumatk3
2009
2010
2011
2012       sqrt_temp = sqrt( airtemp)
2013
2014 ! *** calculate mean free path [ m ]:
2015 !     6.6328e-8 is the sea level value given in table i.2.8
2016 !     on page 10 of u.s. standard atmosphere 1962
2017       lamda = 6.6328e-8 * p0 * airtemp  / ( t0 * airprs )
2018
2019 ! *** calculate dynamic viscosity [ kg m**-1 s**-1 ]:
2020 !     u.s. standard atmosphere 1962 page 14 expression
2021 !     for dynamic viscosity is:
2022 !     dynamic viscosity =  beta * t * sqrt(t) / ( t + s)
2023 !     where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
2024       amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )
2025
2026 ! *** coagulation
2027 !     calculate coagulation coefficients using a method dictated by
2028 !     the value of fastcoag_flag.  if true, the computationally-
2029 !     efficient getcoags routine is used.  if false, the more intensive
2030 !     gauss-hermite numerical quadrature method is used.  see section
2031 !     2.1 of bhave et al. (2004) for further discussion.
2032
2033 ! *** calculate term used in equation a6 of binkowski & shankar (1995)
2034       knc      = two3 * boltz *  airtemp / amu
2035 ! *** calculate terms used in equation a5 of binkowski & shankar (1995)
2036       kfmat    = sqrt( 3.0 * boltz * airtemp / pdensat )
2037       kfmac    = sqrt( 3.0 * boltz * airtemp / pdensac )
2038       kfmatac  = sqrt( 6.0 * boltz * airtemp / ( pdensat + pdensac ) )
2039
2040 ! *** transfer of number to accumulation mode from aitken mode is zero
2041       bacat(1) = 0.0
2042
2043 ! *** calculate intermodal and intramodal coagulation coefficients
2044 !     for zeroth and second moments, and intermodal coagulation
2045 !     coefficient for third moment
2046         call getcoags( lamda, kfmatac, kfmat, kfmac, knc,   &
2047                        dgatk,   dgacc,   sgatk,   sgacc,     &
2048                        xxlsgat,  xxlsgac,     &
2049                        batat(2), batat(1), bacac(2), bacac(1),   &
2050                        batac(2), bacat(2), batac(1), c3ij )
2051
2052 ! convert from the "cmaq" coag rate parameters
2053 ! to the "mirage2" parameters
2054         dumacc2 = ( (dgacc**2) * exp( 2.0*xxlsgac*xxlsgac ) )
2055         dumatk2 = ( (dgatk**2) * exp( 2.0*xxlsgat*xxlsgat ) )
2056         dumatk3 = ( (dgatk**3) * exp( 4.5*xxlsgat*xxlsgat ) )
2057
2058         betaii0  = max( 0.0_r8, batat(1) )
2059         betajj0  = max( 0.0_r8, bacac(1) )
2060         betaij0  = max( 0.0_r8, batac(1) )
2061         betaij3  = max( 0.0_r8, c3ij / dumatk3 )
2062
2063         betajj2  = max( 0.0_r8, bacac(2) / dumacc2 )
2064         betaii2  = max( 0.0_r8, batat(2) / dumatk2 )
2065         betaij2i = max( 0.0_r8, batac(2) / dumatk2 )
2066         betaij2j = max( 0.0_r8, bacat(2) / dumatk2 )
2067
2068
2069       return
2070       end subroutine getcoags_wrapper_f
2071
2072
2073
2074 ! //////////////////////////////////////////////////////////////////
2075 !  subroutine getcoags calculates the coagulation rates using a new
2076 !     approximate algorithm for the 2nd moment.  the 0th and 3rd moments
2077 !     are done by analytic expressions from whitby et al. (1991).  the
2078 !     correction factors are also similar to those from whitby et al.
2079 !     (1991), but are derived from the gauss-hermite numerical
2080 !     quadratures used by binkowski and roselle (2003).
2081 !
2082 !     called from aerostep as:
2083 !     call getcoags( lamda, kfmatac, kfmat, kfmac, knc,
2084 !                    dgat,dgac, sgatk, sgacc, xxlsgat,xxlsgac,
2085 !                    batat(2), batat(1), bacac(2), bacac(1),
2086 !                    batac(2), bacat(2), batac(1), c3ij )
2087 !     where all input and outputs are real*8
2088 !
2089 !  revision history:
2090 !   fsb 08/25/03 coded by dr. francis s. binkowksi
2091 !
2092 !   fsb 08/25/04 added in-line documentation
2093 !
2094 !   rce 04/15/2007
2095 !       code taken from cmaq v4.6 code; converted to f90;
2096 !       added "intent" to subr arguments;
2097 !       renamed "r4" & "r8" variables to "rx4" & "rx8";
2098 !       changed "real*N" declarations to "real(rN)" (N = 4 or 8)
2099 !
2100 !  references:
2101 !   1. whitby, e. r., p. h. mcmurry, u. shankar, and f. s. binkowski,
2102 !   modal aerosol dynamics modeling, rep. 600/3-91/020, atmospheric
2103 !   research and exposure assessment laboratory,
2104 !   u.s. environmental protection agency, research triangle park, n.c.,
2105 !   (ntis pb91-161729/as), 1991
2106 !
2107 !   2. binkowski, f.s. an u. shankar, the regional particulate matter
2108 !   model 1. model decsription and preliminary results, journal of
2109 !   geophysical research, 100, d12, pp 26,191-26,209,
2110 !   december 20, 1995.
2111 !
2112 !   3. binkowski, f.s. and s.j. roselle, models-3 community
2113 !      multiscale air quality (cmaq) model aerosol component 1:
2114 !      model description.  j. geophys. res., vol 108, no d6, 4183
2115 !      doi:10.1029/2001jd001409, 2003.
2116
2117
2118       subroutine getcoags( lamda, kfmatac, kfmat, kfmac, knc,   &
2119                            dgatk, dgacc, sgatk, sgacc, xxlsgat,xxlsgac,   &
2120                            qs11, qn11, qs22, qn22,   &
2121                            qs12, qs21, qn12, qv12 )
2122
2123       implicit none
2124
2125       real(r8), intent(in) ::  lamda     ! mean free path [ m ]
2126
2127 ! *** coefficients for free molecular regime
2128       real(r8), intent(in) ::  kfmat     ! aitken mode
2129       real(r8), intent(in) ::  kfmac     ! accumulation mode
2130       real(r8), intent(in) ::  kfmatac   ! aitken to accumulation mode
2131
2132       real(r8), intent(in) ::  knc   ! coefficient for near continnuum regime
2133
2134 ! *** modal geometric mean diameters: [ m ]
2135       real(r8), intent(in) :: dgatk          ! aitken mode
2136       real(r8), intent(in) :: dgacc          ! accumulation mode
2137
2138 ! *** modal geometric standard deviation
2139       real(r8), intent(in) :: sgatk          ! atken mode
2140       real(r8), intent(in) :: sgacc          ! accumulation mode
2141
2142 ! *** natural log of modal geometric standard deviation
2143       real(r8), intent(in) :: xxlsgat         ! aitken mode
2144       real(r8), intent(in) :: xxlsgac         ! accumulation mode
2145
2146 ! *** coagulation coefficients
2147       real(r8), intent(out) :: qs11, qn11, qs22, qn22,   &
2148                                qs12, qs21, qn12, qv12
2149
2150       integer ibeta, n1, n2a, n2n ! indices for correction factors
2151
2152       real(r8)  i1fm_at
2153       real(r8)  i1nc_at
2154       real(r8)  i1_at
2155
2156       real(r8)  i1fm_ac
2157       real(r8)  i1nc_ac
2158       real(r8)  i1_ac
2159
2160       real(r8)  i1fm
2161       real(r8)  i1nc
2162       real(r8)  i1
2163
2164       real(r8) constii
2165
2166       real(r8)    kngat, kngac
2167       real(r8)    one, two, half
2168        parameter( one = 1.0d0, two = 2.0d0, half = 0.5d0 )
2169       real(r8)    a
2170 !       parameter( a = 2.492d0)
2171       parameter( a = 1.246d0)
2172       real      two3rds
2173        parameter( two3rds = 2.d0 / 3.d0)
2174
2175       real(r8)   sqrttwo  !  sqrt(two)
2176       real(r8)   dlgsqt2  !  1/ln( sqrt( 2 ) )
2177
2178
2179       real(r8)    esat01         ! aitken mode exp( log^2( sigmag )/8 )
2180       real(r8)    esac01         ! accumulation mode exp( log^2( sigmag )/8 )
2181
2182       real(r8)    esat04
2183       real(r8)    esac04
2184
2185       real(r8)    esat05
2186       real(r8)    esac05
2187
2188       real(r8)    esat08
2189       real(r8)    esac08
2190
2191       real(r8)    esat09
2192       real(r8)    esac09
2193
2194       real(r8)    esat16
2195       real(r8)    esac16
2196
2197       real(r8)    esat20
2198       real(r8)    esac20
2199
2200       real(r8)    esat24
2201       real(r8)    esac24
2202
2203       real(r8)    esat25
2204       real(r8)    esac25
2205
2206       real(r8)    esat36
2207       real(r8)    esac36
2208
2209       real(r8)    esat49
2210
2211       real(r8)    esat64
2212       real(r8)    esac64
2213
2214       real(r8)    esat100
2215
2216       real(r8) dgat2, dgac2, dgat3, dgac3
2217       real(r8) sqdgat, sqdgac
2218       real(r8) sqdgat5, sqdgac5
2219       real(r8) sqdgat7
2220       real(r8) r, r2, r3, rx4, r5, r6, rx8
2221       real(r8) ri1, ri2, ri3, ri4
2222       real(r8) rat
2223       real(r8) coagfm0, coagnc0
2224       real(r8) coagfm3, coagnc3
2225       real(r8) coagfm_at, coagfm_ac
2226       real(r8) coagnc_at, coagnc_ac
2227       real(r8) coagatat0
2228       real(r8) coagacac0
2229       real(r8) coagatat2
2230       real(r8) coagacac2
2231       real(r8) coagatac0, coagatac3
2232       real(r8) coagatac2
2233       real(r8) coagacat2
2234       real(r8) xm2at, xm3at, xm2ac, xm3ac
2235
2236 ! *** correction factors for coagulation rates
2237       real(r4), save :: bm0( 10 )          ! m0 intramodal fm - rpm values
2238       real(r4), save :: bm0ij( 10, 10, 10 ) ! m0 intermodal fm
2239       real(r4), save :: bm3i( 10, 10, 10 ) ! m3 intermodal fm- rpm values
2240       real(r4), save :: bm2ii(10) ! m2 intramodal fm
2241       real(r4), save :: bm2iitt(10) ! m2 intramodal total
2242       real(r4), save :: bm2ij(10,10,10) ! m2 intermodal fm i to j
2243       real(r4), save :: bm2ji(10,10,10) ! m2 total intermodal  j from i
2244
2245 ! *** populate the arrays for the correction factors.
2246
2247 ! rpm 0th moment correction factors for unimodal fm coagulation  rates
2248       data      bm0  /   &
2249             0.707106785165097, 0.726148960080488, 0.766430744110958,   &
2250             0.814106389441342, 0.861679526483207, 0.903600509090092,   &
2251             0.936578814219156, 0.960098926735545, 0.975646823342881,   &
2252             0.985397173215326   /
2253
2254
2255 ! fsb new fm correction factors for m0 intermodal coagulation
2256
2257       data (bm0ij (  1,  1,ibeta), ibeta = 1,10) /   &
2258         0.628539,  0.639610,  0.664514,  0.696278,  0.731558,   &
2259         0.768211,  0.804480,  0.838830,  0.870024,  0.897248/
2260       data (bm0ij (  1,  2,ibeta), ibeta = 1,10) /   &
2261         0.639178,  0.649966,  0.674432,  0.705794,  0.740642,   &
2262         0.776751,  0.812323,  0.845827,  0.876076,  0.902324/
2263       data (bm0ij (  1,  3,ibeta), ibeta = 1,10) /   &
2264         0.663109,  0.673464,  0.697147,  0.727637,  0.761425,   &
2265         0.796155,  0.829978,  0.861419,  0.889424,  0.913417/
2266       data (bm0ij (  1,  4,ibeta), ibeta = 1,10) /   &
2267         0.693693,  0.703654,  0.726478,  0.755786,  0.787980,   &
2268         0.820626,  0.851898,  0.880459,  0.905465,  0.926552/
2269       data (bm0ij (  1,  5,ibeta), ibeta = 1,10) /   &
2270         0.727803,  0.737349,  0.759140,  0.786870,  0.816901,   &
2271         0.846813,  0.874906,  0.900060,  0.921679,  0.939614/
2272       data (bm0ij (  1,  6,ibeta), ibeta = 1,10) /   &
2273         0.763461,  0.772483,  0.792930,  0.818599,  0.845905,   &
2274         0.872550,  0.897051,  0.918552,  0.936701,  0.951528/
2275       data (bm0ij (  1,  7,ibeta), ibeta = 1,10) /   &
2276         0.799021,  0.807365,  0.826094,  0.849230,  0.873358,   &
2277         0.896406,  0.917161,  0.935031,  0.949868,  0.961828/
2278       data (bm0ij (  1,  8,ibeta), ibeta = 1,10) /   &
2279         0.833004,  0.840514,  0.857192,  0.877446,  0.898147,   &
2280         0.917518,  0.934627,  0.949106,  0.960958,  0.970403/
2281       data (bm0ij (  1,  9,ibeta), ibeta = 1,10) /   &
2282         0.864172,  0.870734,  0.885153,  0.902373,  0.919640,   &
2283         0.935494,  0.949257,  0.960733,  0.970016,  0.977346/
2284       data (bm0ij (  1, 10,ibeta), ibeta = 1,10) /   &
2285         0.891658,  0.897227,  0.909343,  0.923588,  0.937629,   &
2286         0.950307,  0.961151,  0.970082,  0.977236,  0.982844/
2287       data (bm0ij (  2,  1,ibeta), ibeta = 1,10) /   &
2288         0.658724,  0.670587,  0.697539,  0.731890,  0.769467,   &
2289         0.807391,  0.843410,  0.875847,  0.903700,  0.926645/
2290       data (bm0ij (  2,  2,ibeta), ibeta = 1,10) /   &
2291         0.667070,  0.678820,  0.705538,  0.739591,  0.776758,   &
2292         0.814118,  0.849415,  0.881020,  0.908006,  0.930121/
2293       data (bm0ij (  2,  3,ibeta), ibeta = 1,10) /   &
2294         0.686356,  0.697839,  0.723997,  0.757285,  0.793389,   &
2295         0.829313,  0.862835,  0.892459,  0.917432,  0.937663/
2296       data (bm0ij (  2,  4,ibeta), ibeta = 1,10) /   &
2297         0.711425,  0.722572,  0.747941,  0.780055,  0.814518,   &
2298         0.848315,  0.879335,  0.906290,  0.928658,  0.946526/
2299       data (bm0ij (  2,  5,ibeta), ibeta = 1,10) /   &
2300         0.739575,  0.750307,  0.774633,  0.805138,  0.837408,   &
2301         0.868504,  0.896517,  0.920421,  0.939932,  0.955299/
2302       data (bm0ij (  2,  6,ibeta), ibeta = 1,10) /   &
2303         0.769143,  0.779346,  0.802314,  0.830752,  0.860333,   &
2304         0.888300,  0.913014,  0.933727,  0.950370,  0.963306/
2305       data (bm0ij (  2,  7,ibeta), ibeta = 1,10) /   &
2306         0.798900,  0.808431,  0.829700,  0.855653,  0.882163,   &
2307         0.906749,  0.928075,  0.945654,  0.959579,  0.970280/
2308       data (bm0ij (  2,  8,ibeta), ibeta = 1,10) /   &
2309         0.827826,  0.836542,  0.855808,  0.878954,  0.902174,   &
2310         0.923316,  0.941345,  0.955989,  0.967450,  0.976174/
2311       data (bm0ij (  2,  9,ibeta), ibeta = 1,10) /   &
2312         0.855068,  0.862856,  0.879900,  0.900068,  0.919956,   &
2313         0.937764,  0.952725,  0.964726,  0.974027,  0.981053/
2314       data (bm0ij (  2, 10,ibeta), ibeta = 1,10) /   &
2315         0.879961,  0.886755,  0.901484,  0.918665,  0.935346,   &
2316         0.950065,  0.962277,  0.971974,  0.979432,  0.985033/
2317       data (bm0ij (  3,  1,ibeta), ibeta = 1,10) /   &
2318         0.724166,  0.735474,  0.761359,  0.794045,  0.828702,   &
2319         0.862061,  0.891995,  0.917385,  0.937959,  0.954036/
2320       data (bm0ij (  3,  2,ibeta), ibeta = 1,10) /   &
2321         0.730416,  0.741780,  0.767647,  0.800116,  0.834344,   &
2322         0.867093,  0.896302,  0.920934,  0.940790,  0.956237/
2323       data (bm0ij (  3,  3,ibeta), ibeta = 1,10) /   &
2324         0.745327,  0.756664,  0.782255,  0.814026,  0.847107,   &
2325         0.878339,  0.905820,  0.928699,  0.946931,  0.960977/
2326       data (bm0ij (  3,  4,ibeta), ibeta = 1,10) /   &
2327         0.765195,  0.776312,  0.801216,  0.831758,  0.863079,   &
2328         0.892159,  0.917319,  0.937939,  0.954145,  0.966486/
2329       data (bm0ij (  3,  5,ibeta), ibeta = 1,10) /   &
2330         0.787632,  0.798347,  0.822165,  0.850985,  0.880049,   &
2331         0.906544,  0.929062,  0.947218,  0.961288,  0.971878/
2332       data (bm0ij (  3,  6,ibeta), ibeta = 1,10) /   &
2333         0.811024,  0.821179,  0.843557,  0.870247,  0.896694,   &
2334         0.920365,  0.940131,  0.955821,  0.967820,  0.976753/
2335       data (bm0ij (  3,  7,ibeta), ibeta = 1,10) /   &
2336         0.834254,  0.843709,  0.864356,  0.888619,  0.912245,   &
2337         0.933019,  0.950084,  0.963438,  0.973530,  0.980973/
2338       data (bm0ij (  3,  8,ibeta), ibeta = 1,10) /   &
2339         0.856531,  0.865176,  0.883881,  0.905544,  0.926290,   &
2340         0.944236,  0.958762,  0.969988,  0.978386,  0.984530/
2341       data (bm0ij (  3,  9,ibeta), ibeta = 1,10) /   &
2342         0.877307,  0.885070,  0.901716,  0.920729,  0.938663,   &
2343         0.953951,  0.966169,  0.975512,  0.982442,  0.987477/
2344       data (bm0ij (  3, 10,ibeta), ibeta = 1,10) /   &
2345         0.896234,  0.903082,  0.917645,  0.934069,  0.949354,   &
2346         0.962222,  0.972396,  0.980107,  0.985788,  0.989894/
2347       data (bm0ij (  4,  1,ibeta), ibeta = 1,10) /   &
2348         0.799294,  0.809144,  0.831293,  0.858395,  0.885897,   &
2349         0.911031,  0.932406,  0.949642,  0.963001,  0.973062/
2350       data (bm0ij (  4,  2,ibeta), ibeta = 1,10) /   &
2351         0.804239,  0.814102,  0.836169,  0.862984,  0.890003,   &
2352         0.914535,  0.935274,  0.951910,  0.964748,  0.974381/
2353       data (bm0ij (  4,  3,ibeta), ibeta = 1,10) /   &
2354         0.815910,  0.825708,  0.847403,  0.873389,  0.899185,   &
2355         0.922275,  0.941543,  0.956826,  0.968507,  0.977204/
2356       data (bm0ij (  4,  4,ibeta), ibeta = 1,10) /   &
2357         0.831348,  0.840892,  0.861793,  0.886428,  0.910463,   &
2358         0.931614,  0.948993,  0.962593,  0.972872,  0.980456/
2359       data (bm0ij (  4,  5,ibeta), ibeta = 1,10) /   &
2360         0.848597,  0.857693,  0.877402,  0.900265,  0.922180,   &
2361         0.941134,  0.956464,  0.968298,  0.977143,  0.983611/
2362       data (bm0ij (  4,  6,ibeta), ibeta = 1,10) /   &
2363         0.866271,  0.874764,  0.892984,  0.913796,  0.933407,   &
2364         0.950088,  0.963380,  0.973512,  0.981006,  0.986440/
2365       data (bm0ij (  4,  7,ibeta), ibeta = 1,10) /   &
2366         0.883430,  0.891216,  0.907762,  0.926388,  0.943660,   &
2367         0.958127,  0.969499,  0.978070,  0.984351,  0.988872/
2368       data (bm0ij (  4,  8,ibeta), ibeta = 1,10) /   &
2369         0.899483,  0.906505,  0.921294,  0.937719,  0.952729,   &
2370         0.965131,  0.974762,  0.981950,  0.987175,  0.990912/
2371       data (bm0ij (  4,  9,ibeta), ibeta = 1,10) /   &
2372         0.914096,  0.920337,  0.933373,  0.947677,  0.960579,   &
2373         0.971111,  0.979206,  0.985196,  0.989520,  0.992597/
2374       data (bm0ij (  4, 10,ibeta), ibeta = 1,10) /   &
2375         0.927122,  0.932597,  0.943952,  0.956277,  0.967268,   &
2376         0.976147,  0.982912,  0.987882,  0.991450,  0.993976/
2377       data (bm0ij (  5,  1,ibeta), ibeta = 1,10) /   &
2378         0.865049,  0.872851,  0.889900,  0.909907,  0.929290,   &
2379         0.946205,  0.959991,  0.970706,  0.978764,  0.984692/
2380       data (bm0ij (  5,  2,ibeta), ibeta = 1,10) /   &
2381         0.868989,  0.876713,  0.893538,  0.913173,  0.932080,   &
2382         0.948484,  0.961785,  0.972080,  0.979796,  0.985457/
2383       data (bm0ij (  5,  3,ibeta), ibeta = 1,10) /   &
2384         0.878010,  0.885524,  0.901756,  0.920464,  0.938235,   &
2385         0.953461,  0.965672,  0.975037,  0.982005,  0.987085/
2386       data (bm0ij (  5,  4,ibeta), ibeta = 1,10) /   &
2387         0.889534,  0.896698,  0.912012,  0.929395,  0.945647,   &
2388         0.959366,  0.970227,  0.978469,  0.984547,  0.988950/
2389       data (bm0ij (  5,  5,ibeta), ibeta = 1,10) /   &
2390         0.902033,  0.908713,  0.922848,  0.938648,  0.953186,   &
2391         0.965278,  0.974729,  0.981824,  0.987013,  0.990746/
2392       data (bm0ij (  5,  6,ibeta), ibeta = 1,10) /   &
2393         0.914496,  0.920599,  0.933389,  0.947485,  0.960262,   &
2394         0.970743,  0.978839,  0.984858,  0.989225,  0.992348/
2395       data (bm0ij (  5,  7,ibeta), ibeta = 1,10) /   &
2396         0.926281,  0.931761,  0.943142,  0.955526,  0.966600,   &
2397         0.975573,  0.982431,  0.987485,  0.991128,  0.993718/
2398       data (bm0ij (  5,  8,ibeta), ibeta = 1,10) /   &
2399         0.937029,  0.941877,  0.951868,  0.962615,  0.972112,   &
2400         0.979723,  0.985488,  0.989705,  0.992725,  0.994863/
2401       data (bm0ij (  5,  9,ibeta), ibeta = 1,10) /   &
2402         0.946580,  0.950819,  0.959494,  0.968732,  0.976811,   &
2403         0.983226,  0.988047,  0.991550,  0.994047,  0.995806/
2404       data (bm0ij (  5, 10,ibeta), ibeta = 1,10) /   &
2405         0.954909,  0.958581,  0.966049,  0.973933,  0.980766,   &
2406         0.986149,  0.990166,  0.993070,  0.995130,  0.996577/
2407       data (bm0ij (  6,  1,ibeta), ibeta = 1,10) /   &
2408         0.914182,  0.919824,  0.931832,  0.945387,  0.957999,   &
2409         0.968606,  0.976982,  0.983331,  0.988013,  0.991407/
2410       data (bm0ij (  6,  2,ibeta), ibeta = 1,10) /   &
2411         0.917139,  0.922665,  0.934395,  0.947580,  0.959792,   &
2412         0.970017,  0.978062,  0.984138,  0.988609,  0.991843/
2413       data (bm0ij (  6,  3,ibeta), ibeta = 1,10) /   &
2414         0.923742,  0.928990,  0.940064,  0.952396,  0.963699,   &
2415         0.973070,  0.980381,  0.985866,  0.989878,  0.992768/
2416       data (bm0ij (  6,  4,ibeta), ibeta = 1,10) /   &
2417         0.931870,  0.936743,  0.946941,  0.958162,  0.968318,   &
2418         0.976640,  0.983069,  0.987853,  0.991330,  0.993822/
2419       data (bm0ij (  6,  5,ibeta), ibeta = 1,10) /   &
2420         0.940376,  0.944807,  0.954004,  0.963999,  0.972928,   &
2421         0.980162,  0.985695,  0.989779,  0.992729,  0.994833/
2422       data (bm0ij (  6,  6,ibeta), ibeta = 1,10) /   &
2423         0.948597,  0.952555,  0.960703,  0.969454,  0.977181,   &
2424         0.983373,  0.988067,  0.991507,  0.993977,  0.995730/
2425       data (bm0ij (  6,  7,ibeta), ibeta = 1,10) /   &
2426         0.956167,  0.959648,  0.966763,  0.974326,  0.980933,   &
2427         0.986177,  0.990121,  0.992993,  0.995045,  0.996495/
2428       data (bm0ij (  6,  8,ibeta), ibeta = 1,10) /   &
2429         0.962913,  0.965937,  0.972080,  0.978552,  0.984153,   &
2430         0.988563,  0.991857,  0.994242,  0.995938,  0.997133/
2431       data (bm0ij (  6,  9,ibeta), ibeta = 1,10) /   &
2432         0.968787,  0.971391,  0.976651,  0.982148,  0.986869,   &
2433         0.990560,  0.993301,  0.995275,  0.996675,  0.997657/
2434       data (bm0ij (  6, 10,ibeta), ibeta = 1,10) /   &
2435         0.973822,  0.976047,  0.980523,  0.985170,  0.989134,   &
2436         0.992215,  0.994491,  0.996124,  0.997277,  0.998085/
2437       data (bm0ij (  7,  1,ibeta), ibeta = 1,10) /   &
2438         0.947410,  0.951207,  0.959119,  0.967781,  0.975592,   &
2439         0.981981,  0.986915,  0.990590,  0.993266,  0.995187/
2440       data (bm0ij (  7,  2,ibeta), ibeta = 1,10) /   &
2441         0.949477,  0.953161,  0.960824,  0.969187,  0.976702,   &
2442         0.982831,  0.987550,  0.991057,  0.993606,  0.995434/
2443       data (bm0ij (  7,  3,ibeta), ibeta = 1,10) /   &
2444         0.954008,  0.957438,  0.964537,  0.972232,  0.979095,   &
2445         0.984653,  0.988907,  0.992053,  0.994330,  0.995958/
2446       data (bm0ij (  7,  4,ibeta), ibeta = 1,10) /   &
2447         0.959431,  0.962539,  0.968935,  0.975808,  0.981882,   &
2448         0.986759,  0.990466,  0.993190,  0.995153,  0.996552/
2449       data (bm0ij (  7,  5,ibeta), ibeta = 1,10) /   &
2450         0.964932,  0.967693,  0.973342,  0.979355,  0.984620,   &
2451         0.988812,  0.991974,  0.994285,  0.995943,  0.997119/
2452       data (bm0ij (  7,  6,ibeta), ibeta = 1,10) /   &
2453         0.970101,  0.972517,  0.977428,  0.982612,  0.987110,   &
2454         0.990663,  0.993326,  0.995261,  0.996644,  0.997621/
2455       data (bm0ij (  7,  7,ibeta), ibeta = 1,10) /   &
2456         0.974746,  0.976834,  0.981055,  0.985475,  0.989280,   &
2457         0.992265,  0.994488,  0.996097,  0.997241,  0.998048/
2458       data (bm0ij (  7,  8,ibeta), ibeta = 1,10) /   &
2459         0.978804,  0.980591,  0.984187,  0.987927,  0.991124,   &
2460         0.993617,  0.995464,  0.996795,  0.997739,  0.998403/
2461       data (bm0ij (  7,  9,ibeta), ibeta = 1,10) /   &
2462         0.982280,  0.983799,  0.986844,  0.989991,  0.992667,   &
2463         0.994742,  0.996273,  0.997372,  0.998149,  0.998695/
2464       data (bm0ij (  7, 10,ibeta), ibeta = 1,10) /   &
2465         0.985218,  0.986503,  0.989071,  0.991711,  0.993945,   &
2466         0.995669,  0.996937,  0.997844,  0.998484,  0.998932/
2467       data (bm0ij (  8,  1,ibeta), ibeta = 1,10) /   &
2468         0.968507,  0.970935,  0.975916,  0.981248,  0.985947,   &
2469         0.989716,  0.992580,  0.994689,  0.996210,  0.997297/
2470       data (bm0ij (  8,  2,ibeta), ibeta = 1,10) /   &
2471         0.969870,  0.972210,  0.977002,  0.982119,  0.986619,   &
2472         0.990219,  0.992951,  0.994958,  0.996405,  0.997437/
2473       data (bm0ij (  8,  3,ibeta), ibeta = 1,10) /   &
2474         0.972820,  0.974963,  0.979339,  0.983988,  0.988054,   &
2475         0.991292,  0.993738,  0.995529,  0.996817,  0.997734/
2476       data (bm0ij (  8,  4,ibeta), ibeta = 1,10) /   &
2477         0.976280,  0.978186,  0.982060,  0.986151,  0.989706,   &
2478         0.992520,  0.994636,  0.996179,  0.997284,  0.998069/
2479       data (bm0ij (  8,  5,ibeta), ibeta = 1,10) /   &
2480         0.979711,  0.981372,  0.984735,  0.988263,  0.991309,   &
2481         0.993706,  0.995499,  0.996801,  0.997730,  0.998389/
2482       data (bm0ij (  8,  6,ibeta), ibeta = 1,10) /   &
2483         0.982863,  0.984292,  0.987172,  0.990174,  0.992750,   &
2484         0.994766,  0.996266,  0.997352,  0.998125,  0.998670/
2485       data (bm0ij (  8,  7,ibeta), ibeta = 1,10) /   &
2486         0.985642,  0.986858,  0.989301,  0.991834,  0.993994,   &
2487         0.995676,  0.996923,  0.997822,  0.998460,  0.998910/
2488       data (bm0ij (  8,  8,ibeta), ibeta = 1,10) /   &
2489         0.988029,  0.989058,  0.991116,  0.993240,  0.995043,   &
2490         0.996440,  0.997472,  0.998214,  0.998739,  0.999108/
2491       data (bm0ij (  8,  9,ibeta), ibeta = 1,10) /   &
2492         0.990046,  0.990912,  0.992640,  0.994415,  0.995914,   &
2493         0.997073,  0.997925,  0.998536,  0.998968,  0.999271/
2494       data (bm0ij (  8, 10,ibeta), ibeta = 1,10) /   &
2495         0.991732,  0.992459,  0.993906,  0.995386,  0.996633,   &
2496         0.997592,  0.998296,  0.998799,  0.999154,  0.999403/
2497       data (bm0ij (  9,  1,ibeta), ibeta = 1,10) /   &
2498         0.981392,  0.982893,  0.985938,  0.989146,  0.991928,   &
2499         0.994129,  0.995783,  0.996991,  0.997857,  0.998473/
2500       data (bm0ij (  9,  2,ibeta), ibeta = 1,10) /   &
2501         0.982254,  0.983693,  0.986608,  0.989673,  0.992328,   &
2502         0.994424,  0.995998,  0.997146,  0.997969,  0.998553/
2503       data (bm0ij (  9,  3,ibeta), ibeta = 1,10) /   &
2504         0.984104,  0.985407,  0.988040,  0.990798,  0.993178,   &
2505         0.995052,  0.996454,  0.997474,  0.998204,  0.998722/
2506       data (bm0ij (  9,  4,ibeta), ibeta = 1,10) /   &
2507         0.986243,  0.987386,  0.989687,  0.992087,  0.994149,   &
2508         0.995765,  0.996971,  0.997846,  0.998470,  0.998913/
2509       data (bm0ij (  9,  5,ibeta), ibeta = 1,10) /   &
2510         0.988332,  0.989313,  0.991284,  0.993332,  0.995082,   &
2511         0.996449,  0.997465,  0.998200,  0.998723,  0.999093/
2512       data (bm0ij (  9,  6,ibeta), ibeta = 1,10) /   &
2513         0.990220,  0.991053,  0.992721,  0.994445,  0.995914,   &
2514         0.997056,  0.997902,  0.998513,  0.998947,  0.999253/
2515       data (bm0ij (  9,  7,ibeta), ibeta = 1,10) /   &
2516         0.991859,  0.992561,  0.993961,  0.995403,  0.996626,   &
2517         0.997574,  0.998274,  0.998778,  0.999136,  0.999387/
2518       data (bm0ij (  9,  8,ibeta), ibeta = 1,10) /   &
2519         0.993250,  0.993837,  0.995007,  0.996208,  0.997223,   &
2520         0.998007,  0.998584,  0.998999,  0.999293,  0.999499/
2521       data (bm0ij (  9,  9,ibeta), ibeta = 1,10) /   &
2522         0.994413,  0.994903,  0.995878,  0.996876,  0.997716,   &
2523         0.998363,  0.998839,  0.999180,  0.999421,  0.999591/
2524       data (bm0ij (  9, 10,ibeta), ibeta = 1,10) /   &
2525         0.995376,  0.995785,  0.996597,  0.997425,  0.998121,   &
2526         0.998655,  0.999048,  0.999328,  0.999526,  0.999665/
2527       data (bm0ij ( 10,  1,ibeta), ibeta = 1,10) /   &
2528         0.989082,  0.989991,  0.991819,  0.993723,  0.995357,   &
2529         0.996637,  0.997592,  0.998286,  0.998781,  0.999132/
2530       data (bm0ij ( 10,  2,ibeta), ibeta = 1,10) /   &
2531         0.989613,  0.990480,  0.992224,  0.994039,  0.995594,   &
2532         0.996810,  0.997717,  0.998375,  0.998845,  0.999178/
2533       data (bm0ij ( 10,  3,ibeta), ibeta = 1,10) /   &
2534         0.990744,  0.991523,  0.993086,  0.994708,  0.996094,   &
2535         0.997176,  0.997981,  0.998564,  0.998980,  0.999274/
2536       data (bm0ij ( 10,  4,ibeta), ibeta = 1,10) /   &
2537         0.992041,  0.992716,  0.994070,  0.995470,  0.996662,   &
2538         0.997591,  0.998280,  0.998778,  0.999133,  0.999383/
2539       data (bm0ij ( 10,  5,ibeta), ibeta = 1,10) /   &
2540         0.993292,  0.993867,  0.995015,  0.996199,  0.997205,   &
2541         0.997985,  0.998564,  0.998981,  0.999277,  0.999487/
2542       data (bm0ij ( 10,  6,ibeta), ibeta = 1,10) /   &
2543         0.994411,  0.994894,  0.995857,  0.996847,  0.997685,   &
2544         0.998334,  0.998814,  0.999159,  0.999404,  0.999577/
2545       data (bm0ij ( 10,  7,ibeta), ibeta = 1,10) /   &
2546         0.995373,  0.995776,  0.996577,  0.997400,  0.998094,   &
2547         0.998630,  0.999026,  0.999310,  0.999512,  0.999654/
2548       data (bm0ij ( 10,  8,ibeta), ibeta = 1,10) /   &
2549         0.996181,  0.996516,  0.997181,  0.997861,  0.998435,   &
2550         0.998877,  0.999202,  0.999435,  0.999601,  0.999717/
2551       data (bm0ij ( 10,  9,ibeta), ibeta = 1,10) /   &
2552         0.996851,  0.997128,  0.997680,  0.998242,  0.998715,   &
2553         0.999079,  0.999346,  0.999538,  0.999673,  0.999769/
2554       data (bm0ij ( 10, 10,ibeta), ibeta = 1,10) /   &
2555         0.997402,  0.997632,  0.998089,  0.998554,  0.998945,   &
2556         0.999244,  0.999464,  0.999622,  0.999733,  0.999811/
2557
2558
2559 ! rpm....   3rd moment nuclei mode corr. fac. for bimodal fm coag rate
2560
2561        data (bm3i( 1, 1,ibeta ), ibeta=1,10)/   &
2562        0.70708,0.71681,0.73821,0.76477,0.79350,0.82265,0.85090,0.87717,   &
2563        0.90069,0.92097/
2564        data (bm3i( 1, 2,ibeta ), ibeta=1,10)/   &
2565        0.72172,0.73022,0.74927,0.77324,0.79936,0.82601,0.85199,0.87637,   &
2566        0.89843,0.91774/
2567        data (bm3i( 1, 3,ibeta ), ibeta=1,10)/   &
2568        0.78291,0.78896,0.80286,0.82070,0.84022,0.85997,0.87901,0.89669,   &
2569        0.91258,0.92647/
2570        data (bm3i( 1, 4,ibeta ), ibeta=1,10)/   &
2571        0.87760,0.88147,0.89025,0.90127,0.91291,0.92420,0.93452,0.94355,   &
2572        0.95113,0.95726/
2573        data (bm3i( 1, 5,ibeta ), ibeta=1,10)/   &
2574        0.94988,0.95184,0.95612,0.96122,0.96628,0.97085,0.97467,0.97763,   &
2575        0.97971,0.98089/
2576        data (bm3i( 1, 6,ibeta ), ibeta=1,10)/   &
2577        0.98318,0.98393,0.98551,0.98728,0.98889,0.99014,0.99095,0.99124,   &
2578        0.99100,0.99020/
2579        data (bm3i( 1, 7,ibeta ), ibeta=1,10)/   &
2580        0.99480,0.99504,0.99551,0.99598,0.99629,0.99635,0.99611,0.99550,   &
2581        0.99450,0.99306/
2582        data (bm3i( 1, 8,ibeta ), ibeta=1,10)/   &
2583        0.99842,0.99848,0.99858,0.99861,0.99850,0.99819,0.99762,0.99674,   &
2584        0.99550,0.99388/
2585        data (bm3i( 1, 9,ibeta ), ibeta=1,10)/   &
2586        0.99951,0.99951,0.99949,0.99939,0.99915,0.99872,0.99805,0.99709,   &
2587        0.99579,0.99411/
2588        data (bm3i( 1,10,ibeta ), ibeta=1,10)/   &
2589        0.99984,0.99982,0.99976,0.99962,0.99934,0.99888,0.99818,0.99719,   &
2590        0.99587,0.99417/
2591        data (bm3i( 2, 1,ibeta ), ibeta=1,10)/   &
2592        0.72957,0.73993,0.76303,0.79178,0.82245,0.85270,0.88085,0.90578,   &
2593        0.92691,0.94415/
2594        data (bm3i( 2, 2,ibeta ), ibeta=1,10)/   &
2595        0.72319,0.73320,0.75547,0.78323,0.81307,0.84287,0.87107,0.89651,   &
2596        0.91852,0.93683/
2597        data (bm3i( 2, 3,ibeta ), ibeta=1,10)/   &
2598        0.74413,0.75205,0.76998,0.79269,0.81746,0.84258,0.86685,0.88938,   &
2599        0.90953,0.92695/
2600        data (bm3i( 2, 4,ibeta ), ibeta=1,10)/   &
2601        0.82588,0.83113,0.84309,0.85825,0.87456,0.89072,0.90594,0.91972,   &
2602        0.93178,0.94203/
2603        data (bm3i( 2, 5,ibeta ), ibeta=1,10)/   &
2604        0.91886,0.92179,0.92831,0.93624,0.94434,0.95192,0.95856,0.96409,   &
2605        0.96845,0.97164/
2606        data (bm3i( 2, 6,ibeta ), ibeta=1,10)/   &
2607        0.97129,0.97252,0.97515,0.97818,0.98108,0.98354,0.98542,0.98665,   &
2608        0.98721,0.98709/
2609        data (bm3i( 2, 7,ibeta ), ibeta=1,10)/   &
2610        0.99104,0.99145,0.99230,0.99320,0.99394,0.99439,0.99448,0.99416,   &
2611        0.99340,0.99217/
2612        data (bm3i( 2, 8,ibeta ), ibeta=1,10)/   &
2613        0.99730,0.99741,0.99763,0.99779,0.99782,0.99762,0.99715,0.99636,   &
2614        0.99519,0.99363/
2615        data (bm3i( 2, 9,ibeta ), ibeta=1,10)/   &
2616        0.99917,0.99919,0.99921,0.99915,0.99895,0.99856,0.99792,0.99698,   &
2617        0.99570,0.99404/
2618        data (bm3i( 2,10,ibeta ), ibeta=1,10)/   &
2619        0.99973,0.99973,0.99968,0.99955,0.99928,0.99883,0.99814,0.99716,   &
2620        0.99584,0.99415/
2621        data (bm3i( 3, 1,ibeta ), ibeta=1,10)/   &
2622        0.78358,0.79304,0.81445,0.84105,0.86873,0.89491,0.91805,0.93743,   &
2623        0.95300,0.96510/
2624        data (bm3i( 3, 2,ibeta ), ibeta=1,10)/   &
2625        0.76412,0.77404,0.79635,0.82404,0.85312,0.88101,0.90610,0.92751,   &
2626        0.94500,0.95879/
2627        data (bm3i( 3, 3,ibeta ), ibeta=1,10)/   &
2628        0.74239,0.75182,0.77301,0.79956,0.82809,0.85639,0.88291,0.90658,   &
2629        0.92683,0.94350/
2630        data (bm3i( 3, 4,ibeta ), ibeta=1,10)/   &
2631        0.78072,0.78758,0.80317,0.82293,0.84437,0.86589,0.88643,0.90526,   &
2632        0.92194,0.93625/
2633        data (bm3i( 3, 5,ibeta ), ibeta=1,10)/   &
2634        0.87627,0.88044,0.88981,0.90142,0.91357,0.92524,0.93585,0.94510,   &
2635        0.95285,0.95911/
2636        data (bm3i( 3, 6,ibeta ), ibeta=1,10)/   &
2637        0.95176,0.95371,0.95796,0.96297,0.96792,0.97233,0.97599,0.97880,   &
2638        0.98072,0.98178/
2639        data (bm3i( 3, 7,ibeta ), ibeta=1,10)/   &
2640        0.98453,0.98523,0.98670,0.98833,0.98980,0.99092,0.99160,0.99179,   &
2641        0.99145,0.99058/
2642        data (bm3i( 3, 8,ibeta ), ibeta=1,10)/   &
2643        0.99534,0.99555,0.99597,0.99637,0.99662,0.99663,0.99633,0.99569,   &
2644        0.99465,0.99318/
2645        data (bm3i( 3, 9,ibeta ), ibeta=1,10)/   &
2646        0.99859,0.99864,0.99872,0.99873,0.99860,0.99827,0.99768,0.99679,   &
2647        0.99555,0.99391/
2648        data (bm3i( 3,10,ibeta ), ibeta=1,10)/   &
2649        0.99956,0.99956,0.99953,0.99942,0.99918,0.99875,0.99807,0.99711,   &
2650        0.99580,0.99412/
2651        data (bm3i( 4, 1,ibeta ), ibeta=1,10)/   &
2652        0.84432,0.85223,0.86990,0.89131,0.91280,0.93223,0.94861,0.96172,   &
2653        0.97185,0.97945/
2654        data (bm3i( 4, 2,ibeta ), ibeta=1,10)/   &
2655        0.82299,0.83164,0.85101,0.87463,0.89857,0.92050,0.93923,0.95443,   &
2656        0.96629,0.97529/
2657        data (bm3i( 4, 3,ibeta ), ibeta=1,10)/   &
2658        0.77870,0.78840,0.81011,0.83690,0.86477,0.89124,0.91476,0.93460,   &
2659        0.95063,0.96316/
2660        data (bm3i( 4, 4,ibeta ), ibeta=1,10)/   &
2661        0.76386,0.77233,0.79147,0.81557,0.84149,0.86719,0.89126,0.91275,   &
2662        0.93116,0.94637/
2663        data (bm3i( 4, 5,ibeta ), ibeta=1,10)/   &
2664        0.82927,0.83488,0.84756,0.86346,0.88040,0.89704,0.91257,0.92649,   &
2665        0.93857,0.94874/
2666        data (bm3i( 4, 6,ibeta ), ibeta=1,10)/   &
2667        0.92184,0.92481,0.93136,0.93925,0.94724,0.95462,0.96104,0.96634,   &
2668        0.97048,0.97348/
2669        data (bm3i( 4, 7,ibeta ), ibeta=1,10)/   &
2670        0.97341,0.97457,0.97706,0.97991,0.98260,0.98485,0.98654,0.98760,   &
2671        0.98801,0.98777/
2672        data (bm3i( 4, 8,ibeta ), ibeta=1,10)/   &
2673        0.99192,0.99229,0.99305,0.99385,0.99449,0.99486,0.99487,0.99449,   &
2674        0.99367,0.99239/
2675        data (bm3i( 4, 9,ibeta ), ibeta=1,10)/   &
2676        0.99758,0.99768,0.99787,0.99800,0.99799,0.99777,0.99727,0.99645,   &
2677        0.99527,0.99369/
2678        data (bm3i( 4,10,ibeta ), ibeta=1,10)/   &
2679        0.99926,0.99928,0.99928,0.99921,0.99900,0.99860,0.99795,0.99701,   &
2680        0.99572,0.99405/
2681        data (bm3i( 5, 1,ibeta ), ibeta=1,10)/   &
2682        0.89577,0.90190,0.91522,0.93076,0.94575,0.95876,0.96932,0.97751,   &
2683        0.98367,0.98820/
2684        data (bm3i( 5, 2,ibeta ), ibeta=1,10)/   &
2685        0.87860,0.88547,0.90052,0.91828,0.93557,0.95075,0.96319,0.97292,   &
2686        0.98028,0.98572/
2687        data (bm3i( 5, 3,ibeta ), ibeta=1,10)/   &
2688        0.83381,0.84240,0.86141,0.88425,0.90707,0.92770,0.94510,0.95906,   &
2689        0.96986,0.97798/
2690        data (bm3i( 5, 4,ibeta ), ibeta=1,10)/   &
2691        0.78530,0.79463,0.81550,0.84127,0.86813,0.89367,0.91642,0.93566,   &
2692        0.95125,0.96347/
2693        data (bm3i( 5, 5,ibeta ), ibeta=1,10)/   &
2694        0.79614,0.80332,0.81957,0.84001,0.86190,0.88351,0.90368,0.92169,   &
2695        0.93718,0.95006/
2696        data (bm3i( 5, 6,ibeta ), ibeta=1,10)/   &
2697        0.88192,0.88617,0.89565,0.90728,0.91931,0.93076,0.94107,0.94997,   &
2698        0.95739,0.96333/
2699        data (bm3i( 5, 7,ibeta ), ibeta=1,10)/   &
2700        0.95509,0.95698,0.96105,0.96583,0.97048,0.97460,0.97796,0.98050,   &
2701        0.98218,0.98304/
2702        data (bm3i( 5, 8,ibeta ), ibeta=1,10)/   &
2703        0.98596,0.98660,0.98794,0.98943,0.99074,0.99172,0.99227,0.99235,   &
2704        0.99192,0.99096/
2705        data (bm3i( 5, 9,ibeta ), ibeta=1,10)/   &
2706        0.99581,0.99600,0.99637,0.99672,0.99691,0.99687,0.99653,0.99585,   &
2707        0.99478,0.99329/
2708        data (bm3i( 5,10,ibeta ), ibeta=1,10)/   &
2709        0.99873,0.99878,0.99884,0.99883,0.99869,0.99834,0.99774,0.99684,   &
2710        0.99558,0.99394/
2711        data (bm3i( 6, 1,ibeta ), ibeta=1,10)/   &
2712        0.93335,0.93777,0.94711,0.95764,0.96741,0.97562,0.98210,0.98701,   &
2713        0.99064,0.99327/
2714        data (bm3i( 6, 2,ibeta ), ibeta=1,10)/   &
2715        0.92142,0.92646,0.93723,0.94947,0.96096,0.97069,0.97842,0.98431,   &
2716        0.98868,0.99186/
2717        data (bm3i( 6, 3,ibeta ), ibeta=1,10)/   &
2718        0.88678,0.89351,0.90810,0.92508,0.94138,0.95549,0.96693,0.97578,   &
2719        0.98243,0.98731/
2720        data (bm3i( 6, 4,ibeta ), ibeta=1,10)/   &
2721        0.83249,0.84124,0.86051,0.88357,0.90655,0.92728,0.94477,0.95880,   &
2722        0.96964,0.97779/
2723        data (bm3i( 6, 5,ibeta ), ibeta=1,10)/   &
2724        0.79593,0.80444,0.82355,0.84725,0.87211,0.89593,0.91735,0.93566,   &
2725        0.95066,0.96255/
2726        data (bm3i( 6, 6,ibeta ), ibeta=1,10)/   &
2727        0.84124,0.84695,0.85980,0.87575,0.89256,0.90885,0.92383,0.93704,   &
2728        0.94830,0.95761/
2729        data (bm3i( 6, 7,ibeta ), ibeta=1,10)/   &
2730        0.92721,0.93011,0.93647,0.94406,0.95166,0.95862,0.96460,0.96949,   &
2731        0.97326,0.97595/
2732        data (bm3i( 6, 8,ibeta ), ibeta=1,10)/   &
2733        0.97573,0.97681,0.97913,0.98175,0.98421,0.98624,0.98772,0.98860,   &
2734        0.98885,0.98847/
2735        data (bm3i( 6, 9,ibeta ), ibeta=1,10)/   &
2736        0.99271,0.99304,0.99373,0.99444,0.99499,0.99528,0.99522,0.99477,   &
2737        0.99390,0.99258/
2738        data (bm3i( 6,10,ibeta ), ibeta=1,10)/   &
2739        0.99782,0.99791,0.99807,0.99817,0.99813,0.99788,0.99737,0.99653,   &
2740        0.99533,0.99374/
2741        data (bm3i( 7, 1,ibeta ), ibeta=1,10)/   &
2742        0.95858,0.96158,0.96780,0.97460,0.98073,0.98575,0.98963,0.99252,   &
2743        0.99463,0.99615/
2744        data (bm3i( 7, 2,ibeta ), ibeta=1,10)/   &
2745        0.95091,0.95438,0.96163,0.96962,0.97688,0.98286,0.98751,0.99099,   &
2746        0.99353,0.99536/
2747        data (bm3i( 7, 3,ibeta ), ibeta=1,10)/   &
2748        0.92751,0.93233,0.94255,0.95406,0.96473,0.97366,0.98070,0.98602,   &
2749        0.98994,0.99278/
2750        data (bm3i( 7, 4,ibeta ), ibeta=1,10)/   &
2751        0.88371,0.89075,0.90595,0.92351,0.94028,0.95474,0.96642,0.97544,   &
2752        0.98220,0.98715/
2753        data (bm3i( 7, 5,ibeta ), ibeta=1,10)/   &
2754        0.82880,0.83750,0.85671,0.87980,0.90297,0.92404,0.94195,0.95644,   &
2755        0.96772,0.97625/
2756        data (bm3i( 7, 6,ibeta ), ibeta=1,10)/   &
2757        0.81933,0.82655,0.84279,0.86295,0.88412,0.90449,0.92295,0.93890,   &
2758        0.95215,0.96281/
2759        data (bm3i( 7, 7,ibeta ), ibeta=1,10)/   &
2760        0.89099,0.89519,0.90448,0.91577,0.92732,0.93820,0.94789,0.95616,   &
2761        0.96297,0.96838/
2762        data (bm3i( 7, 8,ibeta ), ibeta=1,10)/   &
2763        0.95886,0.96064,0.96448,0.96894,0.97324,0.97701,0.98004,0.98228,   &
2764        0.98371,0.98435/
2765        data (bm3i( 7, 9,ibeta ), ibeta=1,10)/   &
2766        0.98727,0.98786,0.98908,0.99043,0.99160,0.99245,0.99288,0.99285,   &
2767        0.99234,0.99131/
2768        data (bm3i( 7,10,ibeta ), ibeta=1,10)/   &
2769        0.99621,0.99638,0.99671,0.99700,0.99715,0.99707,0.99670,0.99599,   &
2770        0.99489,0.99338/
2771        data (bm3i( 8, 1,ibeta ), ibeta=1,10)/   &
2772        0.97470,0.97666,0.98064,0.98491,0.98867,0.99169,0.99399,0.99569,   &
2773        0.99691,0.99779/
2774        data (bm3i( 8, 2,ibeta ), ibeta=1,10)/   &
2775        0.96996,0.97225,0.97693,0.98196,0.98643,0.99003,0.99279,0.99482,   &
2776        0.99630,0.99735/
2777        data (bm3i( 8, 3,ibeta ), ibeta=1,10)/   &
2778        0.95523,0.95848,0.96522,0.97260,0.97925,0.98468,0.98888,0.99200,   &
2779        0.99427,0.99590/
2780        data (bm3i( 8, 4,ibeta ), ibeta=1,10)/   &
2781        0.92524,0.93030,0.94098,0.95294,0.96397,0.97317,0.98038,0.98582,   &
2782        0.98981,0.99270/
2783        data (bm3i( 8, 5,ibeta ), ibeta=1,10)/   &
2784        0.87576,0.88323,0.89935,0.91799,0.93583,0.95126,0.96377,0.97345,   &
2785        0.98072,0.98606/
2786        data (bm3i( 8, 6,ibeta ), ibeta=1,10)/   &
2787        0.83078,0.83894,0.85705,0.87899,0.90126,0.92179,0.93950,0.95404,   &
2788        0.96551,0.97430/
2789        data (bm3i( 8, 7,ibeta ), ibeta=1,10)/   &
2790        0.85727,0.86294,0.87558,0.89111,0.90723,0.92260,0.93645,0.94841,   &
2791        0.95838,0.96643/
2792        data (bm3i( 8, 8,ibeta ), ibeta=1,10)/   &
2793        0.93337,0.93615,0.94220,0.94937,0.95647,0.96292,0.96840,0.97283,   &
2794        0.97619,0.97854/
2795        data (bm3i( 8, 9,ibeta ), ibeta=1,10)/   &
2796        0.97790,0.97891,0.98105,0.98346,0.98569,0.98751,0.98879,0.98950,   &
2797        0.98961,0.98912/
2798        data (bm3i( 8,10,ibeta ), ibeta=1,10)/   &
2799        0.99337,0.99367,0.99430,0.99493,0.99541,0.99562,0.99551,0.99501,   &
2800        0.99410,0.99274/
2801        data (bm3i( 9, 1,ibeta ), ibeta=1,10)/   &
2802        0.98470,0.98594,0.98844,0.99106,0.99334,0.99514,0.99650,0.99749,   &
2803        0.99821,0.99872/
2804        data (bm3i( 9, 2,ibeta ), ibeta=1,10)/   &
2805        0.98184,0.98330,0.98624,0.98934,0.99205,0.99420,0.99582,0.99701,   &
2806        0.99787,0.99848/
2807        data (bm3i( 9, 3,ibeta ), ibeta=1,10)/   &
2808        0.97288,0.97498,0.97927,0.98385,0.98789,0.99113,0.99360,0.99541,   &
2809        0.99673,0.99766/
2810        data (bm3i( 9, 4,ibeta ), ibeta=1,10)/   &
2811        0.95403,0.95741,0.96440,0.97202,0.97887,0.98444,0.98872,0.99190,   &
2812        0.99421,0.99586/
2813        data (bm3i( 9, 5,ibeta ), ibeta=1,10)/   &
2814        0.91845,0.92399,0.93567,0.94873,0.96076,0.97079,0.97865,0.98457,   &
2815        0.98892,0.99206/
2816        data (bm3i( 9, 6,ibeta ), ibeta=1,10)/   &
2817        0.86762,0.87533,0.89202,0.91148,0.93027,0.94669,0.96013,0.97062,   &
2818        0.97855,0.98441/
2819        data (bm3i( 9, 7,ibeta ), ibeta=1,10)/   &
2820        0.84550,0.85253,0.86816,0.88721,0.90671,0.92490,0.94083,0.95413,   &
2821        0.96481,0.97314/
2822        data (bm3i( 9, 8,ibeta ), ibeta=1,10)/   &
2823        0.90138,0.90544,0.91437,0.92513,0.93602,0.94615,0.95506,0.96258,   &
2824        0.96868,0.97347/
2825        data (bm3i( 9, 9,ibeta ), ibeta=1,10)/   &
2826        0.96248,0.96415,0.96773,0.97187,0.97583,0.97925,0.98198,0.98394,   &
2827        0.98514,0.98559/
2828        data (bm3i( 9,10,ibeta ), ibeta=1,10)/   &
2829        0.98837,0.98892,0.99005,0.99127,0.99232,0.99306,0.99339,0.99328,   &
2830        0.99269,0.99161/
2831        data (bm3i(10, 1,ibeta ), ibeta=1,10)/   &
2832        0.99080,0.99158,0.99311,0.99471,0.99607,0.99715,0.99795,0.99853,   &
2833        0.99895,0.99925/
2834        data (bm3i(10, 2,ibeta ), ibeta=1,10)/   &
2835        0.98910,0.99001,0.99182,0.99371,0.99533,0.99661,0.99757,0.99826,   &
2836        0.99876,0.99912/
2837        data (bm3i(10, 3,ibeta ), ibeta=1,10)/   &
2838        0.98374,0.98506,0.98772,0.99051,0.99294,0.99486,0.99630,0.99736,   &
2839        0.99812,0.99866/
2840        data (bm3i(10, 4,ibeta ), ibeta=1,10)/   &
2841        0.97238,0.97453,0.97892,0.98361,0.98773,0.99104,0.99354,0.99538,   &
2842        0.99671,0.99765/
2843        data (bm3i(10, 5,ibeta ), ibeta=1,10)/   &
2844        0.94961,0.95333,0.96103,0.96941,0.97693,0.98303,0.98772,0.99119,   &
2845        0.99371,0.99551/
2846        data (bm3i(10, 6,ibeta ), ibeta=1,10)/   &
2847        0.90943,0.91550,0.92834,0.94275,0.95608,0.96723,0.97600,0.98263,   &
2848        0.98751,0.99103/
2849        data (bm3i(10, 7,ibeta ), ibeta=1,10)/   &
2850        0.86454,0.87200,0.88829,0.90749,0.92630,0.94300,0.95687,0.96785,   &
2851        0.97626,0.98254/
2852        data (bm3i(10, 8,ibeta ), ibeta=1,10)/   &
2853        0.87498,0.88048,0.89264,0.90737,0.92240,0.93642,0.94877,0.95917,   &
2854        0.96762,0.97429/
2855        data (bm3i(10, 9,ibeta ), ibeta=1,10)/   &
2856        0.93946,0.94209,0.94781,0.95452,0.96111,0.96704,0.97203,0.97602,   &
2857        0.97900,0.98106/
2858        data (bm3i(10,10,ibeta ), ibeta=1,10)/   &
2859        0.97977,0.98071,0.98270,0.98492,0.98695,0.98858,0.98970,0.99027,   &
2860        0.99026,0.98968/
2861
2862 ! fsb fm correction for intramodal m2 coagulation
2863        data bm2ii /   &
2864         0.707107,  0.720583,  0.745310,  0.748056,  0.696935,   &
2865         0.604164,  0.504622,  0.416559,  0.343394,  0.283641/
2866
2867 ! *** total correction for intramodal m2 coagulation
2868
2869       data bm2iitt /   &
2870         1.000000,  0.907452,  0.680931,  0.409815,  0.196425,   &
2871         0.078814,  0.028473,  0.009800,  0.003322,  0.001129/
2872
2873
2874 ! fsb fm correction for m2 i to j coagulation
2875
2876       data (bm2ij (  1,  1,ibeta), ibeta = 1,10) /   &
2877         0.707107,  0.716828,  0.738240,  0.764827,  0.793610,   &
2878         0.822843,  0.851217,  0.877670,  0.901404,  0.921944/
2879       data (bm2ij (  1,  2,ibeta), ibeta = 1,10) /   &
2880         0.719180,  0.727975,  0.747638,  0.772334,  0.799234,   &
2881         0.826666,  0.853406,  0.878482,  0.901162,  0.920987/
2882       data (bm2ij (  1,  3,ibeta), ibeta = 1,10) /   &
2883         0.760947,  0.767874,  0.783692,  0.803890,  0.826015,   &
2884         0.848562,  0.870498,  0.891088,  0.909823,  0.926400/
2885       data (bm2ij (  1,  4,ibeta), ibeta = 1,10) /   &
2886         0.830926,  0.836034,  0.847708,  0.862528,  0.878521,   &
2887         0.894467,  0.909615,  0.923520,  0.935959,  0.946858/
2888       data (bm2ij (  1,  5,ibeta), ibeta = 1,10) /   &
2889         0.903643,  0.907035,  0.914641,  0.924017,  0.933795,   &
2890         0.943194,  0.951806,  0.959449,  0.966087,  0.971761/
2891       data (bm2ij (  1,  6,ibeta), ibeta = 1,10) /   &
2892         0.954216,  0.956094,  0.960211,  0.965123,  0.970068,   &
2893         0.974666,  0.978750,  0.982277,  0.985268,  0.987775/
2894       data (bm2ij (  1,  7,ibeta), ibeta = 1,10) /   &
2895         0.980546,  0.981433,  0.983343,  0.985568,  0.987751,   &
2896         0.989735,  0.991461,  0.992926,  0.994150,  0.995164/
2897       data (bm2ij (  1,  8,ibeta), ibeta = 1,10) /   &
2898         0.992142,  0.992524,  0.993338,  0.994272,  0.995174,   &
2899         0.995981,  0.996675,  0.997257,  0.997740,  0.998137/
2900       data (bm2ij (  1,  9,ibeta), ibeta = 1,10) /   &
2901         0.996868,  0.997026,  0.997361,  0.997742,  0.998106,   &
2902         0.998430,  0.998705,  0.998935,  0.999125,  0.999280/
2903       data (bm2ij (  1, 10,ibeta), ibeta = 1,10) /   &
2904         0.998737,  0.998802,  0.998939,  0.999094,  0.999241,   &
2905         0.999371,  0.999481,  0.999573,  0.999648,  0.999709/
2906       data (bm2ij (  2,  1,ibeta), ibeta = 1,10) /   &
2907         0.729600,  0.739948,  0.763059,  0.791817,  0.822510,   &
2908         0.852795,  0.881000,  0.905999,  0.927206,  0.944532/
2909       data (bm2ij (  2,  2,ibeta), ibeta = 1,10) /   &
2910         0.727025,  0.737116,  0.759615,  0.787657,  0.817740,   &
2911         0.847656,  0.875801,  0.901038,  0.922715,  0.940643/
2912       data (bm2ij (  2,  3,ibeta), ibeta = 1,10) /   &
2913         0.738035,  0.746779,  0.766484,  0.791340,  0.818324,   &
2914         0.845546,  0.871629,  0.895554,  0.916649,  0.934597/
2915       data (bm2ij (  2,  4,ibeta), ibeta = 1,10) /   &
2916         0.784185,  0.790883,  0.806132,  0.825501,  0.846545,   &
2917         0.867745,  0.888085,  0.906881,  0.923705,  0.938349/
2918       data (bm2ij (  2,  5,ibeta), ibeta = 1,10) /   &
2919         0.857879,  0.862591,  0.873238,  0.886539,  0.900645,   &
2920         0.914463,  0.927360,  0.939004,  0.949261,  0.958125/
2921       data (bm2ij (  2,  6,ibeta), ibeta = 1,10) /   &
2922         0.925441,  0.928304,  0.934645,  0.942324,  0.950181,   &
2923         0.957600,  0.964285,  0.970133,  0.975147,  0.979388/
2924       data (bm2ij (  2,  7,ibeta), ibeta = 1,10) /   &
2925         0.966728,  0.968176,  0.971323,  0.975027,  0.978705,   &
2926         0.982080,  0.985044,  0.987578,  0.989710,  0.991485/
2927       data (bm2ij (  2,  8,ibeta), ibeta = 1,10) /   &
2928         0.986335,  0.986980,  0.988362,  0.989958,  0.991511,   &
2929         0.992912,  0.994122,  0.995143,  0.995992,  0.996693/
2930       data (bm2ij (  2,  9,ibeta), ibeta = 1,10) /   &
2931         0.994547,  0.994817,  0.995391,  0.996046,  0.996677,   &
2932         0.997238,  0.997719,  0.998122,  0.998454,  0.998727/
2933       data (bm2ij (  2, 10,ibeta), ibeta = 1,10) /   &
2934         0.997817,  0.997928,  0.998163,  0.998429,  0.998683,   &
2935         0.998908,  0.999099,  0.999258,  0.999389,  0.999497/
2936       data (bm2ij (  3,  1,ibeta), ibeta = 1,10) /   &
2937         0.783612,  0.793055,  0.814468,  0.841073,  0.868769,   &
2938         0.894963,  0.918118,  0.937527,  0.953121,  0.965244/
2939       data (bm2ij (  3,  2,ibeta), ibeta = 1,10) /   &
2940         0.772083,  0.781870,  0.803911,  0.831238,  0.859802,   &
2941         0.887036,  0.911349,  0.931941,  0.948649,  0.961751/
2942       data (bm2ij (  3,  3,ibeta), ibeta = 1,10) /   &
2943         0.755766,  0.765509,  0.787380,  0.814630,  0.843526,   &
2944         0.871670,  0.897443,  0.919870,  0.938557,  0.953576/
2945       data (bm2ij (  3,  4,ibeta), ibeta = 1,10) /   &
2946         0.763816,  0.772145,  0.790997,  0.814784,  0.840434,   &
2947         0.865978,  0.890034,  0.911671,  0.930366,  0.945963/
2948       data (bm2ij (  3,  5,ibeta), ibeta = 1,10) /   &
2949         0.813597,  0.819809,  0.833889,  0.851618,  0.870640,   &
2950         0.889514,  0.907326,  0.923510,  0.937768,  0.950003/
2951       data (bm2ij (  3,  6,ibeta), ibeta = 1,10) /   &
2952         0.886317,  0.890437,  0.899643,  0.910955,  0.922730,   &
2953         0.934048,  0.944422,  0.953632,  0.961624,  0.968444/
2954       data (bm2ij (  3,  7,ibeta), ibeta = 1,10) /   &
2955         0.944565,  0.946855,  0.951872,  0.957854,  0.963873,   &
2956         0.969468,  0.974438,  0.978731,  0.982372,  0.985424/
2957       data (bm2ij (  3,  8,ibeta), ibeta = 1,10) /   &
2958         0.976358,  0.977435,  0.979759,  0.982467,  0.985125,   &
2959         0.987540,  0.989642,  0.991425,  0.992916,  0.994150/
2960       data (bm2ij (  3,  9,ibeta), ibeta = 1,10) /   &
2961         0.990471,  0.990932,  0.991917,  0.993048,  0.994142,   &
2962         0.995121,  0.995964,  0.996671,  0.997258,  0.997740/
2963       data (bm2ij (  3, 10,ibeta), ibeta = 1,10) /   &
2964         0.996199,  0.996389,  0.996794,  0.997254,  0.997694,   &
2965         0.998086,  0.998420,  0.998699,  0.998929,  0.999117/
2966       data (bm2ij (  4,  1,ibeta), ibeta = 1,10) /   &
2967         0.844355,  0.852251,  0.869914,  0.891330,  0.912823,   &
2968         0.932259,  0.948642,  0.961767,  0.971897,  0.979510/
2969       data (bm2ij (  4,  2,ibeta), ibeta = 1,10) /   &
2970         0.831550,  0.839954,  0.858754,  0.881583,  0.904592,   &
2971         0.925533,  0.943309,  0.957647,  0.968779,  0.977185/
2972       data (bm2ij (  4,  3,ibeta), ibeta = 1,10) /   &
2973         0.803981,  0.813288,  0.834060,  0.859400,  0.885285,   &
2974         0.909286,  0.930084,  0.947193,  0.960714,  0.971078/
2975       data (bm2ij (  4,  4,ibeta), ibeta = 1,10) /   &
2976         0.781787,  0.791080,  0.811931,  0.837749,  0.864768,   &
2977         0.890603,  0.913761,  0.933477,  0.949567,  0.962261/
2978       data (bm2ij (  4,  5,ibeta), ibeta = 1,10) /   &
2979         0.791591,  0.799355,  0.816916,  0.838961,  0.862492,   &
2980         0.885595,  0.907003,  0.925942,  0.942052,  0.955310/
2981       data (bm2ij (  4,  6,ibeta), ibeta = 1,10) /   &
2982         0.844933,  0.850499,  0.863022,  0.878593,  0.895038,   &
2983         0.911072,  0.925939,  0.939227,  0.950765,  0.960550/
2984       data (bm2ij (  4,  7,ibeta), ibeta = 1,10) /   &
2985         0.912591,  0.916022,  0.923607,  0.932777,  0.942151,   &
2986         0.951001,  0.958976,  0.965950,  0.971924,  0.976965/
2987       data (bm2ij (  4,  8,ibeta), ibeta = 1,10) /   &
2988         0.959859,  0.961617,  0.965433,  0.969924,  0.974382,   &
2989         0.978472,  0.982063,  0.985134,  0.987716,  0.989865/
2990       data (bm2ij (  4,  9,ibeta), ibeta = 1,10) /   &
2991         0.983377,  0.984162,  0.985844,  0.987788,  0.989681,   &
2992         0.991386,  0.992860,  0.994104,  0.995139,  0.995991/
2993       data (bm2ij (  4, 10,ibeta), ibeta = 1,10) /   &
2994         0.993343,  0.993672,  0.994370,  0.995169,  0.995937,   &
2995         0.996622,  0.997209,  0.997700,  0.998106,  0.998439/
2996       data (bm2ij (  5,  1,ibeta), ibeta = 1,10) /   &
2997         0.895806,  0.901918,  0.915233,  0.930783,  0.945768,   &
2998         0.958781,  0.969347,  0.977540,  0.983697,  0.988225/
2999       data (bm2ij (  5,  2,ibeta), ibeta = 1,10) /   &
3000         0.885634,  0.892221,  0.906629,  0.923540,  0.939918,   &
3001         0.954213,  0.965873,  0.974951,  0.981794,  0.986840/
3002       data (bm2ij (  5,  3,ibeta), ibeta = 1,10) /   &
3003         0.860120,  0.867858,  0.884865,  0.904996,  0.924724,   &
3004         0.942177,  0.956602,  0.967966,  0.976616,  0.983043/
3005       data (bm2ij (  5,  4,ibeta), ibeta = 1,10) /   &
3006         0.827462,  0.836317,  0.855885,  0.879377,  0.902897,   &
3007         0.924232,  0.942318,  0.956900,  0.968222,  0.976774/
3008       data (bm2ij (  5,  5,ibeta), ibeta = 1,10) /   &
3009         0.805527,  0.814279,  0.833853,  0.857892,  0.882726,   &
3010         0.906095,  0.926690,  0.943938,  0.957808,  0.968615/
3011       data (bm2ij (  5,  6,ibeta), ibeta = 1,10) /   &
3012         0.820143,  0.827223,  0.843166,  0.863002,  0.883905,   &
3013         0.904128,  0.922585,  0.938687,  0.952222,  0.963255/
3014       data (bm2ij (  5,  7,ibeta), ibeta = 1,10) /   &
3015         0.875399,  0.880208,  0.890929,  0.904065,  0.917699,   &
3016         0.930756,  0.942656,  0.953131,  0.962113,  0.969657/
3017       data (bm2ij (  5,  8,ibeta), ibeta = 1,10) /   &
3018         0.934782,  0.937520,  0.943515,  0.950656,  0.957840,   &
3019         0.964516,  0.970446,  0.975566,  0.979905,  0.983534/
3020       data (bm2ij (  5,  9,ibeta), ibeta = 1,10) /   &
3021         0.971369,  0.972679,  0.975505,  0.978797,  0.982029,   &
3022         0.984964,  0.987518,  0.989685,  0.991496,  0.992994/
3023       data (bm2ij (  5, 10,ibeta), ibeta = 1,10) /   &
3024         0.988329,  0.988893,  0.990099,  0.991485,  0.992825,   &
3025         0.994025,  0.995058,  0.995925,  0.996643,  0.997234/
3026       data (bm2ij (  6,  1,ibeta), ibeta = 1,10) /   &
3027         0.933384,  0.937784,  0.947130,  0.957655,  0.967430,   &
3028         0.975639,  0.982119,  0.987031,  0.990657,  0.993288/
3029       data (bm2ij (  6,  2,ibeta), ibeta = 1,10) /   &
3030         0.926445,  0.931227,  0.941426,  0.952975,  0.963754,   &
3031         0.972845,  0.980044,  0.985514,  0.989558,  0.992498/
3032       data (bm2ij (  6,  3,ibeta), ibeta = 1,10) /   &
3033         0.907835,  0.913621,  0.926064,  0.940308,  0.953745,   &
3034         0.965189,  0.974327,  0.981316,  0.986510,  0.990297/
3035       data (bm2ij (  6,  4,ibeta), ibeta = 1,10) /   &
3036         0.879088,  0.886306,  0.901945,  0.920079,  0.937460,   &
3037         0.952509,  0.964711,  0.974166,  0.981265,  0.986484/
3038       data (bm2ij (  6,  5,ibeta), ibeta = 1,10) /   &
3039         0.846500,  0.854862,  0.873189,  0.894891,  0.916264,   &
3040         0.935315,  0.951197,  0.963812,  0.973484,  0.980715/
3041       data (bm2ij (  6,  6,ibeta), ibeta = 1,10) /   &
3042         0.828137,  0.836250,  0.854310,  0.876287,  0.898710,   &
3043         0.919518,  0.937603,  0.952560,  0.964461,  0.973656/
3044       data (bm2ij (  6,  7,ibeta), ibeta = 1,10) /   &
3045         0.848595,  0.854886,  0.868957,  0.886262,  0.904241,   &
3046         0.921376,  0.936799,  0.950096,  0.961172,  0.970145/
3047       data (bm2ij (  6,  8,ibeta), ibeta = 1,10) /   &
3048         0.902919,  0.906922,  0.915760,  0.926427,  0.937312,   &
3049         0.947561,  0.956758,  0.964747,  0.971525,  0.977175/
3050       data (bm2ij (  6,  9,ibeta), ibeta = 1,10) /   &
3051         0.952320,  0.954434,  0.959021,  0.964418,  0.969774,   &
3052         0.974688,  0.979003,  0.982690,  0.985789,  0.988364/
3053       data (bm2ij (  6, 10,ibeta), ibeta = 1,10) /   &
3054         0.979689,  0.980650,  0.982712,  0.985093,  0.987413,   &
3055         0.989502,  0.991308,  0.992831,  0.994098,  0.995142/
3056       data (bm2ij (  7,  1,ibeta), ibeta = 1,10) /   &
3057         0.958611,  0.961598,  0.967817,  0.974620,  0.980752,   &
3058         0.985771,  0.989650,  0.992543,  0.994653,  0.996171/
3059       data (bm2ij (  7,  2,ibeta), ibeta = 1,10) /   &
3060         0.954225,  0.957488,  0.964305,  0.971795,  0.978576,   &
3061         0.984144,  0.988458,  0.991681,  0.994034,  0.995728/
3062       data (bm2ij (  7,  3,ibeta), ibeta = 1,10) /   &
3063         0.942147,  0.946158,  0.954599,  0.963967,  0.972529,   &
3064         0.979612,  0.985131,  0.989271,  0.992301,  0.994487/
3065       data (bm2ij (  7,  4,ibeta), ibeta = 1,10) /   &
3066         0.921821,  0.927048,  0.938140,  0.950598,  0.962118,   &
3067         0.971752,  0.979326,  0.985046,  0.989254,  0.992299/
3068       data (bm2ij (  7,  5,ibeta), ibeta = 1,10) /   &
3069         0.893419,  0.900158,  0.914598,  0.931070,  0.946584,   &
3070         0.959795,  0.970350,  0.978427,  0.984432,  0.988811/
3071       data (bm2ij (  7,  6,ibeta), ibeta = 1,10) /   &
3072         0.863302,  0.871111,  0.888103,  0.907990,  0.927305,   &
3073         0.944279,  0.958245,  0.969211,  0.977540,  0.983720/
3074       data (bm2ij (  7,  7,ibeta), ibeta = 1,10) /   &
3075         0.850182,  0.857560,  0.873890,  0.893568,  0.913408,   &
3076         0.931591,  0.947216,  0.960014,  0.970121,  0.977886/
3077       data (bm2ij (  7,  8,ibeta), ibeta = 1,10) /   &
3078         0.875837,  0.881265,  0.893310,  0.907936,  0.922910,   &
3079         0.936977,  0.949480,  0.960154,  0.968985,  0.976111/
3080       data (bm2ij (  7,  9,ibeta), ibeta = 1,10) /   &
3081         0.926228,  0.929445,  0.936486,  0.944868,  0.953293,   &
3082         0.961108,  0.968028,  0.973973,  0.978974,  0.983118/
3083       data (bm2ij (  7, 10,ibeta), ibeta = 1,10) /   &
3084         0.965533,  0.967125,  0.970558,  0.974557,  0.978484,   &
3085         0.982050,  0.985153,  0.987785,  0.989982,  0.991798/
3086       data (bm2ij (  8,  1,ibeta), ibeta = 1,10) /   &
3087         0.974731,  0.976674,  0.980660,  0.984926,  0.988689,   &
3088         0.991710,  0.994009,  0.995703,  0.996929,  0.997805/
3089       data (bm2ij (  8,  2,ibeta), ibeta = 1,10) /   &
3090         0.972062,  0.974192,  0.978571,  0.983273,  0.987432,   &
3091         0.990780,  0.993333,  0.995218,  0.996581,  0.997557/
3092       data (bm2ij (  8,  3,ibeta), ibeta = 1,10) /   &
3093         0.964662,  0.967300,  0.972755,  0.978659,  0.983921,   &
3094         0.988181,  0.991444,  0.993859,  0.995610,  0.996863/
3095       data (bm2ij (  8,  4,ibeta), ibeta = 1,10) /   &
3096         0.951782,  0.955284,  0.962581,  0.970559,  0.977737,   &
3097         0.983593,  0.988103,  0.991454,  0.993889,  0.995635/
3098       data (bm2ij (  8,  5,ibeta), ibeta = 1,10) /   &
3099         0.931947,  0.936723,  0.946751,  0.957843,  0.967942,   &
3100         0.976267,  0.982734,  0.987571,  0.991102,  0.993642/
3101       data (bm2ij (  8,  6,ibeta), ibeta = 1,10) /   &
3102         0.905410,  0.911665,  0.924950,  0.939908,  0.953798,   &
3103         0.965469,  0.974684,  0.981669,  0.986821,  0.990556/
3104       data (bm2ij (  8,  7,ibeta), ibeta = 1,10) /   &
3105         0.878941,  0.886132,  0.901679,  0.919688,  0.936970,   &
3106         0.951980,  0.964199,  0.973709,  0.980881,  0.986174/
3107       data (bm2ij (  8,  8,ibeta), ibeta = 1,10) /   &
3108         0.871653,  0.878218,  0.892652,  0.909871,  0.927034,   &
3109         0.942592,  0.955836,  0.966604,  0.975065,  0.981545/
3110       data (bm2ij (  8,  9,ibeta), ibeta = 1,10) /   &
3111         0.900693,  0.905239,  0.915242,  0.927232,  0.939335,   &
3112         0.950555,  0.960420,  0.968774,  0.975651,  0.981188/
3113       data (bm2ij (  8, 10,ibeta), ibeta = 1,10) /   &
3114         0.944922,  0.947435,  0.952894,  0.959317,  0.965689,   &
3115         0.971529,  0.976645,  0.981001,  0.984641,  0.987642/
3116       data (bm2ij (  9,  1,ibeta), ibeta = 1,10) /   &
3117         0.984736,  0.985963,  0.988453,  0.991078,  0.993357,   &
3118         0.995161,  0.996519,  0.997512,  0.998226,  0.998734/
3119       data (bm2ij (  9,  2,ibeta), ibeta = 1,10) /   &
3120         0.983141,  0.984488,  0.987227,  0.990119,  0.992636,   &
3121         0.994632,  0.996137,  0.997238,  0.998030,  0.998595/
3122       data (bm2ij (  9,  3,ibeta), ibeta = 1,10) /   &
3123         0.978726,  0.980401,  0.983819,  0.987450,  0.990626,   &
3124         0.993157,  0.995071,  0.996475,  0.997486,  0.998206/
3125       data (bm2ij (  9,  4,ibeta), ibeta = 1,10) /   &
3126         0.970986,  0.973224,  0.977818,  0.982737,  0.987072,   &
3127         0.990546,  0.993184,  0.995124,  0.996523,  0.997521/
3128       data (bm2ij (  9,  5,ibeta), ibeta = 1,10) /   &
3129         0.958579,  0.961700,  0.968149,  0.975116,  0.981307,   &
3130         0.986301,  0.990112,  0.992923,  0.994954,  0.996404/
3131       data (bm2ij (  9,  6,ibeta), ibeta = 1,10) /   &
3132         0.940111,  0.944479,  0.953572,  0.963506,  0.972436,   &
3133         0.979714,  0.985313,  0.989468,  0.992483,  0.994641/
3134       data (bm2ij (  9,  7,ibeta), ibeta = 1,10) /   &
3135         0.916127,  0.921878,  0.934003,  0.947506,  0.959899,   &
3136         0.970199,  0.978255,  0.984314,  0.988755,  0.991960/
3137       data (bm2ij (  9,  8,ibeta), ibeta = 1,10) /   &
3138         0.893848,  0.900364,  0.914368,  0.930438,  0.945700,   &
3139         0.958824,  0.969416,  0.977603,  0.983746,  0.988262/
3140       data (bm2ij (  9,  9,ibeta), ibeta = 1,10) /   &
3141         0.892161,  0.897863,  0.910315,  0.925021,  0.939523,   &
3142         0.952544,  0.963544,  0.972442,  0.979411,  0.984742/
3143       data (bm2ij (  9, 10,ibeta), ibeta = 1,10) /   &
3144         0.922260,  0.925966,  0.934047,  0.943616,  0.953152,   &
3145         0.961893,  0.969506,  0.975912,  0.981167,  0.985394/
3146       data (bm2ij ( 10,  1,ibeta), ibeta = 1,10) /   &
3147         0.990838,  0.991598,  0.993128,  0.994723,  0.996092,   &
3148         0.997167,  0.997969,  0.998552,  0.998969,  0.999265/
3149       data (bm2ij ( 10,  2,ibeta), ibeta = 1,10) /   &
3150         0.989892,  0.990727,  0.992411,  0.994167,  0.995678,   &
3151         0.996864,  0.997751,  0.998396,  0.998858,  0.999186/
3152       data (bm2ij ( 10,  3,ibeta), ibeta = 1,10) /   &
3153         0.987287,  0.988327,  0.990428,  0.992629,  0.994529,   &
3154         0.996026,  0.997148,  0.997965,  0.998551,  0.998967/
3155       data (bm2ij ( 10,  4,ibeta), ibeta = 1,10) /   &
3156         0.982740,  0.984130,  0.986952,  0.989926,  0.992508,   &
3157         0.994551,  0.996087,  0.997208,  0.998012,  0.998584/
3158       data (bm2ij ( 10,  5,ibeta), ibeta = 1,10) /   &
3159         0.975380,  0.977330,  0.981307,  0.985529,  0.989216,   &
3160         0.992147,  0.994358,  0.995975,  0.997136,  0.997961/
3161       data (bm2ij ( 10,  6,ibeta), ibeta = 1,10) /   &
3162         0.963911,  0.966714,  0.972465,  0.978614,  0.984022,   &
3163         0.988346,  0.991620,  0.994020,  0.995747,  0.996974/
3164       data (bm2ij ( 10,  7,ibeta), ibeta = 1,10) /   &
3165         0.947187,  0.951161,  0.959375,  0.968258,  0.976160,   &
3166         0.982540,  0.987409,  0.991000,  0.993592,  0.995441/
3167       data (bm2ij ( 10,  8,ibeta), ibeta = 1,10) /   &
3168         0.926045,  0.931270,  0.942218,  0.954297,  0.965273,   &
3169         0.974311,  0.981326,  0.986569,  0.990394,  0.993143/
3170       data (bm2ij ( 10,  9,ibeta), ibeta = 1,10) /   &
3171         0.908092,  0.913891,  0.926288,  0.940393,  0.953667,   &
3172         0.964987,  0.974061,  0.981038,  0.986253,  0.990078/
3173       data (bm2ij ( 10, 10,ibeta), ibeta = 1,10) /   &
3174         0.911143,  0.915972,  0.926455,  0.938721,  0.950701,   &
3175         0.961370,  0.970329,  0.977549,  0.983197,  0.987518/
3176
3177
3178 ! fsb total correction factor for m2 coagulation j from i
3179
3180       data  (bm2ji( 1, 1,ibeta), ibeta = 1,10) /   &
3181         0.753466,  0.756888,  0.761008,  0.759432,  0.748675,   &
3182         0.726951,  0.693964,  0.650915,  0.600227,  0.545000/
3183       data  (bm2ji( 1, 2,ibeta), ibeta = 1,10) /   &
3184         0.824078,  0.828698,  0.835988,  0.838943,  0.833454,   &
3185         0.817148,  0.789149,  0.750088,  0.701887,  0.647308/
3186       data  (bm2ji( 1, 3,ibeta), ibeta = 1,10) /   &
3187         1.007389,  1.014362,  1.028151,  1.041011,  1.047939,   &
3188         1.045707,  1.032524,  1.007903,  0.972463,  0.927667/
3189       data  (bm2ji( 1, 4,ibeta), ibeta = 1,10) /   &
3190         1.246157,  1.255135,  1.274249,  1.295351,  1.313362,   &
3191         1.325187,  1.329136,  1.324491,  1.311164,  1.289459/
3192       data  (bm2ji( 1, 5,ibeta), ibeta = 1,10) /   &
3193         1.450823,  1.459551,  1.478182,  1.499143,  1.518224,   &
3194         1.533312,  1.543577,  1.548882,  1.549395,  1.545364/
3195       data  (bm2ji( 1, 6,ibeta), ibeta = 1,10) /   &
3196         1.575248,  1.581832,  1.595643,  1.610866,  1.624601,   &
3197         1.635690,  1.643913,  1.649470,  1.652688,  1.653878/
3198       data  (bm2ji( 1, 7,ibeta), ibeta = 1,10) /   &
3199         1.638426,  1.642626,  1.651293,  1.660641,  1.668926,   &
3200         1.675571,  1.680572,  1.684147,  1.686561,  1.688047/
3201       data  (bm2ji( 1, 8,ibeta), ibeta = 1,10) /   &
3202         1.669996,  1.672392,  1.677283,  1.682480,  1.687028,   &
3203         1.690651,  1.693384,  1.695372,  1.696776,  1.697734/
3204       data  (bm2ji( 1, 9,ibeta), ibeta = 1,10) /   &
3205         1.686148,  1.687419,  1.689993,  1.692704,  1.695057,   &
3206         1.696922,  1.698329,  1.699359,  1.700099,  1.700621/
3207       data  (bm2ji( 1,10,ibeta), ibeta = 1,10) /   &
3208         1.694364,  1.695010,  1.696313,  1.697676,  1.698853,   &
3209         1.699782,  1.700482,  1.700996,  1.701366,  1.701631/
3210       data  (bm2ji( 2, 1,ibeta), ibeta = 1,10) /   &
3211         0.783166,  0.779369,  0.768044,  0.747572,  0.716709,   &
3212         0.675422,  0.624981,  0.567811,  0.507057,  0.445975/
3213       data  (bm2ji( 2, 2,ibeta), ibeta = 1,10) /   &
3214         0.848390,  0.847100,  0.840874,  0.826065,  0.800296,   &
3215         0.762625,  0.713655,  0.655545,  0.591603,  0.525571/
3216       data  (bm2ji( 2, 3,ibeta), ibeta = 1,10) /   &
3217         1.039894,  1.043786,  1.049445,  1.049664,  1.039407,   &
3218         1.015322,  0.975983,  0.922180,  0.856713,  0.783634/
3219       data  (bm2ji( 2, 4,ibeta), ibeta = 1,10) /   &
3220         1.345995,  1.356064,  1.376947,  1.398304,  1.412685,   &
3221         1.414611,  1.400652,  1.369595,  1.322261,  1.260993/
3222       data  (bm2ji( 2, 5,ibeta), ibeta = 1,10) /   &
3223         1.675575,  1.689859,  1.720957,  1.756659,  1.788976,   &
3224         1.812679,  1.824773,  1.824024,  1.810412,  1.784630/
3225       data  (bm2ji( 2, 6,ibeta), ibeta = 1,10) /   &
3226         1.919835,  1.933483,  1.962973,  1.996810,  2.028377,   &
3227         2.054172,  2.072763,  2.083963,  2.088190,  2.086052/
3228       data  (bm2ji( 2, 7,ibeta), ibeta = 1,10) /   &
3229         2.064139,  2.074105,  2.095233,  2.118909,  2.140688,   &
3230         2.158661,  2.172373,  2.182087,  2.188330,  2.191650/
3231       data  (bm2ji( 2, 8,ibeta), ibeta = 1,10) /   &
3232         2.144871,  2.150990,  2.163748,  2.177731,  2.190364,   &
3233         2.200712,  2.208687,  2.214563,  2.218716,  2.221502/
3234       data  (bm2ji( 2, 9,ibeta), ibeta = 1,10) /   &
3235         2.189223,  2.192595,  2.199540,  2.207033,  2.213706,   &
3236         2.219125,  2.223297,  2.226403,  2.228660,  2.230265/
3237       data  (bm2ji( 2,10,ibeta), ibeta = 1,10) /   &
3238         2.212595,  2.214342,  2.217912,  2.221723,  2.225082,   &
3239         2.227791,  2.229869,  2.231417,  2.232551,  2.233372/
3240       data  (bm2ji( 3, 1,ibeta), ibeta = 1,10) /   &
3241         0.837870,  0.824476,  0.793119,  0.750739,  0.700950,   &
3242         0.646691,  0.590508,  0.534354,  0.479532,  0.426856/
3243       data  (bm2ji( 3, 2,ibeta), ibeta = 1,10) /   &
3244         0.896771,  0.885847,  0.859327,  0.821694,  0.775312,   &
3245         0.722402,  0.665196,  0.605731,  0.545742,  0.486687/
3246       data  (bm2ji( 3, 3,ibeta), ibeta = 1,10) /   &
3247         1.076089,  1.071727,  1.058845,  1.036171,  1.002539,   &
3248         0.957521,  0.901640,  0.836481,  0.764597,  0.689151/
3249       data  (bm2ji( 3, 4,ibeta), ibeta = 1,10) /   &
3250         1.409571,  1.415168,  1.425346,  1.432021,  1.428632,   &
3251         1.409696,  1.371485,  1.312958,  1.236092,  1.145293/
3252       data  (bm2ji( 3, 5,ibeta), ibeta = 1,10) /   &
3253         1.862757,  1.880031,  1.918394,  1.963456,  2.004070,   &
3254         2.030730,  2.036144,  2.016159,  1.970059,  1.900079/
3255       data  (bm2ji( 3, 6,ibeta), ibeta = 1,10) /   &
3256         2.289741,  2.313465,  2.366789,  2.431612,  2.495597,   &
3257         2.549838,  2.588523,  2.608665,  2.609488,  2.591662/
3258       data  (bm2ji( 3, 7,ibeta), ibeta = 1,10) /   &
3259         2.597157,  2.618731,  2.666255,  2.722597,  2.777531,   &
3260         2.825187,  2.862794,  2.889648,  2.906199,  2.913380/
3261       data  (bm2ji( 3, 8,ibeta), ibeta = 1,10) /   &
3262         2.797975,  2.813116,  2.845666,  2.882976,  2.918289,   &
3263         2.948461,  2.972524,  2.990687,  3.003664,  3.012284/
3264       data  (bm2ji( 3, 9,ibeta), ibeta = 1,10) /   &
3265         2.920832,  2.929843,  2.948848,  2.970057,  2.989632,   &
3266         3.006057,  3.019067,  3.028979,  3.036307,  3.041574/
3267       data  (bm2ji( 3,10,ibeta), ibeta = 1,10) /   &
3268         2.989627,  2.994491,  3.004620,  3.015720,  3.025789,   &
3269         3.034121,  3.040664,  3.045641,  3.049347,  3.052066/
3270       data  (bm2ji( 4, 1,ibeta), ibeta = 1,10) /   &
3271         0.893179,  0.870897,  0.820996,  0.759486,  0.695488,   &
3272         0.634582,  0.579818,  0.532143,  0.490927,  0.454618/
3273       data  (bm2ji( 4, 2,ibeta), ibeta = 1,10) /   &
3274         0.948355,  0.927427,  0.880215,  0.821146,  0.758524,   &
3275         0.697680,  0.641689,  0.591605,  0.546919,  0.506208/
3276       data  (bm2ji( 4, 3,ibeta), ibeta = 1,10) /   &
3277         1.109562,  1.093648,  1.056438,  1.007310,  0.951960,   &
3278         0.894453,  0.837364,  0.781742,  0.727415,  0.673614/
3279       data  (bm2ji( 4, 4,ibeta), ibeta = 1,10) /   &
3280         1.423321,  1.417557,  1.402442,  1.379079,  1.347687,   &
3281         1.308075,  1.259703,  1.201983,  1.134778,  1.058878/
3282       data  (bm2ji( 4, 5,ibeta), ibeta = 1,10) /   &
3283         1.933434,  1.944347,  1.968765,  1.997653,  2.023054,   &
3284         2.036554,  2.029949,  1.996982,  1.934982,  1.845473/
3285       data  (bm2ji( 4, 6,ibeta), ibeta = 1,10) /   &
3286         2.547772,  2.577105,  2.645918,  2.735407,  2.830691,   &
3287         2.917268,  2.981724,  3.013684,  3.007302,  2.961560/
3288       data  (bm2ji( 4, 7,ibeta), ibeta = 1,10) /   &
3289         3.101817,  3.139271,  3.225851,  3.336402,  3.453409,   &
3290         3.563116,  3.655406,  3.724014,  3.766113,  3.781394/
3291       data  (bm2ji( 4, 8,ibeta), ibeta = 1,10) /   &
3292         3.540920,  3.573780,  3.647439,  3.737365,  3.828468,   &
3293         3.911436,  3.981317,  4.036345,  4.076749,  4.103751/
3294       data  (bm2ji( 4, 9,ibeta), ibeta = 1,10) /   &
3295         3.856771,  3.879363,  3.928579,  3.986207,  4.042173,   &
3296         4.091411,  4.132041,  4.164052,  4.188343,  4.206118/
3297       data  (bm2ji( 4,10,ibeta), ibeta = 1,10) /   &
3298         4.053923,  4.067191,  4.095509,  4.127698,  4.158037,   &
3299         4.184055,  4.205135,  4.221592,  4.234115,  4.243463/
3300       data  (bm2ji( 5, 1,ibeta), ibeta = 1,10) /   &
3301         0.935846,  0.906814,  0.843358,  0.768710,  0.695885,   &
3302         0.631742,  0.579166,  0.538471,  0.508410,  0.486863/
3303       data  (bm2ji( 5, 2,ibeta), ibeta = 1,10) /   &
3304         0.988308,  0.959524,  0.896482,  0.821986,  0.748887,   &
3305         0.684168,  0.630908,  0.589516,  0.558676,  0.536056/
3306       data  (bm2ji( 5, 3,ibeta), ibeta = 1,10) /   &
3307         1.133795,  1.107139,  1.048168,  0.977258,  0.906341,   &
3308         0.842477,  0.789093,  0.746731,  0.713822,  0.687495/
3309       data  (bm2ji( 5, 4,ibeta), ibeta = 1,10) /   &
3310         1.405692,  1.385781,  1.340706,  1.284776,  1.227085,   &
3311         1.173532,  1.127008,  1.087509,  1.052712,  1.018960/
3312       data  (bm2ji( 5, 5,ibeta), ibeta = 1,10) /   &
3313         1.884992,  1.879859,  1.868463,  1.854995,  1.841946,   &
3314         1.829867,  1.816972,  1.799319,  1.771754,  1.729406/
3315       data  (bm2ji( 5, 6,ibeta), ibeta = 1,10) /   &
3316         2.592275,  2.612268,  2.661698,  2.731803,  2.815139,   &
3317         2.901659,  2.978389,  3.031259,  3.048045,  3.021122/
3318       data  (bm2ji( 5, 7,ibeta), ibeta = 1,10) /   &
3319         3.390321,  3.435519,  3.545615,  3.698419,  3.876958,   &
3320         4.062790,  4.236125,  4.378488,  4.475619,  4.519170/
3321       data  (bm2ji( 5, 8,ibeta), ibeta = 1,10) /   &
3322         4.161376,  4.216558,  4.346896,  4.519451,  4.711107,   &
3323         4.902416,  5.077701,  5.226048,  5.341423,  5.421764/
3324       data  (bm2ji( 5, 9,ibeta), ibeta = 1,10) /   &
3325         4.843961,  4.892035,  5.001492,  5.138515,  5.281684,   &
3326         5.416805,  5.535493,  5.634050,  5.712063,  5.770996/
3327       data  (bm2ji( 5,10,ibeta), ibeta = 1,10) /   &
3328         5.352093,  5.385119,  5.458056,  5.545311,  5.632162,   &
3329         5.710566,  5.777005,  5.830863,  5.873123,  5.905442/
3330       data  (bm2ji( 6, 1,ibeta), ibeta = 1,10) /   &
3331         0.964038,  0.930794,  0.859433,  0.777776,  0.700566,   &
3332         0.634671,  0.582396,  0.543656,  0.517284,  0.501694/
3333       data  (bm2ji( 6, 2,ibeta), ibeta = 1,10) /   &
3334         1.013416,  0.979685,  0.907197,  0.824135,  0.745552,   &
3335         0.678616,  0.625870,  0.587348,  0.561864,  0.547674/
3336       data  (bm2ji( 6, 3,ibeta), ibeta = 1,10) /   &
3337         1.145452,  1.111457,  1.038152,  0.953750,  0.873724,   &
3338         0.805955,  0.753621,  0.717052,  0.694920,  0.684910/
3339       data  (bm2ji( 6, 4,ibeta), ibeta = 1,10) /   &
3340         1.376547,  1.345004,  1.276415,  1.196704,  1.121091,   &
3341         1.058249,  1.012197,  0.983522,  0.970323,  0.968933/
3342       data  (bm2ji( 6, 5,ibeta), ibeta = 1,10) /   &
3343         1.778801,  1.755897,  1.706074,  1.649008,  1.597602,   &
3344         1.560087,  1.540365,  1.538205,  1.549738,  1.568333/
3345       data  (bm2ji( 6, 6,ibeta), ibeta = 1,10) /   &
3346         2.447603,  2.445172,  2.443762,  2.451842,  2.475877,   &
3347         2.519039,  2.580118,  2.653004,  2.727234,  2.789738/
3348       data  (bm2ji( 6, 7,ibeta), ibeta = 1,10) /   &
3349         3.368490,  3.399821,  3.481357,  3.606716,  3.772101,   &
3350         3.969416,  4.184167,  4.396163,  4.582502,  4.721838/
3351       data  (bm2ji( 6, 8,ibeta), ibeta = 1,10) /   &
3352         4.426458,  4.489861,  4.648250,  4.877510,  5.160698,   &
3353         5.477495,  5.803123,  6.111250,  6.378153,  6.586050/
3354       data  (bm2ji( 6, 9,ibeta), ibeta = 1,10) /   &
3355         5.568061,  5.644988,  5.829837,  6.081532,  6.371214,   &
3356         6.672902,  6.963737,  7.226172,  7.449199,  7.627886/
3357       data  (bm2ji( 6,10,ibeta), ibeta = 1,10) /   &
3358         6.639152,  6.707020,  6.863974,  7.065285,  7.281744,   &
3359         7.492437,  7.683587,  7.847917,  7.983296,  8.090977/
3360       data  (bm2ji( 7, 1,ibeta), ibeta = 1,10) /   &
3361         0.980853,  0.945724,  0.871244,  0.787311,  0.708818,   &
3362         0.641987,  0.588462,  0.547823,  0.518976,  0.500801/
3363       data  (bm2ji( 7, 2,ibeta), ibeta = 1,10) /   &
3364         1.026738,  0.990726,  0.914306,  0.828140,  0.747637,   &
3365         0.679351,  0.625127,  0.584662,  0.556910,  0.540749/
3366       data  (bm2ji( 7, 3,ibeta), ibeta = 1,10) /   &
3367         1.146496,  1.108808,  1.028695,  0.938291,  0.854101,   &
3368         0.783521,  0.728985,  0.690539,  0.667272,  0.657977/
3369       data  (bm2ji( 7, 4,ibeta), ibeta = 1,10) /   &
3370         1.344846,  1.306434,  1.224543,  1.132031,  1.046571,   &
3371         0.976882,  0.926488,  0.896067,  0.884808,  0.891027/
3372       data  (bm2ji( 7, 5,ibeta), ibeta = 1,10) /   &
3373         1.670227,  1.634583,  1.558421,  1.472939,  1.396496,   &
3374         1.339523,  1.307151,  1.300882,  1.319622,  1.360166/
3375       data  (bm2ji( 7, 6,ibeta), ibeta = 1,10) /   &
3376         2.224548,  2.199698,  2.148284,  2.095736,  2.059319,   &
3377         2.050496,  2.075654,  2.136382,  2.229641,  2.347958/
3378       data  (bm2ji( 7, 7,ibeta), ibeta = 1,10) /   &
3379         3.104483,  3.105947,  3.118398,  3.155809,  3.230427,   &
3380         3.350585,  3.519071,  3.731744,  3.976847,  4.235616/
3381       data  (bm2ji( 7, 8,ibeta), ibeta = 1,10) /   &
3382         4.288426,  4.331456,  4.447024,  4.633023,  4.891991,   &
3383         5.221458,  5.610060,  6.036467,  6.471113,  6.880462/
3384       data  (bm2ji( 7, 9,ibeta), ibeta = 1,10) /   &
3385         5.753934,  5.837061,  6.048530,  6.363800,  6.768061,   &
3386         7.241280,  7.755346,  8.276666,  8.771411,  9.210826/
3387       data  (bm2ji( 7,10,ibeta), ibeta = 1,10) /   &
3388         7.466219,  7.568810,  7.819032,  8.168340,  8.582973,   &
3389         9.030174,  9.478159,  9.899834, 10.275940, 10.595910/
3390       data  (bm2ji( 8, 1,ibeta), ibeta = 1,10) /   &
3391         0.990036,  0.954782,  0.880531,  0.797334,  0.719410,   &
3392         0.652220,  0.596923,  0.552910,  0.519101,  0.494529/
3393       data  (bm2ji( 8, 2,ibeta), ibeta = 1,10) /   &
3394         1.032428,  0.996125,  0.919613,  0.833853,  0.753611,   &
3395         0.684644,  0.628260,  0.583924,  0.550611,  0.527407/
3396       data  (bm2ji( 8, 3,ibeta), ibeta = 1,10) /   &
3397         1.141145,  1.102521,  1.021017,  0.929667,  0.844515,   &
3398         0.772075,  0.714086,  0.670280,  0.639824,  0.621970/
3399       data  (bm2ji( 8, 4,ibeta), ibeta = 1,10) /   &
3400         1.314164,  1.273087,  1.186318,  1.089208,  0.999476,   &
3401         0.924856,  0.867948,  0.829085,  0.807854,  0.803759/
3402       data  (bm2ji( 8, 5,ibeta), ibeta = 1,10) /   &
3403         1.580611,  1.538518,  1.449529,  1.350459,  1.260910,   &
3404         1.190526,  1.143502,  1.121328,  1.124274,  1.151974/
3405       data  (bm2ji( 8, 6,ibeta), ibeta = 1,10) /   &
3406         2.016773,  1.977721,  1.895727,  1.806974,  1.732891,   &
3407         1.685937,  1.673026,  1.697656,  1.761039,  1.862391/
3408       data  (bm2ji( 8, 7,ibeta), ibeta = 1,10) /   &
3409         2.750093,  2.723940,  2.672854,  2.628264,  2.612250,   &
3410         2.640406,  2.723211,  2.866599,  3.071893,  3.335217/
3411       data  (bm2ji( 8, 8,ibeta), ibeta = 1,10) /   &
3412         3.881905,  3.887143,  3.913667,  3.981912,  4.111099,   &
3413         4.316575,  4.608146,  4.988157,  5.449592,  5.974848/
3414       data  (bm2ji( 8, 9,ibeta), ibeta = 1,10) /   &
3415         5.438870,  5.492742,  5.640910,  5.886999,  6.241641,   &
3416         6.710609,  7.289480,  7.960725,  8.693495,  9.446644/
3417       data  (bm2ji( 8,10,ibeta), ibeta = 1,10) /   &
3418         7.521152,  7.624621,  7.892039,  8.300444,  8.839787,   &
3419         9.493227, 10.231770, 11.015642, 11.799990, 12.542260/
3420       data  (bm2ji( 9, 1,ibeta), ibeta = 1,10) /   &
3421         0.994285,  0.960012,  0.887939,  0.807040,  0.730578,   &
3422         0.663410,  0.606466,  0.559137,  0.520426,  0.489429/
3423       data  (bm2ji( 9, 2,ibeta), ibeta = 1,10) /   &
3424         1.033505,  0.998153,  0.923772,  0.840261,  0.761383,   &
3425         0.692242,  0.633873,  0.585709,  0.546777,  0.516215/
3426       data  (bm2ji( 9, 3,ibeta), ibeta = 1,10) /   &
3427         1.132774,  1.094907,  1.015161,  0.925627,  0.841293,   &
3428         0.767888,  0.706741,  0.657439,  0.619135,  0.591119/
3429       data  (bm2ji( 9, 4,ibeta), ibeta = 1,10) /   &
3430         1.286308,  1.245273,  1.158809,  1.061889,  0.971208,   &
3431         0.893476,  0.830599,  0.782561,  0.748870,  0.729198/
3432       data  (bm2ji( 9, 5,ibeta), ibeta = 1,10) /   &
3433         1.511105,  1.467141,  1.374520,  1.271162,  1.175871,   &
3434         1.096887,  1.037243,  0.997820,  0.978924,  0.980962/
3435       data  (bm2ji( 9, 6,ibeta), ibeta = 1,10) /   &
3436         1.857468,  1.812177,  1.717002,  1.612197,  1.519171,   &
3437         1.448660,  1.405871,  1.393541,  1.413549,  1.467532/
3438       data  (bm2ji( 9, 7,ibeta), ibeta = 1,10) /   &
3439         2.430619,  2.388452,  2.301326,  2.210241,  2.139724,   &
3440         2.104571,  2.114085,  2.174696,  2.291294,  2.467500/
3441       data  (bm2ji( 9, 8,ibeta), ibeta = 1,10) /   &
3442         3.385332,  3.357690,  3.306611,  3.269804,  3.274462,   &
3443         3.340862,  3.484609,  3.717740,  4.048748,  4.481588/
3444       data  (bm2ji( 9, 9,ibeta), ibeta = 1,10) /   &
3445         4.850497,  4.858280,  4.896008,  4.991467,  5.171511,   &
3446         5.459421,  5.873700,  6.426128,  7.119061,  7.942603/
3447       data  (bm2ji( 9,10,ibeta), ibeta = 1,10) /   &
3448         6.957098,  7.020164,  7.197272,  7.499331,  7.946554,   &
3449         8.555048,  9.330503, 10.263610, 11.327454, 12.478332/
3450       data  (bm2ji(10, 1,ibeta), ibeta = 1,10) /   &
3451         0.994567,  0.961842,  0.892854,  0.814874,  0.740198,   &
3452         0.673303,  0.615105,  0.565139,  0.522558,  0.486556/
3453       data  (bm2ji(10, 2,ibeta), ibeta = 1,10) /   &
3454         1.031058,  0.997292,  0.926082,  0.845571,  0.768501,   &
3455         0.699549,  0.639710,  0.588538,  0.545197,  0.508894/
3456       data  (bm2ji(10, 3,ibeta), ibeta = 1,10) /   &
3457         1.122535,  1.086287,  1.009790,  0.923292,  0.840626,   &
3458         0.766982,  0.703562,  0.650004,  0.605525,  0.569411/
3459       data  (bm2ji(10, 4,ibeta), ibeta = 1,10) /   &
3460         1.261142,  1.221555,  1.137979,  1.043576,  0.953745,   &
3461         0.874456,  0.807292,  0.752109,  0.708326,  0.675477/
3462       data  (bm2ji(10, 5,ibeta), ibeta = 1,10) /   &
3463         1.456711,  1.413432,  1.322096,  1.219264,  1.122319,   &
3464         1.038381,  0.969743,  0.916811,  0.879544,  0.858099/
3465       data  (bm2ji(10, 6,ibeta), ibeta = 1,10) /   &
3466         1.741792,  1.695157,  1.596897,  1.487124,  1.385734,   &
3467         1.301670,  1.238638,  1.198284,  1.181809,  1.190689/
3468       data  (bm2ji(10, 7,ibeta), ibeta = 1,10) /   &
3469         2.190197,  2.141721,  2.040226,  1.929245,  1.832051,   &
3470         1.760702,  1.721723,  1.719436,  1.757705,  1.840677/
3471       data  (bm2ji(10, 8,ibeta), ibeta = 1,10) /   &
3472         2.940764,  2.895085,  2.801873,  2.707112,  2.638603,   &
3473         2.613764,  2.644686,  2.741255,  2.912790,  3.168519/
3474       data  (bm2ji(10, 9,ibeta), ibeta = 1,10) /   &
3475         4.186191,  4.155844,  4.101953,  4.069102,  4.089886,   &
3476         4.189530,  4.389145,  4.707528,  5.161567,  5.765283/
3477       data  (bm2ji(10,10,ibeta), ibeta = 1,10) /   &
3478         6.119526,  6.127611,  6.171174,  6.286528,  6.508738,   &
3479         6.869521,  7.396912,  8.113749,  9.034683, 10.162190/
3480
3481 ! *** end of data statements.
3482
3483
3484 ! *** start calculations:
3485
3486       constii = abs( half * ( two ) ** two3rds - one )
3487       sqrttwo = sqrt(two)
3488       dlgsqt2 = one / log( sqrttwo )
3489
3490          esat01   = exp( 0.125d0 * xxlsgat * xxlsgat )
3491          esac01   = exp( 0.125d0 * xxlsgac * xxlsgac )
3492
3493          esat04  = esat01 ** 4
3494          esac04  = esac01 ** 4
3495
3496          esat05  = esat04 * esat01
3497          esac05  = esac04 * esac01
3498
3499          esat08  = esat04 * esat04
3500          esac08  = esac04 * esac04
3501
3502          esat09  = esat08 * esat01
3503          esac09  = esac08 * esac01
3504
3505          esat16  = esat08 * esat08
3506          esac16  = esac08 * esac08
3507
3508          esat20  = esat16 * esat04
3509          esac20  = esac16 * esac04
3510
3511          esat24  = esat20 * esat04
3512          esac24  = esac20 * esac04
3513
3514          esat25  = esat20 * esat05
3515          esac25  = esac20 * esac05
3516
3517          esat36  = esat20 * esat16
3518          esac36  = esac20 * esac16
3519
3520          esat49  = esat24 * esat25
3521
3522          esat64  = esat20 * esat20 * esat24
3523          esac64  = esac20 * esac20 * esac24
3524
3525          esat100 = esat64 * esat36
3526
3527          dgat2   = dgatk * dgatk
3528          dgat3   = dgatk * dgatk * dgatk
3529          dgac2   = dgacc * dgacc
3530          dgac3   = dgacc * dgacc * dgacc
3531
3532          sqdgat  = sqrt( dgatk )
3533          sqdgac  = sqrt( dgacc )
3534          sqdgat5 = dgat2 * sqdgat
3535          sqdgac5 = dgac2 * sqdgac
3536          sqdgat7 = dgat3 * sqdgat
3537
3538          xm2at = dgat2 * esat16
3539          xm3at = dgat3 * esat36
3540
3541          xm2ac = dgac2 * esac16
3542          xm3ac = dgac3 * esac36
3543
3544 ! *** for the free molecular regime:  page h.3 of whitby et al. (1991)
3545
3546          r       = sqdgac / sqdgat
3547          r2      = r * r
3548          r3      = r2 * r
3549          rx4     = r2 * r2
3550          r5      = r3 * r2
3551          r6      = r3 * r3
3552          rx8      = rx4 * rx4
3553          ri1     = one / r
3554          ri2     = one / r2
3555          ri3     = one / r3
3556          ri4     = ri2 * ri2
3557          kngat   = two * lamda / dgatk
3558          kngac   = two * lamda / dgacc
3559
3560
3561 ! *** calculate ratio of geometric mean diameters
3562          rat = dgacc / dgatk
3563 ! *** trap subscripts for bm0 and bm0i, between 1 and 10
3564 !     see page h.5 of whitby et al. (1991)
3565
3566       n2n = max( 1, min( 10,   &
3567             nint( 4.0 * ( sgatk - 0.75d0 ) ) ) )
3568
3569       n2a = max( 1, min( 10,   &
3570             nint( 4.0 * ( sgacc - 0.75d0 ) ) ) )
3571
3572       n1  = max( 1, min( 10,   &
3573              1 + nint( dlgsqt2 * log( rat ) ) ) )
3574
3575 ! *** intermodal coagulation
3576
3577
3578 ! *** set up for zeroeth moment
3579
3580 ! *** near-continuum form:  equation h.10a of whitby et al. (1991)
3581
3582          coagnc0 = knc * (   &
3583           two + a * ( kngat * ( esat04 + r2 * esat16 * esac04 )   &
3584                     + kngac * ( esac04 + ri2 * esac16 * esat04 ) )   &
3585                     + ( r2 + ri2 ) * esat04 * esac04  )
3586
3587
3588 ! *** free-molecular form:  equation h.7a of whitby et al. (1991)
3589
3590          coagfm0 = kfmatac * sqdgat * bm0ij(n1,n2n,n2a) * (   &
3591                    esat01 + r * esac01 + two * r2 * esat01 * esac04   &
3592                  + rx4 * esat09 * esac16 + ri3 * esat16 * esac09   &
3593                  + two * ri1 * esat04 + esac01  )
3594
3595
3596 ! *** loss to accumulation mode
3597
3598 ! *** harmonic mean
3599
3600       coagatac0 = coagnc0 * coagfm0 / ( coagnc0 + coagfm0 )
3601
3602       qn12 = coagatac0
3603
3604
3605 ! *** set up for second moment
3606 !      the second moment equations are new and begin with equations a1
3607 !     through a4 of binkowski and shankar (1995). after some algebraic
3608 !     rearrangement and application of the extended mean value theorem
3609 !     of integral calculus, equations are obtained that can be solved
3610 !     analytically with correction factors as has been done by
3611 !     whitby et al. (1991)
3612
3613 ! *** the term ( dp1 + dp2 ) ** (2/3) in equations a3 and a4 of
3614 !     binkowski and shankar (1995) is approximated by
3615 !     (dgat ** 3 + dgac **3 ) ** 2/3
3616
3617 ! *** near-continuum form
3618
3619       i1nc = knc * dgat2 * (   &
3620              two * esat16   &
3621            + r2 * esat04 * esac04   &
3622            + ri2 * esat36 * esac04   &
3623            + a * kngat * (   &
3624                  esat04   &
3625            +     ri2 * esat16 * esac04   &
3626            +     ri4 * esat36 * esac16   &
3627            +     r2 * esac04 )  )
3628
3629
3630
3631
3632 ! *** free-molecular form
3633
3634        i1fm =  kfmatac * sqdgat5 * bm2ij(n1,n2n,n2a) * (   &
3635                esat25   &
3636             +  two * r2 * esat09 * esac04   &
3637             +  rx4 * esat01 * esac16   &
3638             +  ri3 * esat64 * esac09   &
3639             +  two * ri1 * esat36 * esac01   &
3640             +  r * esat16 * esac01  )
3641
3642
3643
3644 ! *** loss to accumulation mode
3645
3646 ! *** harmonic mean
3647
3648       i1 = ( i1fm * i1nc ) / ( i1fm + i1nc )
3649
3650       coagatac2 = i1
3651
3652       qs12 = coagatac2
3653
3654
3655 ! *** gain by accumulation mode
3656
3657       coagacat2 = ( ( one + r6 ) ** two3rds - rx4 ) * i1
3658
3659       qs21 = coagacat2 * bm2ji(n1,n2n,n2a)
3660
3661 ! *** set up for third moment
3662
3663 ! *** near-continuum form: equation h.10b of whitby et al. (1991)
3664
3665       coagnc3 = knc * dgat3 * (   &
3666                 two * esat36   &
3667               + a * kngat * ( esat16 + r2 * esat04 * esac04 )   &
3668               + a * kngac * ( esat36 * esac04 + ri2 * esat64 * esac16 )   &
3669               + r2 * esat16 * esac04 + ri2 * esat64 * esac04 )
3670
3671
3672 ! *** free_molecular form: equation h.7b of whitby et al. (1991)
3673
3674       coagfm3 = kfmatac * sqdgat7 * bm3i( n1, n2n, n2a ) * (   &
3675                esat49   &
3676               +  r * esat36  * esac01   &
3677               + two * r2 * esat25  * esac04   &
3678               + rx4 * esat09  * esac16   &
3679               + ri3 * esat100 * esac09   &
3680               + two * ri1 * esat64  * esac01 )
3681
3682 ! *** gain by accumulation mode = loss from aitken mode
3683
3684 ! *** harmonic mean
3685
3686       coagatac3 = coagnc3 * coagfm3 / ( coagnc3 + coagfm3 )
3687
3688       qv12 = coagatac3
3689
3690 ! *** intramodal coagulation
3691
3692 ! *** zeroeth moment
3693
3694 ! *** aitken mode
3695
3696 ! *** near-continuum form: equation h.12a of whitby et al. (1991)
3697
3698       coagnc_at = knc * (one + esat08 + a * kngat * (esat20 + esat04))
3699
3700 ! *** free-molecular form: equation h.11a of whitby et al. (1991)
3701
3702       coagfm_at = kfmat * sqdgat * bm0(n2n) *   &
3703                  ( esat01 + esat25 + two * esat05 )
3704
3705
3706 ! *** harmonic mean
3707
3708       coagatat0 = coagfm_at * coagnc_at / ( coagfm_at + coagnc_at )
3709
3710       qn11 = coagatat0
3711
3712
3713 ! *** accumulation mode
3714
3715 ! *** near-continuum form: equation h.12a of whitby et al. (1991)
3716
3717       coagnc_ac = knc * (one + esac08 + a * kngac * (esac20 + esac04))
3718
3719 ! *** free-molecular form: equation h.11a of whitby et al. (1991)
3720
3721       coagfm_ac = kfmac * sqdgac * bm0(n2a) *   &
3722                    ( esac01 + esac25 + two * esac05 )
3723
3724 ! *** harmonic mean
3725
3726       coagacac0 = coagfm_ac * coagnc_ac / ( coagfm_ac + coagnc_ac )
3727
3728       qn22 = coagacac0
3729
3730
3731 ! *** set up for second moment
3732 !      the second moment equations are new and begin with 3.11a on page
3733 !     45 of whitby et al. (1991). after some algebraic rearrangement and
3734 !     application of the extended mean value theorem of integral calculus
3735 !     equations are obtained that can be solved analytically with
3736 !     correction factors as has been done by whitby et al. (1991)
3737
3738 ! *** aitken mode
3739
3740 ! *** near-continuum
3741
3742       i1nc_at = knc * dgat2 * (   &
3743              two * esat16   &
3744            + esat04 * esat04   &
3745            + esat36 * esat04   &
3746            + a * kngat * (   &
3747                 two * esat04   &
3748            +     esat16 * esat04   &
3749            +     esat36 * esat16 )  )
3750
3751 ! *** free- molecular form
3752
3753        i1fm_at =  kfmat * sqdgat5 * bm2ii(n2n) * (   &
3754                esat25   &
3755             +  two * esat09 * esat04   &
3756             +  esat01 * esat16   &
3757             +  esat64 * esat09   &
3758             +  two * esat36 * esat01   &
3759             +  esat16 * esat01  )
3760
3761       i1_at = ( i1nc_at * i1fm_at ) / ( i1nc_at + i1fm_at  )
3762
3763       coagatat2 = constii * i1_at
3764
3765       qs11 = coagatat2 * bm2iitt(n2n)
3766
3767 ! *** accumulation mode
3768
3769 ! *** near-continuum
3770
3771       i1nc_ac = knc * dgac2 * (   &
3772              two * esac16   &
3773            + esac04 * esac04   &
3774            + esac36 * esac04   &
3775            + a * kngac * (   &
3776                 two * esac04   &
3777            +     esac16 * esac04   &
3778            +     esac36 * esac16 )  )
3779
3780 ! *** free- molecular form
3781
3782        i1fm_ac =  kfmac * sqdgac5 * bm2ii(n2a) * (   &
3783                esac25   &
3784             +  two * esac09 * esac04   &
3785             +  esac01 * esac16   &
3786             +  esac64 * esac09   &
3787             +  two * esac36 * esac01   &
3788             +  esac16 * esac01  )
3789
3790       i1_ac = ( i1nc_ac * i1fm_ac ) / ( i1nc_ac + i1fm_ac  )
3791
3792       coagacac2 = constii * i1_ac
3793
3794       qs22 = coagacac2 * bm2iitt(n2a)
3795
3796
3797       return
3798
3799       end  subroutine getcoags
3800
3801 !----------------------------------------------------------------------
3802 !----------------------------------------------------------------------
3803
3804    end module modal_aero_coag
3805
3806
3807