#define WRF_PORT
#define MODAL_AERO
! module_cam_mam_coag.F
! adapted from cam3 modal_aero_coag.F90 by r.c.easter, august 2010
! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
! 2010-08-04 rce notes
! > change pcnstxx to a subr argument because gas_pcnst is
! not a constant/parameter in wrf-chem
!
!--------------------------------------------------------------
#include "MODAL_AERO_CPP_DEFINES.h"
! modal_aero_coag.F90
!----------------------------------------------------------------------
!BOP
!
! !MODULE: modal_aero_coag --- modal aerosol coagulation
!
! !INTERFACE:
module modal_aero_coag
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use shr_kind_mod, only: r4 => shr_kind_r4
#ifndef WRF_PORT
use chem_mods, only: gas_pcnst
#else
use module_cam_support, only: gas_pcnst => gas_pcnst_modal_aero
#endif
use modal_aero_data, only: maxd_aspectype
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public modal_aero_coag_sub, modal_aero_coag_init
! !PUBLIC DATA MEMBERS:
#ifndef WRF_PORT
integer, parameter :: pcnstxx = gas_pcnst
#endif
#if ( defined MODAL_AERO_7MODE )
integer, parameter, public :: pair_option_acoag = 3
#elif ( defined MODAL_AERO_3MODE )
integer, parameter, public :: pair_option_acoag = 1
#endif
! specifies pairs of modes for which coagulation is calculated
! 1 -- [aitken-->accum]
! 2 -- [aitken-->accum], and [pcarbon-->accum]
! 3 -- [aitken-->accum], [pcarbon-->accum],
! and [aitken-->pcarbon--(aging)-->accum]
! other -- do no coag
integer, parameter, public :: maxpair_acoag = 10
integer, parameter, public :: maxspec_acoag = maxd_aspectype
integer, public :: npair_acoag
integer, public :: modefrm_acoag(maxpair_acoag)
integer, public :: modetoo_acoag(maxpair_acoag)
integer, public :: modetooeff_acoag(maxpair_acoag)
integer, public :: nspecfrm_acoag(maxpair_acoag)
integer, public :: lspecfrm_acoag(maxspec_acoag,maxpair_acoag)
integer, public :: lspectoo_acoag(maxspec_acoag,maxpair_acoag)
! !DESCRIPTION: This module implements ...
!
! !REVISION HISTORY:
!
! RCE 07.04.13: Adapted from MIRAGE2 code
!
!EOP
!----------------------------------------------------------------------
!BOC
! list private module data here
!EOC
!----------------------------------------------------------------------
contains
!----------------------------------------------------------------------
!BOP
! !ROUTINE: modal_aero_coag_sub --- ...
!
! !INTERFACE:
subroutine modal_aero_coag_sub( &
lchnk, ncol, nstep, &
loffset, deltat_main, &
latndx, lonndx, &
t, pmid, pdel, &
q, &
dgncur_a, dgncur_awet, &
wetdens_a &
#ifdef WRF_PORT
, pcnstxx &
#endif
)
!----------------------------------------------------------------------
! Authors: R. Easter
!----------------------------------------------------------------------
! !USES:
#ifndef WRF_PORT
use mo_constants, only: pi
#endif
use modal_aero_data
use modal_aero_gasaerexch, only: n_so4_monolayers_pcage, &
soa_equivso4_factor
#ifndef WRF_PORT
use abortutils, only: endrun
use cam_history, only: outfld, fieldname_len
use chem_mods, only: adv_mass
use constituents, only: pcnst, cnst_name
#else
use module_cam_support, only: endrun, outfld, fieldname_len, pcnst =>pcnst_runtime
use constituents, only: cnst_name
use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
#endif
use physconst, only: gravit, mwdry, r_universal
#ifndef WRF_PORT
use ppgrid, only: pcols, pver
use spmd_utils, only: iam, masterproc
#else
use module_cam_support, only: pcols, pver, iam, masterproc
#endif
implicit none
! !PARAMETERS:
#ifdef WRF_PORT
integer, intent(in) :: pcnstxx ! number of species in "incoming q array"
#endif
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of columns in chunk
integer, intent(in) :: nstep ! model step
integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
integer, intent(in) :: latndx(pcols), lonndx(pcols)
real(r8), intent(in) :: deltat_main ! model timestep (s)
real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
real(r8), intent(inout) :: q(ncol,pver,pcnstxx)
! tracer mixing ratio (TMR) array
! *** MUST BE mol/mol-air or #/mol-air
! *** NOTE ncol & pcnstxx dimensions
real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode)
! dry geo. mean dia. (m) of number distrib.
real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
! wet geo. mean dia. (m) of number distrib.
real(r8), intent(in) :: wetdens_a(pcols,pver,ntot_amode)
! density of wet aerosol (kg/m3)
! !DESCRIPTION:
! computes changes due to coagulation involving
! aitken mode (modeptr_aitken) with accumulation mode (modeptr_accum)
! this version will
! compute changes to mass and number, but not to surface area
! calculates coagulation rate coefficients using either
! new CMAQ V4.6 fast method
! older cmaq slow method (direct gauss-hermite quadrature)
!
! !REVISION HISTORY:
! RCE 07.04.15: Adapted from MIRAGE2 code and CMAQ V4.6 code
!
!EOP
!----------------------------------------------------------------------
!BOC
! local variables
integer :: i, iok, ipair, ip_aitacc, ip_aitpca, ip_pcaacc, iq
integer :: idomode(ntot_amode), iselfcoagdone(ntot_amode)
integer :: jfreqcoag
integer :: k
integer :: l, l1, l2, la, lmz, lsfrm, lstoo, lunout
integer :: modefrm, modetoo, mait, macc, mpca
integer :: n, nfreqcoag
integer, save :: nerr = 0 ! number of errors for entire run
integer, save :: nerrmax = 9999 ! maximum number of errors before abort
integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1
logical, parameter :: fastcoag_flag = .true. ! selects coag rate-coef method
real(r8) :: aircon
real(r8) :: deltat, deltatinv_main
real(r8) :: dr_so4_monolayers_pcage
real(r8) :: dryvol_a(pcols,pver,ntot_amode)
real(r8) :: dumexp, dumloss, dumprod
real(r8) :: dumsfc_frm_old, dumsfc_frm_new
real(r8) :: dum_m2v
real(r8) :: fac_m2v_aitage(maxd_aspectype), fac_m2v_pcarbon(maxd_aspectype)
real(r8) :: fac_volsfc_pcarbon
real(r8) :: lnsg_frm, lnsg_too
real(r8) :: sg_frm, sg_too
real(r8) :: tmpa, tmpb, tmpc, tmpf, tmpg, tmph, tmpn
real(r8) :: tmp1, tmp2
real(r8) :: tmp_qold
real(r8) :: v2ncur_a_tmp
real(r8) :: vol_core, vol_shell
real(r8) :: wetdens_frm, wetdens_too, wetdgnum_frm, wetdgnum_too
real(r8) :: xbetaij0, xbetaij2i, xbetaij2j, xbetaij3, &
xbetaii0, xbetaii2, xbetajj0, xbetajj2
real(r8) :: xferamt, xferfracvol, xferfrac_pcage, xferfrac_max
real(r8) :: xnumbconc(ntot_amode)
real(r8) :: xnumbconcavg(ntot_amode), xnumbconcnew(ntot_amode)
real(r8) :: ybetaij0(maxpair_acoag), ybetaij3(maxpair_acoag)
real(r8) :: ybetaii0(maxpair_acoag), ybetajj0(maxpair_acoag)
real(r8) :: dqdt(ncol,pver,pcnstxx) ! TMR "dq/dt" array - NOTE dims
logical :: dotend(pcnst) ! identifies the species that
! tendencies are computed for
real(r8) :: qsrflx(pcols)
character(len=fieldname_len) :: tmpname
character(len=fieldname_len+3) :: fieldname
! begin
! check if any coagulation pairs exist
if (npair_acoag <= 0) return
!--------------------------------------------------------------------------------
if (ldiag1 > 0) then
if (nstep <= 3) then
do i = 1, ncol
if (lonndx(i) /= 37) cycle
if (latndx(i) /= 23) cycle
if (nstep > 3) cycle
write( *, '(/a,i7,i5,2(2x,2i5))' ) &
'*** modal_aero_coag_sub -- nstep, iam, lat, lon, pcols, ncol =', &
nstep, iam, latndx(i), lonndx(i), pcols, ncol
end do
end if
! if (ncol /= -999888777) return
if (nstep > 3) call endrun( 'modal_aero_coag_sub -- nstep>3 testing halt' )
end if ! (ldiag1 > 0)
!--------------------------------------------------------------------------------
dotend(:) = .false.
dqdt(1:ncol,:,:) = 0.0
lunout = 6
!
! determine if coagulation will be done on this time-step
! currently coagulation is done every 3 hours
!
! deltat = 3600.0*3.0
deltat = deltat_main
nfreqcoag = max( 1, nint( deltat/deltat_main ) )
jfreqcoag = nfreqcoag/2
xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8) ! 1-eps
if (nfreqcoag .gt. 1) then
if ( mod(nstep,nfreqcoag) .ne. jfreqcoag ) return
end if
!
! set idomode
!
idomode(:) = 0
do ipair = 1, npair_acoag
idomode(modefrm_acoag(ipair)) = 1
idomode(modetoo_acoag(ipair)) = 1
end do
!
! other init
!
macc = modeptr_accum
mait = modeptr_aitken
mpca = modeptr_pcarbon
fac_m2v_aitage(:) = 0.0
fac_m2v_pcarbon(:) = 0.0
if (pair_option_acoag == 3) then
! following ipair definitions MUST BE CONSISTENT with
! the coding in modal_aero_coag_init for pair_option_acoag == 3
ip_aitacc = 1
ip_pcaacc = 2
ip_aitpca = 3
! use 1 mol (bi-)sulfate = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10
ipair = ip_aitpca
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair)
if (lsfrm == lptr_so4_a_amode(mait)) then
fac_m2v_aitage(iq) = specmw_so4_amode / specdens_so4_amode
else if (lsfrm == lptr_nh4_a_amode(mait)) then
fac_m2v_aitage(iq) = specmw_nh4_amode / specdens_nh4_amode
else if (lsfrm == lptr_soa_a_amode(mait)) then
fac_m2v_aitage(iq) = soa_equivso4_factor* &
(specmw_soa_amode / specdens_soa_amode)
! for soa, the soa_equivso4_factor converts the soa volume into an
! so4(+nh4) volume that has same hygroscopicity contribution as soa
! this allows aging calculations to be done in terms of the amount
! of (equivalent) so4(+nh4) in the shell
! (see modal_aero_gasaerexch)
end if
end do
do l = 1, nspec_amode(mpca)
l2 = lspectype_amode(l,mpca)
! fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
! [m3-AP/kmol-AP] = [kg-AP/kmol-AP] / [kg-AP/m3-AP]
fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
end do
fac_volsfc_pcarbon = exp( 2.5*(alnsg_amode(mpca)**2) )
else
ip_aitacc = -999888777
ip_pcaacc = -999888777
ip_aitpca = -999888777
end if
!
! loop over levels and columns to calc the coagulation
!
! integrate coagulation changes over deltat = nfreqcoag*deltat_main
! then compute tendencies as
! dqdt = (q(t+deltat) - q(t))/deltat_main
! because tendencies are applied (in physics_update) over deltat_main
!
deltat = nfreqcoag*deltat_main
deltatinv_main = 1.0_r8/(deltat_main*(1.0_r8 + 1.0e-15_r8))
main_k: do k = 1, pver
main_i: do i = 1, ncol
! air molar density (kmol/m3)
aircon = (pmid(i,k)/(r_universal*t(i,k)))
! calculate number conc. (#/m3) for modes doing coagulation
do n = 1, ntot_amode
if (idomode(n) .gt. 0) then
xnumbconc(n) = q(i,k,numptr_amode(n)-loffset)*aircon
xnumbconc(n) = max( 0.0_r8, xnumbconc(n) )
end if
iselfcoagdone(n) = 0
end do
!
! calculate coagulation rates for each pair
!
main_ipair1: do ipair = 1, npair_acoag
modefrm = modefrm_acoag(ipair)
modetoo = modetoo_acoag(ipair)
!
! compute coagulation rates using cmaq "fast" method
! (based on E. Whitby's approximation approach)
! here subr. arguments are all in mks unit
!
call getcoags_wrapper_f( &
t(i,k), pmid(i,k), &
dgncur_awet(i,k,modefrm), dgncur_awet(i,k,modetoo), &
sigmag_amode(modefrm), sigmag_amode(modetoo), &
alnsg_amode(modefrm), alnsg_amode(modetoo), &
wetdens_a(i,k,modefrm), wetdens_a(i,k,modetoo), &
xbetaij0, xbetaij2i, xbetaij2j, xbetaij3, &
xbetaii0, xbetaii2, xbetajj0, xbetajj2 )
! test diagnostics begin --------------------------------------------
if (ldiag2 > 0) then
if (nstep <= 3) then
if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
if ((mod(k-1,5) == 0) .or. (k>=23)) then
wetdgnum_frm = dgncur_awet(i,k,modefrm)
wetdgnum_too = dgncur_awet(i,k,modetoo)
wetdens_frm = wetdens_a(i,k,modefrm)
wetdens_too = wetdens_a(i,k,modetoo)
sg_frm = sigmag_amode(modefrm)
sg_too = sigmag_amode(modetoo)
lnsg_frm = alnsg_amode(modefrm)
lnsg_too = alnsg_amode(modetoo)
call getcoags_wrapper_f( &
t(i,k), pmid(i,k), &
wetdgnum_frm, wetdgnum_too, &
sg_frm, sg_too, &
lnsg_frm, lnsg_too, &
wetdens_frm, wetdens_too, &
xbetaij0, xbetaij2i, xbetaij2j, xbetaij3, &
xbetaii0, xbetaii2, xbetajj0, xbetajj2 )
write(lunout,9801)
write(lunout,9810) 'nstep,lat,lon,k,ipair ', &
nstep, latndx(i), lonndx(i), k, ipair
write(lunout,9820) 'tk, pmb, aircon, pdel ', &
t(i,k), pmid(i,k)*1.0e-2, aircon, pdel(i,k)*1.0e-2
write(lunout,9820) 'wetdens-cgs, sg f/t', &
wetdens_frm*1.0e-3, wetdens_too*1.0e-3, &
sg_frm, sg_too
write(lunout,9820) 'dgnwet-um, dgndry-um f/t', &
1.0e6*wetdgnum_frm, 1.0e6*wetdgnum_too, &
1.0e6*dgncur_a(i,k,modefrm), 1.0e6*dgncur_a(i,k,modetoo)
write(lunout,9820) 'xbeta ij0, ij3, ii0, jj0', &
xbetaij0, xbetaij3, xbetaii0, xbetajj0
write(lunout,9820) 'xbeta ij2i & j, ii2, jj2', &
xbetaij2i, xbetaij2j, xbetaii2, xbetajj2
write(lunout,9820) 'numbii, numbjj, deltat ', &
xnumbconc(modefrm), xnumbconc(modetoo), deltat
write(lunout,9820) 'loss ij3, ii0, jj0 ', &
(xbetaij3*xnumbconc(modetoo)*deltat), &
(xbetaij0*xnumbconc(modetoo)*deltat+ &
xbetaii0*xnumbconc(modefrm)*deltat), &
(xbetajj0*xnumbconc(modetoo)*deltat)
9801 format( / 72x, 'ACOAG' )
9810 format( 'ACOAG ', a, 2i8, 3i7, 3(1pe15.6) )
9820 format( 'ACOAG ', a, 4(1pe15.6) )
9830 format( 'ACOAG ', a, i1, a, 4(1pe15.6) )
end if
end if
end if
end if ! (ldiag2 > 0)
! test diagnostics end ----------------------------------------------
ybetaij0(ipair) = xbetaij0
ybetaij3(ipair) = xbetaij3
ybetaii0(ipair) = xbetaii0
ybetajj0(ipair) = xbetajj0
end do main_ipair1
if ( (pair_option_acoag == 1) .or. &
(pair_option_acoag == 2) ) then
!
! calculate number and mass changes for pair_option_acoag == 1,2
!
main_ipair2: do ipair = 1, npair_acoag
modefrm = modefrm_acoag(ipair)
modetoo = modetoo_acoag(ipair)
! calculate number changes
! apply self-coagulation losses only once to a mode (when iselfcoagdone=0)
! first calc change to "too" mode
! next calc change to "frm" mode, using average number conc of "too"
if ( (mprognum_amode(modetoo) > 0) .and. &
(iselfcoagdone(modetoo) <= 0) ) then
iselfcoagdone(modetoo) = 1
tmpn = xnumbconc(modetoo)
xnumbconcnew(modetoo) = tmpn/(1.0_r8 + deltat*ybetajj0(ipair)*tmpn)
xnumbconcavg(modetoo) = 0.5_r8*(xnumbconcnew(modetoo) + tmpn)
lstoo = numptr_amode(modetoo) - loffset
q(i,k,lstoo) = xnumbconcnew(modetoo)/aircon
dqdt(i,k,lstoo) = (xnumbconcnew(modetoo)-tmpn)*deltatinv_main/aircon
end if
if ( (mprognum_amode(modefrm) > 0) .and. &
(iselfcoagdone(modefrm) <= 0) ) then
iselfcoagdone(modefrm) = 1
tmpn = xnumbconc(modefrm)
tmpa = deltat*ybetaij0(ipair)*xnumbconcavg(modetoo)
tmpb = deltat*ybetaii0(ipair)
tmpc = tmpa + tmpb*tmpn
if (abs(tmpc) < 0.01_r8) then
xnumbconcnew(modefrm) = tmpn*exp(-tmpc)
else if (abs(tmpa) < 0.001_r8) then
xnumbconcnew(modefrm) = &
exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
else
tmpf = tmpb*tmpn/tmpc
tmpg = exp(-tmpa)
tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
xnumbconcnew(modefrm) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
end if
xnumbconcavg(modefrm) = 0.5_r8*(xnumbconcnew(modefrm) + tmpn)
lsfrm = numptr_amode(modefrm) - loffset
q(i,k,lsfrm) = xnumbconcnew(modefrm)/aircon
dqdt(i,k,lsfrm) = (xnumbconcnew(modefrm)-tmpn)*deltatinv_main/aircon
end if
! calculate mass changes
! xbetaij3*xnumbconc(modetoo) = first order loss rate for modefrm volume
! xferfracvol = fraction of modefrm volume transferred to modetoo over deltat
dumloss = ybetaij3(ipair)*xnumbconcavg(modetoo)
xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair) - loffset
lstoo = lspectoo_acoag(iq,ipair) - loffset
if (lsfrm > 0) then
xferamt = q(i,k,lsfrm)*xferfracvol
dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
if (lstoo > 0) then
dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
q(i,k,lstoo) = q(i,k,lstoo) + xferamt
end if
end if
end do
end do main_ipair2
else if (pair_option_acoag == 3) then
!
! calculate number and mass changes for pair_option_acoag == 3
!
! calculate number changes to accum mode
if (mprognum_amode(macc) > 0) then
tmpn = xnumbconc(macc)
xnumbconcnew(macc) = tmpn/(1.0_r8 + deltat*ybetajj0(ip_aitacc)*tmpn)
xnumbconcavg(macc) = 0.5_r8*(xnumbconcnew(macc) + tmpn)
lstoo = numptr_amode(macc) - loffset
q(i,k,lstoo) = xnumbconcnew(macc)/aircon
dqdt(i,k,lstoo) = (xnumbconcnew(macc)-tmpn)*deltatinv_main/aircon
end if
! calculate number changes to primary carbon mode
modefrm = modeptr_pcarbon
if (mprognum_amode(mpca) > 0) then
tmpn = xnumbconc(mpca)
tmpa = deltat*ybetaij0(ip_pcaacc)*xnumbconcavg(macc)
tmpb = deltat*ybetaii0(ip_pcaacc)
tmpc = tmpa + tmpb*tmpn
if (abs(tmpc) < 0.01_r8) then
xnumbconcnew(mpca) = tmpn*exp(-tmpc)
else if (abs(tmpa) < 0.001_r8) then
xnumbconcnew(mpca) = &
exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
else
tmpf = tmpb*tmpn/tmpc
tmpg = exp(-tmpa)
tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
xnumbconcnew(mpca) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
end if
xnumbconcavg(mpca) = 0.5_r8*(xnumbconcnew(mpca) + tmpn)
lsfrm = numptr_amode(mpca) - loffset
q(i,k,lsfrm) = xnumbconcnew(mpca)/aircon
dqdt(i,k,lsfrm) = (xnumbconcnew(mpca)-tmpn)*deltatinv_main/aircon
end if
! calculate number changes to aitken mode
if (mprognum_amode(mait) > 0) then
tmpn = xnumbconc(mait)
tmpa = deltat*( ybetaij0(ip_aitacc)*xnumbconcavg(macc) &
+ ybetaij0(ip_aitpca)*xnumbconcavg(mpca) )
tmpb = deltat*ybetaii0(ip_aitacc)
tmpc = tmpa + tmpb*tmpn
if (abs(tmpc) < 0.01_r8) then
xnumbconcnew(mait) = tmpn*exp(-tmpc)
else if (abs(tmpa) < 0.001_r8) then
xnumbconcnew(mait) = &
exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
else
tmpf = tmpb*tmpn/tmpc
tmpg = exp(-tmpa)
tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
xnumbconcnew(mait) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
end if
xnumbconcavg(mait) = 0.5_r8*(xnumbconcnew(mait) + tmpn)
lsfrm = numptr_amode(mait) - loffset
q(i,k,lsfrm) = xnumbconcnew(mait)/aircon
dqdt(i,k,lsfrm) = (xnumbconcnew(mait)-tmpn)*deltatinv_main/aircon
end if
! calculate mass changes from aitken-->accum direct coagulation and
! aitken-->pcarbon-->accum coagulation/aging
! also calc volume of shell material (so4 & nh4 from aitken-->pcarbon)
dumloss = ybetaij3(ip_aitacc)*xnumbconcavg(macc) &
+ ybetaij3(ip_aitpca)*xnumbconcavg(mpca)
tmpa = ybetaij3(ip_aitpca)*xnumbconcavg(mpca)/max( dumloss, 1.0e-37_r8 )
xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
vol_shell = 0.0
ipair = ip_aitacc
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair) - loffset
lstoo = lspectoo_acoag(iq,ipair) - loffset
if (lsfrm > 0) then
xferamt = q(i,k,lsfrm)*xferfracvol
dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
if (lstoo > 0) then
dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
q(i,k,lstoo) = q(i,k,lstoo) + xferamt
end if
vol_shell = vol_shell + xferamt*tmpa*fac_m2v_aitage(iq)
end if
end do
! now calculate aging transfer fraction for pcarbon-->accum
! this duplicates the code in modal_aero_gasaerexch
vol_core = 0.0
do l = 1, nspec_amode(mpca)
vol_core = vol_core + &
q(i,k,lmassptr_amode(l,mpca)-loffset)*fac_m2v_pcarbon(l)
end do
tmp1 = vol_shell*dgncur_a(i,k,mpca)*fac_volsfc_pcarbon
tmp2 = 6.0_r8*dr_so4_monolayers_pcage*vol_core
tmp2 = max( tmp2, 0.0_r8 )
if (tmp1 >= tmp2) then
xferfrac_pcage = xferfrac_max
else
xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
end if
! calculate mass changes from pcarbon-->accum by direct coagulation
! and aging
dumloss = ybetaij3(ip_pcaacc)*xnumbconcavg(macc)
xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
xferfracvol = xferfracvol + xferfrac_pcage
xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
ipair = ip_pcaacc
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair) - loffset
lstoo = lspectoo_acoag(iq,ipair) - loffset
if (lsfrm > 0) then
xferamt = q(i,k,lsfrm)*xferfracvol
dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
if (lstoo > 0) then
dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
q(i,k,lstoo) = q(i,k,lstoo) + xferamt
end if
end if
end do
lsfrm = numptr_amode(mpca) - loffset
lstoo = numptr_amode(macc) - loffset
if (lsfrm > 0) then
xferamt = q(i,k,lsfrm)*xferfrac_pcage
dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
if (lstoo > 0) then
dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
q(i,k,lstoo) = q(i,k,lstoo) + xferamt
end if
end if
else ! (pair_option_acoag /= 1,2,3) then
write(lunout,*) '*** modal_aero_coag_sub error'
write(lunout,*) ' cannot do _coag_sub error pair_option_acoag =', &
pair_option_acoag
call endrun( 'modal_aero_coag_sub error' )
end if ! (pair_option_acoag == ...)
! test diagnostics begin --------------------------------------------
if (ldiag3 > 0) then
if (nstep <= 3) then
if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
if ((mod(k-1,5) == 0) .or. (k>=23)) then
if (pair_option_acoag == 3) then
write(*,*)
write(lunout,9820) 'xnumbconcavg ait,acc,pca', &
xnumbconcavg(mait), xnumbconcavg(macc), xnumbconcavg(mpca)
write(lunout,9820) 'vshell, core ', &
vol_shell, vol_core
write(lunout,9820) 'dr_mono, dgn ', &
dr_so4_monolayers_pcage, dgncur_a(i,k,mpca)
write(lunout,9820) 'tmp1, tmp2 ', tmp1, tmp2
write(lunout,9820) 'xferfrac_age ', xferfrac_pcage
end if
do ipair = 1, npair_acoag
modefrm = modefrm_acoag(ipair)
modetoo = modetoo_acoag(ipair)
if (npair_acoag > 1) then
write(lunout,*)
write(lunout,9810) 'ipair = ', ipair
end if
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair) - loffset
lstoo = lspectoo_acoag(iq,ipair) - loffset
if (lsfrm > 0) then
tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
! write(lunout,9820) 'm1 frm dqdt/q0,dqdt,q0/1', &
write(lunout,9830) 'm', iq, &
' frm dqdt/q0,dqdt,q0/1', &
dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
end if
if (lstoo > 0) then
tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
write(lunout,9830) 'm', iq, &
' too dqdt/q0,dqdt,q0/1', &
dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
end if
end do ! iq
lsfrm = numptr_amode(modefrm) - loffset
lstoo = numptr_amode(modetoo) - loffset
if (lsfrm > 0) then
tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
write(lunout,9820) 'n frm dqdt/q0,dqdt,q0/1', &
dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
end if
if (lstoo > 0) then
tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
write(lunout,9820) 'n too dqdt/q0,dqdt,q0/1', &
dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
end if
end do ! ipair
end if
end if
end if
end if ! (ldiag3 > 0)
! test diagnostics end ----------------------------------------------
end do main_i
end do main_k
! set dotend's
do ipair = 1, npair_acoag
modefrm = modefrm_acoag(ipair)
modetoo = modetoo_acoag(ipair)
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair) - loffset
lstoo = lspectoo_acoag(iq,ipair) - loffset
if (lsfrm > 0) dotend(lsfrm) = .true.
if (lstoo > 0) dotend(lstoo) = .true.
end do
if (mprognum_amode(modefrm) > 0) then
lsfrm = numptr_amode(modefrm) - loffset
if (lsfrm > 0) dotend(lsfrm) = .true.
end if
if (mprognum_amode(modetoo) > 0) then
lstoo = numptr_amode(modetoo) - loffset
if (lstoo > 0) dotend(lstoo) = .true.
end if
end do
! do history file column-tendency fields
do l = loffset+1, pcnst
lmz = l - loffset
if ( .not. dotend(lmz) ) cycle
qsrflx(:) = 0.0_r8
do k = 1, pver
do i = 1, ncol
qsrflx(i) = qsrflx(i) + dqdt(i,k,lmz)*pdel(i,k)
end do
end do
qsrflx(:) = qsrflx(:)*(adv_mass(lmz)/(gravit*mwdry))
fieldname = trim(cnst_name(l)) // '_sfcoag1'
call outfld( fieldname, qsrflx, pcols, lchnk )
! if (( masterproc ) .and. (nstep < 1)) &
! write(*,'(2(a,2x),1p,e11.3)') &
! 'modal_aero_coag_sub outfld', fieldname, adv_mass(lmz)
end do ! l = ...
return
!EOC
end subroutine modal_aero_coag_sub
!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine modal_aero_coag_init
!
! computes pointers for species transfer during coagulation
!
use modal_aero_data
use modal_aero_gasaerexch, only: &
modefrm_pcage, nspecfrm_pcage, lspecfrm_pcage, lspectoo_pcage
#ifndef WRF_PORT
use abortutils, only: endrun
use cam_history, only: addfld, add_default, fieldname_len, phys_decomp
use constituents, only: pcnst, cnst_name
use spmd_utils, only: masterproc
use phys_control, only: phys_getopts
#else
use module_cam_support, only: endrun, addfld, add_default, fieldname_len, phys_decomp, &
pcnst =>pcnst_runtime, masterproc
use constituents, only: cnst_name
#endif
implicit none
! local variables
integer :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
integer :: l, lsfrm, lstoo, lunout
integer :: m, mfrm, mtoo, mtef
integer :: nsamefrm, nsametoo, nspec
character(len=fieldname_len) :: tmpname
character(len=fieldname_len+3) :: fieldname
character(128) :: long_name
character(8) :: unit
logical :: dotend(pcnst)
logical :: history_aerosol ! Output the MAM aerosol tendencies
!-----------------------------------------------------------------------
#ifndef WRF_PORT
call phys_getopts( history_aerosol_out = history_aerosol )
#else
history_aerosol =.false.
#endif
lunout = 6
!
! define "from mode" and "to mode" for each coagulation pairing
! currently just a2-->a1 coagulation
!
if (pair_option_acoag == 1) then
npair_acoag = 1
modefrm_acoag(1) = modeptr_aitken
modetoo_acoag(1) = modeptr_accum
modetooeff_acoag(1) = modeptr_accum
else if (pair_option_acoag == 2) then
npair_acoag = 2
modefrm_acoag(1) = modeptr_aitken
modetoo_acoag(1) = modeptr_accum
modetooeff_acoag(1) = modeptr_accum
modefrm_acoag(2) = modeptr_pcarbon
modetoo_acoag(2) = modeptr_accum
modetooeff_acoag(2) = modeptr_accum
else if (pair_option_acoag == 3) then
npair_acoag = 3
modefrm_acoag(1) = modeptr_aitken
modetoo_acoag(1) = modeptr_accum
modetooeff_acoag(1) = modeptr_accum
modefrm_acoag(2) = modeptr_pcarbon
modetoo_acoag(2) = modeptr_accum
modetooeff_acoag(2) = modeptr_accum
modefrm_acoag(3) = modeptr_aitken
modetoo_acoag(3) = modeptr_pcarbon
modetooeff_acoag(3) = modeptr_accum
if (modefrm_pcage <= 0) then
write(*,*) '*** modal_aero_coag_init error'
write(*,*) ' pair_option_acoag, modefrm_pcage mismatch'
write(*,*) ' pair_option_acoag, modefrm_pcage =', &
pair_option_acoag, modefrm_pcage
call endrun( 'modal_aero_coag_init error' )
end if
else
npair_acoag = 0
return
end if
!
! define species involved in each coagulation pairing
! (include aerosol water)
!
aa_ipair: do ipair = 1, npair_acoag
mfrm = modefrm_acoag(ipair)
mtoo = modetoo_acoag(ipair)
mtef = modetooeff_acoag(ipair)
if ( (mfrm < 1) .or. (mfrm > ntot_amode) .or. &
(mtoo < 1) .or. (mtoo > ntot_amode) .or. &
(mtef < 1) .or. (mtef > ntot_amode) ) then
write(*,*) '*** modal_aero_coag_init error'
write(*,*) ' ipair, ntot_amode =', ipair, ntot_amode
write(*,*) ' mfrm, mtoo, mtef =', mfrm, mtoo, mtef
call endrun( 'modal_aero_coag_init error' )
end if
mtoo = mtef ! effective modetoo
nspec = 0
aa_iqfrm: do iqfrm = 1, nspec_amode(mfrm)
lsfrm = lmassptr_amode(iqfrm,mfrm)
if ((lsfrm .lt. 1) .or. (lsfrm .gt. pcnst)) cycle aa_iqfrm
! find "too" species having same lspectype_amode as the "frm" species
! several species in a mode may have the same lspectype_amode, so also
! use the ordering as a criterion (e.g., 1st <--> 1st, 2nd <--> 2nd)
iqfrm_aa = 1
iqtoo_aa = 1
if (iqfrm .gt. nspec_amode(mfrm)) then
iqfrm_aa = nspec_amode(mfrm) + 1
iqtoo_aa = nspec_amode(mtoo) + 1
end if
nsamefrm = 0
do iq = iqfrm_aa, iqfrm
if ( lspectype_amode(iq ,mfrm) .eq. &
lspectype_amode(iqfrm,mfrm) ) then
nsamefrm = nsamefrm + 1
end if
end do
nsametoo = 0
lstoo = 0
do iqtoo = iqtoo_aa, nspec_amode(mtoo)
if ( lspectype_amode(iqtoo,mtoo) .eq. &
lspectype_amode(iqfrm,mfrm) ) then
nsametoo = nsametoo + 1
if (nsametoo .eq. nsamefrm) then
lstoo = lmassptr_amode(iqtoo,mtoo)
exit
end if
end if
end do
nspec = nspec + 1
lspecfrm_acoag(nspec,ipair) = lsfrm
lspectoo_acoag(nspec,ipair) = lstoo
end do aa_iqfrm
! lsfrm = lwaterptr_amode(mfrm)
! if ((lsfrm .ge. 1) .and. (lsfrm .le. pcnst)) then
! lstoo = lwaterptr_amode(mtoo)
! if ((lstoo .lt. 1) .or. (lstoo .gt. pcnst)) lstoo = 0
! nspec = nspec + 1
! lspecfrm_acoag(nspec,ipair) = lsfrm
! lspectoo_acoag(nspec,ipair) = lstoo
! end if
nspecfrm_acoag(ipair) = nspec
end do aa_ipair
!
! output results
!
if ( masterproc ) then
write(lunout,9310)
do ipair = 1, npair_acoag
mfrm = modefrm_acoag(ipair)
mtoo = modetoo_acoag(ipair)
mtef = modetooeff_acoag(ipair)
write(lunout,9320) ipair, mfrm, mtoo, mtef
do iq = 1, nspecfrm_acoag(ipair)
lsfrm = lspecfrm_acoag(iq,ipair)
lstoo = lspectoo_acoag(iq,ipair)
if (lstoo .gt. 0) then
write(lunout,9330) lsfrm, cnst_name(lsfrm), &
lstoo, cnst_name(lstoo)
else
write(lunout,9340) lsfrm, cnst_name(lsfrm)
end if
end do
end do ! ipair = ...
write(lunout,*)
end if ! ( masterproc )
9310 format( / 'subr. modal_aero_coag_init' )
9320 format( 'pair', i3, 5x, 'mode', i3, &
' ---> mode', i3, ' eff', i3 )
9330 format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
9340 format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )
!
! create history file column-tendency fields
!
dotend(:) = .false.
do ipair = 1, npair_acoag
do iq = 1, nspecfrm_acoag(ipair)
l = lspecfrm_acoag(iq,ipair)
if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
l = lspectoo_acoag(iq,ipair)
if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
end do
m = modefrm_acoag(ipair)
if ((m > 0) .and. (m <= ntot_amode)) then
l = numptr_amode(m)
if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
end if
m = modetoo_acoag(ipair)
if ((m > 0) .and. (m <= ntot_amode)) then
l = numptr_amode(m)
if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
end if
end do ! ipair = ...
if (pair_option_acoag == 3) then
do iq = 1, nspecfrm_pcage
lsfrm = lspecfrm_pcage(iq)
lstoo = lspectoo_pcage(iq)
if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
dotend(lsfrm) = .true.
if ((lstoo > 0) .and. (lstoo <= pcnst)) then
dotend(lstoo) = .true.
end if
end if
end do
end if
do l = 1, pcnst
if ( .not. dotend(l) ) cycle
tmpname = cnst_name(l)
unit = 'kg/m2/s'
do m = 1, ntot_amode
if (l == numptr_amode(m)) unit = '#/m2/s'
end do
fieldname = trim(tmpname) // '_sfcoag1'
long_name = trim(tmpname) // ' modal_aero coagulation column tendency'
call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
if ( history_aerosol ) then
call add_default( fieldname, 1, ' ' )
endif
if ( masterproc ) write(*,'(3(a,2x))') &
'modal_aero_coag_init addfld', fieldname, unit
end do ! l = ...
return
end subroutine modal_aero_coag_init
!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine calc_coag_coef4( &
dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj, &
wetddrydi, wetddrydj, &
temp, presscgs, lunerr, lunout, iok, &
betaij0, betaij2i, betaij2j, betaij3, &
betaii0, betaii2, betajj0, betajj2 )
!
! computes the following coagulation rate "coefficients"
! for lognormal modes
!
! self coagulation
! dNi/dt = - betaii0*Ni*Ni
! dNj/dt = - betajj0*Nj*Nj
! dSi/dt = - betaii2*Si*Ni
! dSj/dt = - betajj2*Sj*Nj
!
! modei-modej coagulation
! dNi/dt = - betaij0*Ni*Nj
! dNj/dt = 0.
! dSi/dt = - betaij2i*Si*Nj
! dSj/dt = + betaij2j*Si*Nj
! dVi/dt = - betaij3*Vi*Nj == - dVj/dt
!
! Ni, Nj are number concentrations (particles/cm3)
! Si, Sj are surface concentrations (cm2/cm3)
! Vi, Vj are volume concentrations (cm3/cm3)
! rates are (N,S,V)/s
!
! mode i is the smaller mode (e.g., Aitken)
! mode j is the larger mode (e.g., accumulation)
!
! input arguments
! dgni, dgnj = DRY median diameters for number distribution (cm)
! alnsgi, alnsgj = ln of geometric standard deviations (dimensionless)
! rhopi, rhopj = WET density (g/cm3)
! wetddrydi, wetddrydj = (WET median diameter)/(DRY median diameter)
! temp = air temperature (K)
! presscgs = air pressure (dynes/cm2)
! lunerr, lunout = logical unit for error or diagnostic output
!
! output arguments
! iok = status flag (+1 = success, 0/negative = failure)
! beta--- = see above
!
implicit none
! arguments
integer, intent(in) :: lunerr, lunout
integer, intent(out) :: iok
real(r4), intent(in) :: &
dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj, &
wetddrydi, wetddrydj, &
temp, presscgs
real(r4), intent(out) :: &
betaij0, betaij2i, betaij2j, betaij3, &
betaii0, betaii2, betajj0, betajj2
! local variables
real(r4) airprs, airtemp
real(r4) dgacc, dgatk, pdensac, pdensat
real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij
real(r4) dp2bar_mks_i, dp2bar_mks_j, dp3bar_mks_i
! check for reasonable inputs
iok = -1
if ((dgni .lt. 1.e-8) .or. (dgni .gt. 1.)) then
write(lunerr,9100) 'dgni', dgni
return
else if ((dgnj .lt. 1.e-8) .or. (dgnj .gt. 1.)) then
write(lunerr,9100) 'dgnj', dgnj
return
! else if (dgni .gt. dgnj) then
! write(lunerr,9110) dgni, dgnj
! return
else if ((alnsgi .lt. 0.0) .or. (alnsgi .gt. 2.3)) then
write(lunerr,9100) 'alnsgi', alnsgi
return
else if ((alnsgj .lt. 0.0) .or. (alnsgj .gt. 2.3)) then
write(lunerr,9100) 'alnsgj', alnsgj
return
else if ((rhopi .lt. 0.01) .or. (rhopi .gt. 100.)) then
write(lunerr,9100) 'rhopi', rhopi
return
else if ((rhopj .lt. 0.01) .or. (rhopj .gt. 100.)) then
write(lunerr,9100) 'rhopj', rhopj
return
else if ((temp .lt. 10.) .or. (temp .gt. 370.)) then
write(lunerr,9100) 'temp', temp
return
else if ((presscgs .lt. 1.e2) .or. (presscgs .gt. 1.5e6)) then
write(lunerr,9100) 'presscgs', presscgs
return
end if
9100 format( '*** subr. calc_coag_coef_4 - bad value for ', a, &
' = ', 1pe15.5 )
!9110 format( '*** subr. calc_coag_coef_4 - dgni > dgnj -- ',
! + 2(1pe15.5) )
iok = +1
!
! define mks variables
!
airtemp = temp
! dyne/cm2 to pa
airprs = presscgs * 1.0e-1
! cm to m
dgatk = dgni * 1.0e-2
dgacc = dgnj * 1.0e-2
! g/cm3 to kg/m3
pdensat = rhopi * 1.0e+3
pdensac = rhopj * 1.0e+3
!
! call interace to binkowski routines
!
call bink_coag_rates( &
dgatk, dgacc, alnsgi, alnsgj, pdensat, pdensac, &
wetddrydi, wetddrydj, &
airtemp, airprs, lunerr, iok, &
batat, batac, bacat, bacac, c3ij )
!
! transfer the batat, batac, bacat, bacac values
! to the beta--- variables
!
! self coagulation, number
! dNi/dt = - betaii0*Ni*Ni
! dNXi/dt = - batat(1)*NXi*NXi
! where Ni is (particles/cm3), NXi is (particles/m3)
! the first-order loss rates are just (1/s) and thus are equal
! d[ln(Ni)]/dt = d[ln(NXi)]/dt
! so
! betaii0 = batat(1) * (NXi/Ni)
! conversion factors are
! 1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
! (NXi/Ni) = 1.0e6
!
betaii0 = batat(1) * 1.0e-14
betajj0 = bacac(1) * 1.0e-14
! self coagulation, surface
! dSi/dt = - betaii2*Si*Ni
! dM2Xi/dt = - batat(1)*NXi*NXi
! where Si is surface in (cm2/cm3) and M2Xi is (surface/pi) in (m2/m3)
! the first-order loss rates are just (1/s) and thus are equal
! betaii2 = batat(2) * (NXi/M2Xi) * (NXi/Ni)
! conversion factors are
! 1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
! (NXi/Ni) = 1.0e6
! (M2Xi/NXi) = dp2bar_mks_i/pi, where dp2bar_mks_i is the average
! Dp**2 in m**2
!
dp2bar_mks_j = (dgnj**2) * exp( 2.0*alnsgj*alnsgj ) * 1.0e-4
dp2bar_mks_i = (dgni**2) * exp( 2.0*alnsgi*alnsgi ) * 1.0e-4
dp3bar_mks_i = (dgni**3) * exp( 4.5*alnsgi*alnsgi ) * 1.0e-6
betaii2 = (batat(2) / dp2bar_mks_i) * 1.0e-14
betajj2 = (bacac(2) / dp2bar_mks_j) * 1.0e-14
! modei-modej coagulation, number
! dNi/dt = - betaij0*Ni*Nj
! dNj/dt = 0.
! conversions as for self coagulation
!
betaij0 = batac(1) * 1.0e-14
! modei-modej coagulation, surface
! dSi/dt = - betaij2i*Si*Nj
! dSj/dt = + betaij2j*Si*Nj
! conversions as for self coagulation
!
betaij2i = (batac(2) / dp2bar_mks_i) * 1.0e-14
betaij2j = (bacat(2) / dp2bar_mks_i) * 1.0e-14
! modei-modej coagulation, volume
! dVi/dt = - betaij3*Vi*Nj
! dM3Xi/dt = - c3ij*NXi*NXj
! where Vi is (cm3/cm3) and M3Xi is (volume*6/pi) in (m3/m3)
! the first-order loss rates are just (1/s) and thus are equal
! betaij3 = c3ij * (NXi/M3Xi) * (NXj/Nj)
! conversion factors are
! 1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
! (NXj/Nj) = 1.0e6
! (M3Xi/NXi) = dp3bar_mks_i*(6/pi), where dp3bar_mks_i is the average
! Dp**3 in m**3
!
betaij3 = ( c3ij / dp3bar_mks_i ) * 1.0e-14
return
end subroutine calc_coag_coef4
!-----------------------------------------------------------------------
subroutine bink_coag_rates( &
dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac, &
wetddrydat, wetddrydac, &
airtemp, airprs, lunerr, iok, &
batat, batac, bacat, bacac, c3ij )
!
! provides interface to F. Binkowski's intracoag and intercoag
!
! this code was "cut" from F. Binkowski's aero_info_ae3.f & aero_subs3_ae3.f
!
! computes the following coagulation rate "coefficients"
!
! self coagulation
! dNXi/dt = - (1.e-20*batat(1))*NXi*NXi
! dNXj/dt = - (1.e-20*bacac(1))*NXj*NXj
! dM2Xi/dt = - (1.e-20*batat(2))*NXi*NXi
! dM2Xj/dt = - (1.e-20*bacac(2))*NXj*NXj
!
! modei-modej coagulation
! dNXi/dt = - (1.e-20*batac(1))*NXi*NXj
! dNXj/dt = 0.
! dM2Xi/dt = - (1.e-20*batac(2))*NXi*NXj
! dM2Xj/dt = + (1.e-20*bacat(2))*NXi*NXj
! dM3Xi/dt = - (1.e-20*c3ij)*NXi*NXj == - dM3Xj/dt
!
! NXi, NXj are number concentrations (particles/m3)
! M2Xi, M2Xj are 2nd moment (surface/pi) concentrations (m2/m3)
! M3Xi, M3Xj are 2rd moment (volume*6/pi) concentrations (m3/m3)
!
! in the above, mode i is aitken mode, mode j is accumulation mode
!
! *** note that the batat, batac, bacat, bacac, c3ij
! all must be multiplied by 1.e-20 to undo the 1.e+20 scaling
! in the real*4 versions of intracoag and intercoag
!
! input arguments
! dgatk, dgacc = DRY median diameters for number distribution (m)
! for aitken (atk or at) and accumulation (acc or ac) modes
! xxlsgat, xxlsgac = ln of geometric standard deviations (dimensionless)
! pdensat, pdensac = WET density (kg/m3)
! wetddrydat, wetddrydac = (WET median diameter)/(DRY median diameter)
! airtemp = air temperature (K)
! airprs = air pressure (Pa)
! lunerr = logical unit for error output
!
! output arguments
! iok = status flag (+1 = success, 0/negative = failure)
! batat, batac, bacat, bacac, c3ij = see above
!
implicit none
! arguments
integer lunerr, iok
real(r4) dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac, &
wetddrydat, wetddrydac, &
airtemp, airprs
real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij
! *** modal geometric mean diameters: [ m ]
! real dgatk ! nuclei mode
! real dgacc ! accumulation mode
! *** log of modal geometric standard deviation
! real xxlsgat ! aitken mode
! real xxlsgac ! accumulation mode
! *** average modal particle densities [ kg/m**3 ]
! real pdensat ! nuclei mode
! real pdensac ! accumulation mode
! real airtemp ! air temperature [ k ]
! real airprs ! air pressure in [ pa ]
! local variables
real(r4) two3
parameter( two3 = 2.0/3.0 )
real(r4) avo ! avogadro's constant [ 1/mol ]
parameter ( avo = 6.0221367e23 )
real(r4) rgasuniv ! universal gas constant [ j/mol-k ]
parameter ( rgasuniv = 8.314510 )
real(r4) boltz ! boltzmann's constant [ j / k]
parameter ( boltz = rgasuniv / avo )
real(r4) p0 ! starting standard surface pressure [ pa ]
parameter ( p0 = 101325.0 )
real(r4) t0 ! starting standard surface temperature [ k ]
parameter ( t0 = 288.15 )
real(r4) xlm ! atmospheric mean free path [ m ]
real(r4) amu ! atmospheric dynamic viscosity [ kg/m s ]
real(r4) kfm, knc, lamda, sqrt_temp
! fsb calculate the square root of the ambient
! temperature for later use
! *** calculate mean free path [ m ]:
! *** 6.6328e-8 is the sea level values given in table i.2.8
! *** on page 10 of u.s. standard atmosphere 1962
xlm = 6.6328e-8 * p0 * airtemp / ( t0 * airprs )
! *** calculate dynamic viscosity [ kg m**-1 s**-1 ]:
! *** u.s. standard atmosphere 1962 page 14 expression
! for dynamic viscosity is:
! dynamic viscosity = beta * t * sqrt(t) / ( t + s)
! where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
sqrt_temp = sqrt( airtemp)
amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )
! *** coagulation
! *** set up coagulation rates
! *** moment independent factors
knc = two3 * boltz * airtemp / amu
lamda = xlm
! *** calculate the coagulation coefficients for use
! in the aitken (nuclei) & accumulation modes
! *** with gauss-hermite numerical quadrature
! using 10 abscissas.
! *** aitken - aitken mode coagulation
kfm = sqrt( 3.0 * boltz * airtemp / pdensat )
call intracoag_gh(lamda, &
kfm, knc, &
dgatk, &
xxlsgat, &
wetddrydat, &
batat(2), &
batat(1) )
! *** accumulation - accumulation mode coagulation
kfm = sqrt( 3.0 * boltz * airtemp / pdensac )
call intracoag_gh(lamda, &
kfm, knc, &
dgacc, &
xxlsgac, &
wetddrydac, &
bacac(2), &
bacac(1) )
! *** aitken accumulation mode coagulation
bacat(1) = 0.0 ! not used
kfm = sqrt( 6.0 * boltz * airtemp / &
( pdensat + pdensac ) )
call intercoag_gh(lamda, &
kfm, knc, &
dgatk, dgacc , &
xxlsgat, xxlsgac , &
wetddrydat, wetddrydac, &
batac(2), &
bacat(2), &
batac(1), &
c3ij )
! c30atac = c3ij * cblk( vac0 ) * cblk( vat0 )
! loss = c30atac / cblk( vat3 )
return
end subroutine bink_coag_rates
! ---------------------------------------------------------------------
! fsb subrs to do gauss-hermite numerical quadrature
subroutine intracoag_gh( lamda, kfm, knc, &
dg, xlnsig, wetddryd, &
quads11, quadn11)
! fsb this version is for intramodal coagulation for number
! and second moment
! fsb this version runs in real*4 arithmetic
! *** this version calculates the coagulation coefficients
! using the harmonic mean approach for both fm and nc cases.
! *** does gauss-hermite quadrature for intra-modal
! coagulation integrals for 2nd moment
! for a lognormal distribution
! defined by dg,xlnsig,
! *** dg and xlnsig are the geometric mean diameters (meters)
! and the logarithms of the
! geometric standard deviations (dimensionless)
! whose meaning is defined below at the end of the routine
! ghxi, ghwi are the gauss-hermite weights and n is one-half the
! number of abscissas, since an even number of abscissas is used
!
!.......................................................................
! (following comments added by rc_easter)
!
! computes the following coagulation rate "coefficients"
! for intramodal coagulation
! dNX/dt = - (1.e-20*quadn11)*NX1*NX1
! dM2X/dt = - (1.e-20*quads11)*NX1*NX1
!
! NX is the mode's number concentration (particles/m3)
! M2X is the mode's 2nd moment (surface/pi) concentration (m2/m3)
!
! input arguments
! lamda = mean free path (m)
! kfm, knc = constants used in free-molecular and near-continuum
! calculations (see subr bink_coag_rates)
! dg = DRY median diameter for number distribution (m)
! xlnsig = ln of geometric standard deviation (dimensionless)
! wetddryd = (WET median diameter)/(DRY median diameter)
!
! output arguments
! quads11, quadn11 = coagulation rate coefficients (see above)
!
! *** note that the quads11, quadn11
! are all scaled by 1.e+20 to avoid underflow, and they
! must be multiplied by 1.e-20 when they are applied
!
! *** wetddryd1 was added because MIRAGE treats N/S/V of
! the DRY size distribution, so the quadrature should be done
! over the DRY size distribution. However, the coagulation kernel
! must be computed using actual (WET) particle sizes
!
!.......................................................................
!
implicit none
integer i,j
real(r4) lamda ! mean free path
real(r4) kfm, knc
real(r4) dg, xlnsig, wetddryd
real(r4) quads11, quadn11
real(r4) pi
parameter( pi = 3.14159265358979)
real(r4) two3rds
parameter( two3rds = 2.0d0 / 3.0d0 )
real(r4) sqrt2
parameter(sqrt2 = 1.41421356237309 )
real(r4) sum1sfm, sum2sfm, sum1nfm, sum2nfm
real(r4) sum1snc, sum2snc, sum1nnc, sum2nnc
real(r4) xi, wxi, xf, dp1p,dp1m,dp1psq,dp1msq, dp1pwet, dp1mwet
real(r4) v1p,v1m, a2p,a2m,v2p,v2m
real(r4) yi,wyi,yf,dp2p,dp2m,dp2psq,dp2msq, dp2pwet, dp2mwet
real(r4) dspp,dsmp,dspm, dsmm
real(r4) bppfm,bmpfm,bpmfm,bmmfm
real(r4) bppnc,bmpnc,bpmnc,bmmnc
real(r4) xx1, xx2
real(r4) xbsfm, xbsnc, xbnfm, xbnnc
real(r4) betafm, betanc
real(r4) a ! approx cunningham corr. factor
parameter( a = 1.246d0 )
real(r4) twoa
parameter( twoa = 2.0d0 * a )
! *** has a fixed number of gauss-herimite abscissas (n)
integer n ! one-half the number of abscissas
parameter ( n = 5 )
real(r4), save :: ghxi(n) ! gauss-hermite abscissas
real(r4), save :: ghwi(n) ! gauss-hermite weights
! ** values from table 25.10 (page 924) of abramowitz and stegun,
! handbook of mathematical functions, national bureau of standards,
! december 1965.
! breaks in number to facilitate comparison with printed table
! *** tests show that 10 point is adquate.
data ghxi/0.342901327223705, &
1.036610829789514, &
1.756683649299882, &
2.532731674232790, &
3.436159118837738/
data ghwi/6.108626337353d-1, &
2.401386110823d-1, &
3.387439445548d-2, &
1.343645746781d-3, &
7.640432855233d-6/
! *** the following expressions are from binkowski & shanker
! jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
! december 20, 1995
! *** for free molecular eq. a5
betafm(xx1, xx2) = kfm * &
sqrt(1.0 / xx1**3 + 1.0 / xx2**3 ) * (xx1 + xx2)**2
! *** for near continuum eq. a6
betanc(xx1, xx2) = knc * (xx1 + xx2) * &
( 1.0d0 / xx1 + 1.0d0 / xx2 + &
twoa * lamda * ( 1.0d0 / xx1 ** 2 &
+ 1.0d0 / xx2 **2 ) )
sum1sfm = 0.0
sum1snc = 0.0
sum1nfm = 0.0
sum1nnc = 0.0
do 1 i=1,n
sum2sfm = 0.0
sum2snc = 0.0
sum2nfm = 0.0
sum2nnc = 0.0
xi = ghxi(i)
wxi = ghwi(i)
xf = exp( sqrt2 * xi *xlnsig)
dp1p = dg*xf
dp1m = dg/xf
dp1psq = dp1p*dp1p
dp1msq = dp1m*dp1m
v1p = dp1p*dp1psq
v1m = dp1m*dp1msq
dp1pwet = dp1p * wetddryd
dp1mwet = dp1m * wetddryd
do 11 j=1,n
yi = ghxi(j)
wyi = ghwi(j)
yf = exp( sqrt2 * yi * xlnsig)
dp2p = dg*yf
dp2m = dg/yf
dp2psq = dp2p*dp2p
dp2msq = dp2m*dp2m
a2p = dp2psq
a2m = dp2msq
v2p = dp2p*dp2psq
v2m =dp2m*dp2msq
dspp = 0.5*(v1p+v2p)**two3rds - a2p
dsmp = 0.5*(v1m+v2p)**two3rds - a2p
dspm = 0.5*(v1p+v2m)**two3rds - a2m
dsmm = 0.5*(v1m+v2m)**two3rds - a2m
dp2pwet = dp2p * wetddryd
dp2mwet = dp2m * wetddryd
! scale by 1.0e+20 to avoid underflow
bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20
bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20
sum2sfm = sum2sfm + wyi*(dspp * bppfm + dspm * bpmfm &
+ dsmp * bmpfm + dsmm * bmmfm )
sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)
sum2snc = sum2snc + wyi*(dspp * bppnc + dspm * bpmnc &
+ dsmp * bmpnc + dsmm * bmmnc )
sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)
11 continue
sum1sfm = sum1sfm + wxi * sum2sfm
sum1nfm = sum1nfm + wxi * sum2nfm
sum1snc = sum1snc + wxi * sum2snc
sum1nnc = sum1nnc + wxi * sum2nnc
1 continue
xbsfm = -sum1sfm / pi
xbsnc = -sum1snc / pi
! quads11 = xbsfm * xbsnc / ( xbsfm + xbsnc )
quads11 = ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) ) &
* min(xbsfm,xbsnc)
! *** quads11 is the intra-modal coagulation term for 2nd moment
xbnfm = 0.5 * sum1nfm / pi
xbnnc = 0.5 * sum1nnc / pi
! quadn11 = xbnfm * xbnnc / ( xbnfm + xbnnc )
quadn11 = ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) ) &
* min(xbnfm,xbnnc)
! *** quadn11 is the intra-modal coagulation term for number
return
end subroutine intracoag_gh
! ---------------------------------------------------------------------
subroutine intercoag_gh( lamda, kfm, knc, dg1, dg2, &
xlnsig1, xlnsig2, &
wetddryd1, wetddryd2, &
quads12, quads21, quadn12, quadv12 )
! fsb this version is for intermodal coagulation for number,
! second, and third moments
! fsb this version runs in real*4 arithmetic
! *** this version calculates the coagulation coefficients
! using the harmonic mean approach for both fm and nc cases.
! *** does gauss-hermite quadrature for inter-modal
! coagulation integrals for 2nd moment
! for two lognormal distributions
! defined by dg1,xlnsig1, dg2,xlnsig2
! *** dg and xlnsig are the geometric mean diameters (meters)
! and the logarithms of the
! geometric standard deviations (dimensionless)
! whose meaning is defined below at the end of the routine
! ghxi, ghwi are the gauss-hermite weights and n is one-half the
! number of abscissas, since an even number of abscissas is used
!
!.......................................................................
! (following comments added by rc_easter)
!
! computes the following coagulation rate "coefficients"
! for mode1-mode2 coagulation
! dNX1/dt = - (1.e-20*quadn12)*NX1*NX2
! dNX2/dt = 0.
! dM2X1/dt = - (1.e-20*quads12)*NX1*NX2
! dM2X2/dt = + (1.e-20*quads21)*NX1*NX2
! dM3X1/dt = - (1.e-20*quads12)*NX1*NX2 == - dM3X2/dt
!
! NX1, NX2 are number concentrations (particles/m3)
! M2X1, M2X2 are 2nd moment (surface/pi) concentrations (m2/m3)
! M3X1, M3X2 are 2rd moment (volume*6/pi) concentrations (m3/m3)
!
! in the above, mode 1 is aitken mode, mode 2 is accumulation mode
!
! input arguments
! lamda = mean free path (m)
! kfm, knc = constants used in free-molecular and near-continuum
! calculations (see subr bink_coag_rates)
! dg1, dg2 = DRY median diameters for number distribution (m)
! for aitken and accumulation modes
! xlnsig1, xlnsig2 = ln of geometric standard deviations (dimensionless)
! wetddryd1, wetddryd2 = (WET median diameter)/(DRY median diameter)
!
! output arguments
! quads12, quads21, quadn12, quadv12 = coagulation rate
! coefficients (see above)
!
! *** note that the quads12, quads21, quadn12, quadv12
! are all scaled by 1.e+20 to avoid underflow, and they
! must be multiplied by 1.e-20 when they are applied
!
! *** wetddryd1/2 were added because MIRAGE treats N/S/V of
! the DRY size distribution, so the quadrature should be done
! over the DRY size distribution. However, the coagulation kernel
! must be computed using actual (WET) particle sizes
!
!.......................................................................
implicit none
integer i,j
real(r4) lamda ! mean free path
real(r4) kfm, knc
real(r4) dg1, xlnsig1, dg2, xlnsig2
real(r4) wetddryd1, wetddryd2
real(r4) quads12, quads21, quadn12, quadv12
real(r4) pi
parameter( pi = 3.14159265358979)
real(r4) two3rds
parameter( two3rds = 2.0d0 / 3.0d0 )
real(r4) sqrt2
parameter(sqrt2 = 1.41421356237309 )
real(r4) sum1s12fm, sum1s21fm, sum2s12fm, sum2s21fm
real(r4) sum1nfm, sum2nfm
real(r4) sum1vfm, sum2vfm
real(r4) sum1s12nc, sum1s21nc, sum2s12nc, sum2s21nc
real(r4) sum1nnc, sum2nnc
real(r4) sum1vnc, sum2vnc
real(r4) xi, wxi,xf, dp1p, dp1m, dp1psq, dp1msq, dp1pwet, dp1mwet
real(r4) a1p, a1m, v1p, v1m
real(r4) a2p, a2m, v2p, v2m
real(r4) yi, wyi, yf, dp2p, dp2m, dp2psq, dp2msq, dp2pwet, dp2mwet
real(r4) dspp, dsmp, dspm, dsmm
real(r4) bppfm, bmpfm, bpmfm, bmmfm
real(r4) bppnc, bmpnc, bpmnc, bmmnc
real(r4) xx1, xx2
real(r4) xbsfm, xbsnc, xbnfm, xbnnc, xbvfm, xbvnc
real(r4) betafm, betanc
real(r4) a ! approx cunningham corr. factor
parameter( a = 1.246d0 )
real(r4) twoa
parameter( twoa = 2.0d0 * a )
! *** has a fixed number of gauss-herimite abscissas (n)
integer n ! one-half the number of abscissas
parameter ( n = 5 )
real(r4), save :: ghxi(n) ! gauss-hermite abscissas
real(r4), save :: ghwi(n) ! gauss-hermite weights
! ** values from table 25.10 (page 924) of abramowitz and stegun,
! handbook of mathematical functions, national bureau of standards,
! december 1965.
! breaks in number to facilitate comparison with printed table
! *** tests show that 10 point is adquate.
data ghxi/0.342901327223705, &
1.036610829789514, &
1.756683649299882, &
2.532731674232790, &
3.436159118837738/
data ghwi/6.108626337353d-1, &
2.401386110823d-1, &
3.387439445548d-2, &
1.343645746781d-3, &
7.640432855233d-6/
! *** the following expressions are from binkowski & shanker
! jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
! december 20, 1995
! *** for free molecular eq. a5
betafm(xx1, xx2) = kfm * &
sqrt(1.0 / xx1**3 + 1.0 / xx2**3 ) * (xx1 + xx2)**2
! *** for near continuum eq. a6
betanc(xx1, xx2) = knc * (xx1 + xx2) * &
( 1.0d0 / xx1 + 1.0d0 / xx2 + &
twoa * lamda * ( 1.0d0 / xx1 ** 2 &
+ 1.0d0 / xx2 **2 ) )
sum1s12fm = 0.0
sum1s12nc = 0.0
sum1s21fm = 0.0
sum1s21nc = 0.0
sum1vnc = 0.0
sum1vfm = 0.0
sum1nfm = 0.0
sum1nnc = 0.0
do 1 i=1,n
sum2s12fm = 0.0
sum2s12nc = 0.0
sum2s21fm = 0.0
sum2s21nc = 0.0
sum2nfm = 0.0
sum2nnc = 0.0
sum2vnc = 0.0
sum2vfm = 0.0
xi = ghxi(i)
wxi = ghwi(i)
xf = exp( sqrt2 * xi *xlnsig1)
dp1p = dg1*xf
dp1m = dg1/xf
dp1psq = dp1p*dp1p
dp1msq = dp1m*dp1m
a1p = dp1psq
a1m = dp1msq
v1p = dp1p*dp1psq
v1m = dp1m*dp1msq
dp1pwet = dp1p * wetddryd1
dp1mwet = dp1m * wetddryd1
do 11 j=1,n
yi = ghxi(j)
wyi = ghwi(j)
yf = exp( sqrt2 * yi * xlnsig2)
dp2p = dg2*yf
dp2m = dg2/yf
dp2psq = dp2p*dp2p
dp2msq = dp2m*dp2m
a2p = dp2psq
a2m = dp2msq
v2p = dp2p*dp2psq
v2m = dp2m*dp2msq
dspp = (v1p+v2p)**two3rds - a2p
dsmp = (v1m+v2p)**two3rds - a2p
dspm = (v1p+v2m)**two3rds - a2m
dsmm = (v1m+v2m)**two3rds - a2m
dp2pwet = dp2p * wetddryd2
dp2mwet = dp2m * wetddryd2
! scale by 1.0e+20 to avoid underflow
bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20
bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20
sum2s12fm = sum2s12fm + wyi*(a1p * bppfm + a1p * bpmfm &
+ a1m * bmpfm + a1m * bmmfm )
sum2s21fm = sum2s21fm + wyi*(dspp * bppfm + dspm * bpmfm &
+ dsmp * bmpfm + dsmm * bmmfm )
sum2s12nc = sum2s12nc + wyi*(a1p * bppnc + a1p * bpmnc &
+ a1m * bmpnc + a1m * bmmnc )
sum2s21nc = sum2s21nc + wyi*(dspp * bppnc + dspm * bpmnc &
+ dsmp * bmpnc + dsmm * bmmnc )
sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)
sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)
sum2vfm = sum2vfm + wyi*(v1p*(bppfm + bpmfm) + &
v1m*(bmpfm + bmmfm) )
sum2vnc = sum2vnc + wyi*(v1p*(bppnc + bpmnc) + &
v1m*(bmpnc + bmmnc) )
11 continue
sum1s12fm = sum1s12fm + wxi * sum2s12fm
sum1s21fm = sum1s21fm + wxi * sum2s21fm
sum1nfm = sum1nfm + wxi * sum2nfm
sum1vfm = sum1vfm + wxi * sum2vfm
sum1s12nc = sum1s12nc + wxi * sum2s12nc
sum1s21nc = sum1s21nc + wxi * sum2s21nc
sum1nnc = sum1nnc + wxi * sum2nnc
sum1vnc = sum1vnc + wxi * sum2vnc
1 continue
! *** second moment intermodal coagulation coefficients
! fsb note: the transfer of second moment is not symmetric.
! see equations a3 & a4 of binkowski & shankar (1995)
! *** to accumulation mode from aitken mode
xbsfm = sum1s21fm / pi
xbsnc = sum1s21nc / pi
! quads21 = xbsfm * xbsnc / ( xbsfm + xbsnc )
quads21 = ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) ) &
* min(xbsfm,xbsnc)
! *** from aitken mode to accumulation mode
xbsfm = sum1s12fm / pi
xbsnc = sum1s12nc / pi
! quads12 = xbsfm * xbsnc / ( xbsfm + xbsnc )
quads12 = ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) ) &
* min(xbsfm,xbsnc)
xbnfm = sum1nfm / pi
xbnnc = sum1nnc / pi
! quadn12 = xbnfm * xbnnc / ( xbnfm + xbnnc )
quadn12 = ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) ) &
* min(xbnfm,xbnnc)
! *** quadn12 is the intermodal coagulation coefficient for number
xbvfm = sum1vfm / pi
xbvnc = sum1vnc / pi
! quadv12 = xbvfm * xbvnc / ( xbvfm + xbvnc )
quadv12 = ( max(xbvfm,xbvnc) / ( xbvfm + xbvnc ) ) &
* min(xbvfm,xbvnc)
! *** quadv12 is the intermodal coagulation coefficient for 3rd moment
return
end subroutine intercoag_gh
!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine getcoags_wrapper_f( &
airtemp, airprs, &
dgatk, dgacc, &
sgatk, sgacc, &
xxlsgat, xxlsgac, &
pdensat, pdensac, &
betaij0, betaij2i, betaij2j, betaij3, &
betaii0, betaii2, betajj0, betajj2 )
!
! interface to subr. getcoags
!
! interface code adapted from subr. aeroproc of cmaq v4.6,
! with some of the parameter values from module aero_info_ae4
!
implicit none
! *** arguments
real(r8), intent(in) :: airtemp ! air temperature [ k ]
real(r8), intent(in) :: airprs ! air pressure in [ pa ]
real(r8), intent(in) :: dgatk ! aitken mode geometric mean diameter [m]
real(r8), intent(in) :: dgacc ! accumulation mode geometric mean diam [m]
real(r8), intent(in) :: sgatk ! aitken mode geometric standard deviation
real(r8), intent(in) :: sgacc ! accumulation mode geometric standard deviation
real(r8), intent(in) :: xxlsgat ! natural log of geometric standard
real(r8), intent(in) :: xxlsgac ! deviations
real(r8), intent(in) :: pdensat ! aitken mode particle density [ kg / m**3 ]
real(r8), intent(in) :: pdensac ! accumulation mode density [ kg / m**3 ]
real(r8), intent(out) :: betaij0, betaij2i, betaij2j, betaij3, &
betaii0, betaii2, betajj0, betajj2
! *** local parameters
real(r8), parameter :: p0 = 101325.0 ! standard surface pressure [ pa ]
real(r8), parameter :: t0 = 288.15 ! standard surface temperature [ k ]
real(r8), parameter :: avo = 6.0221367e23 ! avogadro's constant [ 1/mol ]
real(r8), parameter :: rgasuniv = 8.314510 ! universal gas constant [ j/mol-k ]
real(r8), parameter :: boltz = rgasuniv/avo ! boltzmann's constant [ j / k]
real(r8), parameter :: two3 = 2.0/3.0
! *** local variables
real(r8) amu ! atmospheric dynamic viscosity [ kg/m s ]
real(r8) sqrt_temp ! square root of ambient temperature
real(r8) lamda ! mean free path [ m ]
! *** intramodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
real(r8) batat( 2 ) ! aitken mode
real(r8) bacac( 2 ) ! accumulation mode
! *** intermodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
real(r8) batac( 2 ) ! aitken to accumulation
real(r8) bacat( 2 ) ! accumulation from aitken
! *** intermodal coagulation rate [ m**3/s ] ( 3rd moment )
real(r8) c3ij ! aitken to accumulation
! *** 3rd moment intermodal transfer rate by coagulation
real(r8) c30atac ! aitken to accumulation
! *** near continnuum regime (independent of mode)
real(r8) knc ! knc = two3 * boltz * airtemp / amu
! *** free molecular regime (depends upon modal density)
real(r8) kfmat ! kfmat = sqrt(3.0*boltz*airtemp/pdensat)
real(r8) kfmac ! kfmac = sqrt(3.0*boltz*airtemp/pdensac)
real(r8) kfmatac ! kfmatac = sqrt( 6.0 * boltz * airtemp /
! ( pdensat + pdensac ) )
real(r8) dumacc2, dumatk2, dumatk3
sqrt_temp = sqrt( airtemp)
! *** calculate mean free path [ m ]:
! 6.6328e-8 is the sea level value given in table i.2.8
! on page 10 of u.s. standard atmosphere 1962
lamda = 6.6328e-8 * p0 * airtemp / ( t0 * airprs )
! *** calculate dynamic viscosity [ kg m**-1 s**-1 ]:
! u.s. standard atmosphere 1962 page 14 expression
! for dynamic viscosity is:
! dynamic viscosity = beta * t * sqrt(t) / ( t + s)
! where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )
! *** coagulation
! calculate coagulation coefficients using a method dictated by
! the value of fastcoag_flag. if true, the computationally-
! efficient getcoags routine is used. if false, the more intensive
! gauss-hermite numerical quadrature method is used. see section
! 2.1 of bhave et al. (2004) for further discussion.
! *** calculate term used in equation a6 of binkowski & shankar (1995)
knc = two3 * boltz * airtemp / amu
! *** calculate terms used in equation a5 of binkowski & shankar (1995)
kfmat = sqrt( 3.0 * boltz * airtemp / pdensat )
kfmac = sqrt( 3.0 * boltz * airtemp / pdensac )
kfmatac = sqrt( 6.0 * boltz * airtemp / ( pdensat + pdensac ) )
! *** transfer of number to accumulation mode from aitken mode is zero
bacat(1) = 0.0
! *** calculate intermodal and intramodal coagulation coefficients
! for zeroth and second moments, and intermodal coagulation
! coefficient for third moment
call getcoags( lamda, kfmatac, kfmat, kfmac, knc, &
dgatk, dgacc, sgatk, sgacc, &
xxlsgat, xxlsgac, &
batat(2), batat(1), bacac(2), bacac(1), &
batac(2), bacat(2), batac(1), c3ij )
! convert from the "cmaq" coag rate parameters
! to the "mirage2" parameters
dumacc2 = ( (dgacc**2) * exp( 2.0*xxlsgac*xxlsgac ) )
dumatk2 = ( (dgatk**2) * exp( 2.0*xxlsgat*xxlsgat ) )
dumatk3 = ( (dgatk**3) * exp( 4.5*xxlsgat*xxlsgat ) )
betaii0 = max( 0.0_r8, batat(1) )
betajj0 = max( 0.0_r8, bacac(1) )
betaij0 = max( 0.0_r8, batac(1) )
betaij3 = max( 0.0_r8, c3ij / dumatk3 )
betajj2 = max( 0.0_r8, bacac(2) / dumacc2 )
betaii2 = max( 0.0_r8, batat(2) / dumatk2 )
betaij2i = max( 0.0_r8, batac(2) / dumatk2 )
betaij2j = max( 0.0_r8, bacat(2) / dumatk2 )
return
end subroutine getcoags_wrapper_f
! //////////////////////////////////////////////////////////////////
! subroutine getcoags calculates the coagulation rates using a new
! approximate algorithm for the 2nd moment. the 0th and 3rd moments
! are done by analytic expressions from whitby et al. (1991). the
! correction factors are also similar to those from whitby et al.
! (1991), but are derived from the gauss-hermite numerical
! quadratures used by binkowski and roselle (2003).
!
! called from aerostep as:
! call getcoags( lamda, kfmatac, kfmat, kfmac, knc,
! dgat,dgac, sgatk, sgacc, xxlsgat,xxlsgac,
! batat(2), batat(1), bacac(2), bacac(1),
! batac(2), bacat(2), batac(1), c3ij )
! where all input and outputs are real*8
!
! revision history:
! fsb 08/25/03 coded by dr. francis s. binkowksi
!
! fsb 08/25/04 added in-line documentation
!
! rce 04/15/2007
! code taken from cmaq v4.6 code; converted to f90;
! added "intent" to subr arguments;
! renamed "r4" & "r8" variables to "rx4" & "rx8";
! changed "real*N" declarations to "real(rN)" (N = 4 or 8)
!
! references:
! 1. whitby, e. r., p. h. mcmurry, u. shankar, and f. s. binkowski,
! modal aerosol dynamics modeling, rep. 600/3-91/020, atmospheric
! research and exposure assessment laboratory,
! u.s. environmental protection agency, research triangle park, n.c.,
! (ntis pb91-161729/as), 1991
!
! 2. binkowski, f.s. an u. shankar, the regional particulate matter
! model 1. model decsription and preliminary results, journal of
! geophysical research, 100, d12, pp 26,191-26,209,
! december 20, 1995.
!
! 3. binkowski, f.s. and s.j. roselle, models-3 community
! multiscale air quality (cmaq) model aerosol component 1:
! model description. j. geophys. res., vol 108, no d6, 4183
! doi:10.1029/2001jd001409, 2003.
subroutine getcoags( lamda, kfmatac, kfmat, kfmac, knc, &
dgatk, dgacc, sgatk, sgacc, xxlsgat,xxlsgac, &
qs11, qn11, qs22, qn22, &
qs12, qs21, qn12, qv12 )
implicit none
real(r8), intent(in) :: lamda ! mean free path [ m ]
! *** coefficients for free molecular regime
real(r8), intent(in) :: kfmat ! aitken mode
real(r8), intent(in) :: kfmac ! accumulation mode
real(r8), intent(in) :: kfmatac ! aitken to accumulation mode
real(r8), intent(in) :: knc ! coefficient for near continnuum regime
! *** modal geometric mean diameters: [ m ]
real(r8), intent(in) :: dgatk ! aitken mode
real(r8), intent(in) :: dgacc ! accumulation mode
! *** modal geometric standard deviation
real(r8), intent(in) :: sgatk ! atken mode
real(r8), intent(in) :: sgacc ! accumulation mode
! *** natural log of modal geometric standard deviation
real(r8), intent(in) :: xxlsgat ! aitken mode
real(r8), intent(in) :: xxlsgac ! accumulation mode
! *** coagulation coefficients
real(r8), intent(out) :: qs11, qn11, qs22, qn22, &
qs12, qs21, qn12, qv12
integer ibeta, n1, n2a, n2n ! indices for correction factors
real(r8) i1fm_at
real(r8) i1nc_at
real(r8) i1_at
real(r8) i1fm_ac
real(r8) i1nc_ac
real(r8) i1_ac
real(r8) i1fm
real(r8) i1nc
real(r8) i1
real(r8) constii
real(r8) kngat, kngac
real(r8) one, two, half
parameter( one = 1.0d0, two = 2.0d0, half = 0.5d0 )
real(r8) a
! parameter( a = 2.492d0)
parameter( a = 1.246d0)
real two3rds
parameter( two3rds = 2.d0 / 3.d0)
real(r8) sqrttwo ! sqrt(two)
real(r8) dlgsqt2 ! 1/ln( sqrt( 2 ) )
real(r8) esat01 ! aitken mode exp( log^2( sigmag )/8 )
real(r8) esac01 ! accumulation mode exp( log^2( sigmag )/8 )
real(r8) esat04
real(r8) esac04
real(r8) esat05
real(r8) esac05
real(r8) esat08
real(r8) esac08
real(r8) esat09
real(r8) esac09
real(r8) esat16
real(r8) esac16
real(r8) esat20
real(r8) esac20
real(r8) esat24
real(r8) esac24
real(r8) esat25
real(r8) esac25
real(r8) esat36
real(r8) esac36
real(r8) esat49
real(r8) esat64
real(r8) esac64
real(r8) esat100
real(r8) dgat2, dgac2, dgat3, dgac3
real(r8) sqdgat, sqdgac
real(r8) sqdgat5, sqdgac5
real(r8) sqdgat7
real(r8) r, r2, r3, rx4, r5, r6, rx8
real(r8) ri1, ri2, ri3, ri4
real(r8) rat
real(r8) coagfm0, coagnc0
real(r8) coagfm3, coagnc3
real(r8) coagfm_at, coagfm_ac
real(r8) coagnc_at, coagnc_ac
real(r8) coagatat0
real(r8) coagacac0
real(r8) coagatat2
real(r8) coagacac2
real(r8) coagatac0, coagatac3
real(r8) coagatac2
real(r8) coagacat2
real(r8) xm2at, xm3at, xm2ac, xm3ac
! *** correction factors for coagulation rates
real(r4), save :: bm0( 10 ) ! m0 intramodal fm - rpm values
real(r4), save :: bm0ij( 10, 10, 10 ) ! m0 intermodal fm
real(r4), save :: bm3i( 10, 10, 10 ) ! m3 intermodal fm- rpm values
real(r4), save :: bm2ii(10) ! m2 intramodal fm
real(r4), save :: bm2iitt(10) ! m2 intramodal total
real(r4), save :: bm2ij(10,10,10) ! m2 intermodal fm i to j
real(r4), save :: bm2ji(10,10,10) ! m2 total intermodal j from i
! *** populate the arrays for the correction factors.
! rpm 0th moment correction factors for unimodal fm coagulation rates
data bm0 / &
0.707106785165097, 0.726148960080488, 0.766430744110958, &
0.814106389441342, 0.861679526483207, 0.903600509090092, &
0.936578814219156, 0.960098926735545, 0.975646823342881, &
0.985397173215326 /
! fsb new fm correction factors for m0 intermodal coagulation
data (bm0ij ( 1, 1,ibeta), ibeta = 1,10) / &
0.628539, 0.639610, 0.664514, 0.696278, 0.731558, &
0.768211, 0.804480, 0.838830, 0.870024, 0.897248/
data (bm0ij ( 1, 2,ibeta), ibeta = 1,10) / &
0.639178, 0.649966, 0.674432, 0.705794, 0.740642, &
0.776751, 0.812323, 0.845827, 0.876076, 0.902324/
data (bm0ij ( 1, 3,ibeta), ibeta = 1,10) / &
0.663109, 0.673464, 0.697147, 0.727637, 0.761425, &
0.796155, 0.829978, 0.861419, 0.889424, 0.913417/
data (bm0ij ( 1, 4,ibeta), ibeta = 1,10) / &
0.693693, 0.703654, 0.726478, 0.755786, 0.787980, &
0.820626, 0.851898, 0.880459, 0.905465, 0.926552/
data (bm0ij ( 1, 5,ibeta), ibeta = 1,10) / &
0.727803, 0.737349, 0.759140, 0.786870, 0.816901, &
0.846813, 0.874906, 0.900060, 0.921679, 0.939614/
data (bm0ij ( 1, 6,ibeta), ibeta = 1,10) / &
0.763461, 0.772483, 0.792930, 0.818599, 0.845905, &
0.872550, 0.897051, 0.918552, 0.936701, 0.951528/
data (bm0ij ( 1, 7,ibeta), ibeta = 1,10) / &
0.799021, 0.807365, 0.826094, 0.849230, 0.873358, &
0.896406, 0.917161, 0.935031, 0.949868, 0.961828/
data (bm0ij ( 1, 8,ibeta), ibeta = 1,10) / &
0.833004, 0.840514, 0.857192, 0.877446, 0.898147, &
0.917518, 0.934627, 0.949106, 0.960958, 0.970403/
data (bm0ij ( 1, 9,ibeta), ibeta = 1,10) / &
0.864172, 0.870734, 0.885153, 0.902373, 0.919640, &
0.935494, 0.949257, 0.960733, 0.970016, 0.977346/
data (bm0ij ( 1, 10,ibeta), ibeta = 1,10) / &
0.891658, 0.897227, 0.909343, 0.923588, 0.937629, &
0.950307, 0.961151, 0.970082, 0.977236, 0.982844/
data (bm0ij ( 2, 1,ibeta), ibeta = 1,10) / &
0.658724, 0.670587, 0.697539, 0.731890, 0.769467, &
0.807391, 0.843410, 0.875847, 0.903700, 0.926645/
data (bm0ij ( 2, 2,ibeta), ibeta = 1,10) / &
0.667070, 0.678820, 0.705538, 0.739591, 0.776758, &
0.814118, 0.849415, 0.881020, 0.908006, 0.930121/
data (bm0ij ( 2, 3,ibeta), ibeta = 1,10) / &
0.686356, 0.697839, 0.723997, 0.757285, 0.793389, &
0.829313, 0.862835, 0.892459, 0.917432, 0.937663/
data (bm0ij ( 2, 4,ibeta), ibeta = 1,10) / &
0.711425, 0.722572, 0.747941, 0.780055, 0.814518, &
0.848315, 0.879335, 0.906290, 0.928658, 0.946526/
data (bm0ij ( 2, 5,ibeta), ibeta = 1,10) / &
0.739575, 0.750307, 0.774633, 0.805138, 0.837408, &
0.868504, 0.896517, 0.920421, 0.939932, 0.955299/
data (bm0ij ( 2, 6,ibeta), ibeta = 1,10) / &
0.769143, 0.779346, 0.802314, 0.830752, 0.860333, &
0.888300, 0.913014, 0.933727, 0.950370, 0.963306/
data (bm0ij ( 2, 7,ibeta), ibeta = 1,10) / &
0.798900, 0.808431, 0.829700, 0.855653, 0.882163, &
0.906749, 0.928075, 0.945654, 0.959579, 0.970280/
data (bm0ij ( 2, 8,ibeta), ibeta = 1,10) / &
0.827826, 0.836542, 0.855808, 0.878954, 0.902174, &
0.923316, 0.941345, 0.955989, 0.967450, 0.976174/
data (bm0ij ( 2, 9,ibeta), ibeta = 1,10) / &
0.855068, 0.862856, 0.879900, 0.900068, 0.919956, &
0.937764, 0.952725, 0.964726, 0.974027, 0.981053/
data (bm0ij ( 2, 10,ibeta), ibeta = 1,10) / &
0.879961, 0.886755, 0.901484, 0.918665, 0.935346, &
0.950065, 0.962277, 0.971974, 0.979432, 0.985033/
data (bm0ij ( 3, 1,ibeta), ibeta = 1,10) / &
0.724166, 0.735474, 0.761359, 0.794045, 0.828702, &
0.862061, 0.891995, 0.917385, 0.937959, 0.954036/
data (bm0ij ( 3, 2,ibeta), ibeta = 1,10) / &
0.730416, 0.741780, 0.767647, 0.800116, 0.834344, &
0.867093, 0.896302, 0.920934, 0.940790, 0.956237/
data (bm0ij ( 3, 3,ibeta), ibeta = 1,10) / &
0.745327, 0.756664, 0.782255, 0.814026, 0.847107, &
0.878339, 0.905820, 0.928699, 0.946931, 0.960977/
data (bm0ij ( 3, 4,ibeta), ibeta = 1,10) / &
0.765195, 0.776312, 0.801216, 0.831758, 0.863079, &
0.892159, 0.917319, 0.937939, 0.954145, 0.966486/
data (bm0ij ( 3, 5,ibeta), ibeta = 1,10) / &
0.787632, 0.798347, 0.822165, 0.850985, 0.880049, &
0.906544, 0.929062, 0.947218, 0.961288, 0.971878/
data (bm0ij ( 3, 6,ibeta), ibeta = 1,10) / &
0.811024, 0.821179, 0.843557, 0.870247, 0.896694, &
0.920365, 0.940131, 0.955821, 0.967820, 0.976753/
data (bm0ij ( 3, 7,ibeta), ibeta = 1,10) / &
0.834254, 0.843709, 0.864356, 0.888619, 0.912245, &
0.933019, 0.950084, 0.963438, 0.973530, 0.980973/
data (bm0ij ( 3, 8,ibeta), ibeta = 1,10) / &
0.856531, 0.865176, 0.883881, 0.905544, 0.926290, &
0.944236, 0.958762, 0.969988, 0.978386, 0.984530/
data (bm0ij ( 3, 9,ibeta), ibeta = 1,10) / &
0.877307, 0.885070, 0.901716, 0.920729, 0.938663, &
0.953951, 0.966169, 0.975512, 0.982442, 0.987477/
data (bm0ij ( 3, 10,ibeta), ibeta = 1,10) / &
0.896234, 0.903082, 0.917645, 0.934069, 0.949354, &
0.962222, 0.972396, 0.980107, 0.985788, 0.989894/
data (bm0ij ( 4, 1,ibeta), ibeta = 1,10) / &
0.799294, 0.809144, 0.831293, 0.858395, 0.885897, &
0.911031, 0.932406, 0.949642, 0.963001, 0.973062/
data (bm0ij ( 4, 2,ibeta), ibeta = 1,10) / &
0.804239, 0.814102, 0.836169, 0.862984, 0.890003, &
0.914535, 0.935274, 0.951910, 0.964748, 0.974381/
data (bm0ij ( 4, 3,ibeta), ibeta = 1,10) / &
0.815910, 0.825708, 0.847403, 0.873389, 0.899185, &
0.922275, 0.941543, 0.956826, 0.968507, 0.977204/
data (bm0ij ( 4, 4,ibeta), ibeta = 1,10) / &
0.831348, 0.840892, 0.861793, 0.886428, 0.910463, &
0.931614, 0.948993, 0.962593, 0.972872, 0.980456/
data (bm0ij ( 4, 5,ibeta), ibeta = 1,10) / &
0.848597, 0.857693, 0.877402, 0.900265, 0.922180, &
0.941134, 0.956464, 0.968298, 0.977143, 0.983611/
data (bm0ij ( 4, 6,ibeta), ibeta = 1,10) / &
0.866271, 0.874764, 0.892984, 0.913796, 0.933407, &
0.950088, 0.963380, 0.973512, 0.981006, 0.986440/
data (bm0ij ( 4, 7,ibeta), ibeta = 1,10) / &
0.883430, 0.891216, 0.907762, 0.926388, 0.943660, &
0.958127, 0.969499, 0.978070, 0.984351, 0.988872/
data (bm0ij ( 4, 8,ibeta), ibeta = 1,10) / &
0.899483, 0.906505, 0.921294, 0.937719, 0.952729, &
0.965131, 0.974762, 0.981950, 0.987175, 0.990912/
data (bm0ij ( 4, 9,ibeta), ibeta = 1,10) / &
0.914096, 0.920337, 0.933373, 0.947677, 0.960579, &
0.971111, 0.979206, 0.985196, 0.989520, 0.992597/
data (bm0ij ( 4, 10,ibeta), ibeta = 1,10) / &
0.927122, 0.932597, 0.943952, 0.956277, 0.967268, &
0.976147, 0.982912, 0.987882, 0.991450, 0.993976/
data (bm0ij ( 5, 1,ibeta), ibeta = 1,10) / &
0.865049, 0.872851, 0.889900, 0.909907, 0.929290, &
0.946205, 0.959991, 0.970706, 0.978764, 0.984692/
data (bm0ij ( 5, 2,ibeta), ibeta = 1,10) / &
0.868989, 0.876713, 0.893538, 0.913173, 0.932080, &
0.948484, 0.961785, 0.972080, 0.979796, 0.985457/
data (bm0ij ( 5, 3,ibeta), ibeta = 1,10) / &
0.878010, 0.885524, 0.901756, 0.920464, 0.938235, &
0.953461, 0.965672, 0.975037, 0.982005, 0.987085/
data (bm0ij ( 5, 4,ibeta), ibeta = 1,10) / &
0.889534, 0.896698, 0.912012, 0.929395, 0.945647, &
0.959366, 0.970227, 0.978469, 0.984547, 0.988950/
data (bm0ij ( 5, 5,ibeta), ibeta = 1,10) / &
0.902033, 0.908713, 0.922848, 0.938648, 0.953186, &
0.965278, 0.974729, 0.981824, 0.987013, 0.990746/
data (bm0ij ( 5, 6,ibeta), ibeta = 1,10) / &
0.914496, 0.920599, 0.933389, 0.947485, 0.960262, &
0.970743, 0.978839, 0.984858, 0.989225, 0.992348/
data (bm0ij ( 5, 7,ibeta), ibeta = 1,10) / &
0.926281, 0.931761, 0.943142, 0.955526, 0.966600, &
0.975573, 0.982431, 0.987485, 0.991128, 0.993718/
data (bm0ij ( 5, 8,ibeta), ibeta = 1,10) / &
0.937029, 0.941877, 0.951868, 0.962615, 0.972112, &
0.979723, 0.985488, 0.989705, 0.992725, 0.994863/
data (bm0ij ( 5, 9,ibeta), ibeta = 1,10) / &
0.946580, 0.950819, 0.959494, 0.968732, 0.976811, &
0.983226, 0.988047, 0.991550, 0.994047, 0.995806/
data (bm0ij ( 5, 10,ibeta), ibeta = 1,10) / &
0.954909, 0.958581, 0.966049, 0.973933, 0.980766, &
0.986149, 0.990166, 0.993070, 0.995130, 0.996577/
data (bm0ij ( 6, 1,ibeta), ibeta = 1,10) / &
0.914182, 0.919824, 0.931832, 0.945387, 0.957999, &
0.968606, 0.976982, 0.983331, 0.988013, 0.991407/
data (bm0ij ( 6, 2,ibeta), ibeta = 1,10) / &
0.917139, 0.922665, 0.934395, 0.947580, 0.959792, &
0.970017, 0.978062, 0.984138, 0.988609, 0.991843/
data (bm0ij ( 6, 3,ibeta), ibeta = 1,10) / &
0.923742, 0.928990, 0.940064, 0.952396, 0.963699, &
0.973070, 0.980381, 0.985866, 0.989878, 0.992768/
data (bm0ij ( 6, 4,ibeta), ibeta = 1,10) / &
0.931870, 0.936743, 0.946941, 0.958162, 0.968318, &
0.976640, 0.983069, 0.987853, 0.991330, 0.993822/
data (bm0ij ( 6, 5,ibeta), ibeta = 1,10) / &
0.940376, 0.944807, 0.954004, 0.963999, 0.972928, &
0.980162, 0.985695, 0.989779, 0.992729, 0.994833/
data (bm0ij ( 6, 6,ibeta), ibeta = 1,10) / &
0.948597, 0.952555, 0.960703, 0.969454, 0.977181, &
0.983373, 0.988067, 0.991507, 0.993977, 0.995730/
data (bm0ij ( 6, 7,ibeta), ibeta = 1,10) / &
0.956167, 0.959648, 0.966763, 0.974326, 0.980933, &
0.986177, 0.990121, 0.992993, 0.995045, 0.996495/
data (bm0ij ( 6, 8,ibeta), ibeta = 1,10) / &
0.962913, 0.965937, 0.972080, 0.978552, 0.984153, &
0.988563, 0.991857, 0.994242, 0.995938, 0.997133/
data (bm0ij ( 6, 9,ibeta), ibeta = 1,10) / &
0.968787, 0.971391, 0.976651, 0.982148, 0.986869, &
0.990560, 0.993301, 0.995275, 0.996675, 0.997657/
data (bm0ij ( 6, 10,ibeta), ibeta = 1,10) / &
0.973822, 0.976047, 0.980523, 0.985170, 0.989134, &
0.992215, 0.994491, 0.996124, 0.997277, 0.998085/
data (bm0ij ( 7, 1,ibeta), ibeta = 1,10) / &
0.947410, 0.951207, 0.959119, 0.967781, 0.975592, &
0.981981, 0.986915, 0.990590, 0.993266, 0.995187/
data (bm0ij ( 7, 2,ibeta), ibeta = 1,10) / &
0.949477, 0.953161, 0.960824, 0.969187, 0.976702, &
0.982831, 0.987550, 0.991057, 0.993606, 0.995434/
data (bm0ij ( 7, 3,ibeta), ibeta = 1,10) / &
0.954008, 0.957438, 0.964537, 0.972232, 0.979095, &
0.984653, 0.988907, 0.992053, 0.994330, 0.995958/
data (bm0ij ( 7, 4,ibeta), ibeta = 1,10) / &
0.959431, 0.962539, 0.968935, 0.975808, 0.981882, &
0.986759, 0.990466, 0.993190, 0.995153, 0.996552/
data (bm0ij ( 7, 5,ibeta), ibeta = 1,10) / &
0.964932, 0.967693, 0.973342, 0.979355, 0.984620, &
0.988812, 0.991974, 0.994285, 0.995943, 0.997119/
data (bm0ij ( 7, 6,ibeta), ibeta = 1,10) / &
0.970101, 0.972517, 0.977428, 0.982612, 0.987110, &
0.990663, 0.993326, 0.995261, 0.996644, 0.997621/
data (bm0ij ( 7, 7,ibeta), ibeta = 1,10) / &
0.974746, 0.976834, 0.981055, 0.985475, 0.989280, &
0.992265, 0.994488, 0.996097, 0.997241, 0.998048/
data (bm0ij ( 7, 8,ibeta), ibeta = 1,10) / &
0.978804, 0.980591, 0.984187, 0.987927, 0.991124, &
0.993617, 0.995464, 0.996795, 0.997739, 0.998403/
data (bm0ij ( 7, 9,ibeta), ibeta = 1,10) / &
0.982280, 0.983799, 0.986844, 0.989991, 0.992667, &
0.994742, 0.996273, 0.997372, 0.998149, 0.998695/
data (bm0ij ( 7, 10,ibeta), ibeta = 1,10) / &
0.985218, 0.986503, 0.989071, 0.991711, 0.993945, &
0.995669, 0.996937, 0.997844, 0.998484, 0.998932/
data (bm0ij ( 8, 1,ibeta), ibeta = 1,10) / &
0.968507, 0.970935, 0.975916, 0.981248, 0.985947, &
0.989716, 0.992580, 0.994689, 0.996210, 0.997297/
data (bm0ij ( 8, 2,ibeta), ibeta = 1,10) / &
0.969870, 0.972210, 0.977002, 0.982119, 0.986619, &
0.990219, 0.992951, 0.994958, 0.996405, 0.997437/
data (bm0ij ( 8, 3,ibeta), ibeta = 1,10) / &
0.972820, 0.974963, 0.979339, 0.983988, 0.988054, &
0.991292, 0.993738, 0.995529, 0.996817, 0.997734/
data (bm0ij ( 8, 4,ibeta), ibeta = 1,10) / &
0.976280, 0.978186, 0.982060, 0.986151, 0.989706, &
0.992520, 0.994636, 0.996179, 0.997284, 0.998069/
data (bm0ij ( 8, 5,ibeta), ibeta = 1,10) / &
0.979711, 0.981372, 0.984735, 0.988263, 0.991309, &
0.993706, 0.995499, 0.996801, 0.997730, 0.998389/
data (bm0ij ( 8, 6,ibeta), ibeta = 1,10) / &
0.982863, 0.984292, 0.987172, 0.990174, 0.992750, &
0.994766, 0.996266, 0.997352, 0.998125, 0.998670/
data (bm0ij ( 8, 7,ibeta), ibeta = 1,10) / &
0.985642, 0.986858, 0.989301, 0.991834, 0.993994, &
0.995676, 0.996923, 0.997822, 0.998460, 0.998910/
data (bm0ij ( 8, 8,ibeta), ibeta = 1,10) / &
0.988029, 0.989058, 0.991116, 0.993240, 0.995043, &
0.996440, 0.997472, 0.998214, 0.998739, 0.999108/
data (bm0ij ( 8, 9,ibeta), ibeta = 1,10) / &
0.990046, 0.990912, 0.992640, 0.994415, 0.995914, &
0.997073, 0.997925, 0.998536, 0.998968, 0.999271/
data (bm0ij ( 8, 10,ibeta), ibeta = 1,10) / &
0.991732, 0.992459, 0.993906, 0.995386, 0.996633, &
0.997592, 0.998296, 0.998799, 0.999154, 0.999403/
data (bm0ij ( 9, 1,ibeta), ibeta = 1,10) / &
0.981392, 0.982893, 0.985938, 0.989146, 0.991928, &
0.994129, 0.995783, 0.996991, 0.997857, 0.998473/
data (bm0ij ( 9, 2,ibeta), ibeta = 1,10) / &
0.982254, 0.983693, 0.986608, 0.989673, 0.992328, &
0.994424, 0.995998, 0.997146, 0.997969, 0.998553/
data (bm0ij ( 9, 3,ibeta), ibeta = 1,10) / &
0.984104, 0.985407, 0.988040, 0.990798, 0.993178, &
0.995052, 0.996454, 0.997474, 0.998204, 0.998722/
data (bm0ij ( 9, 4,ibeta), ibeta = 1,10) / &
0.986243, 0.987386, 0.989687, 0.992087, 0.994149, &
0.995765, 0.996971, 0.997846, 0.998470, 0.998913/
data (bm0ij ( 9, 5,ibeta), ibeta = 1,10) / &
0.988332, 0.989313, 0.991284, 0.993332, 0.995082, &
0.996449, 0.997465, 0.998200, 0.998723, 0.999093/
data (bm0ij ( 9, 6,ibeta), ibeta = 1,10) / &
0.990220, 0.991053, 0.992721, 0.994445, 0.995914, &
0.997056, 0.997902, 0.998513, 0.998947, 0.999253/
data (bm0ij ( 9, 7,ibeta), ibeta = 1,10) / &
0.991859, 0.992561, 0.993961, 0.995403, 0.996626, &
0.997574, 0.998274, 0.998778, 0.999136, 0.999387/
data (bm0ij ( 9, 8,ibeta), ibeta = 1,10) / &
0.993250, 0.993837, 0.995007, 0.996208, 0.997223, &
0.998007, 0.998584, 0.998999, 0.999293, 0.999499/
data (bm0ij ( 9, 9,ibeta), ibeta = 1,10) / &
0.994413, 0.994903, 0.995878, 0.996876, 0.997716, &
0.998363, 0.998839, 0.999180, 0.999421, 0.999591/
data (bm0ij ( 9, 10,ibeta), ibeta = 1,10) / &
0.995376, 0.995785, 0.996597, 0.997425, 0.998121, &
0.998655, 0.999048, 0.999328, 0.999526, 0.999665/
data (bm0ij ( 10, 1,ibeta), ibeta = 1,10) / &
0.989082, 0.989991, 0.991819, 0.993723, 0.995357, &
0.996637, 0.997592, 0.998286, 0.998781, 0.999132/
data (bm0ij ( 10, 2,ibeta), ibeta = 1,10) / &
0.989613, 0.990480, 0.992224, 0.994039, 0.995594, &
0.996810, 0.997717, 0.998375, 0.998845, 0.999178/
data (bm0ij ( 10, 3,ibeta), ibeta = 1,10) / &
0.990744, 0.991523, 0.993086, 0.994708, 0.996094, &
0.997176, 0.997981, 0.998564, 0.998980, 0.999274/
data (bm0ij ( 10, 4,ibeta), ibeta = 1,10) / &
0.992041, 0.992716, 0.994070, 0.995470, 0.996662, &
0.997591, 0.998280, 0.998778, 0.999133, 0.999383/
data (bm0ij ( 10, 5,ibeta), ibeta = 1,10) / &
0.993292, 0.993867, 0.995015, 0.996199, 0.997205, &
0.997985, 0.998564, 0.998981, 0.999277, 0.999487/
data (bm0ij ( 10, 6,ibeta), ibeta = 1,10) / &
0.994411, 0.994894, 0.995857, 0.996847, 0.997685, &
0.998334, 0.998814, 0.999159, 0.999404, 0.999577/
data (bm0ij ( 10, 7,ibeta), ibeta = 1,10) / &
0.995373, 0.995776, 0.996577, 0.997400, 0.998094, &
0.998630, 0.999026, 0.999310, 0.999512, 0.999654/
data (bm0ij ( 10, 8,ibeta), ibeta = 1,10) / &
0.996181, 0.996516, 0.997181, 0.997861, 0.998435, &
0.998877, 0.999202, 0.999435, 0.999601, 0.999717/
data (bm0ij ( 10, 9,ibeta), ibeta = 1,10) / &
0.996851, 0.997128, 0.997680, 0.998242, 0.998715, &
0.999079, 0.999346, 0.999538, 0.999673, 0.999769/
data (bm0ij ( 10, 10,ibeta), ibeta = 1,10) / &
0.997402, 0.997632, 0.998089, 0.998554, 0.998945, &
0.999244, 0.999464, 0.999622, 0.999733, 0.999811/
! rpm.... 3rd moment nuclei mode corr. fac. for bimodal fm coag rate
data (bm3i( 1, 1,ibeta ), ibeta=1,10)/ &
0.70708,0.71681,0.73821,0.76477,0.79350,0.82265,0.85090,0.87717, &
0.90069,0.92097/
data (bm3i( 1, 2,ibeta ), ibeta=1,10)/ &
0.72172,0.73022,0.74927,0.77324,0.79936,0.82601,0.85199,0.87637, &
0.89843,0.91774/
data (bm3i( 1, 3,ibeta ), ibeta=1,10)/ &
0.78291,0.78896,0.80286,0.82070,0.84022,0.85997,0.87901,0.89669, &
0.91258,0.92647/
data (bm3i( 1, 4,ibeta ), ibeta=1,10)/ &
0.87760,0.88147,0.89025,0.90127,0.91291,0.92420,0.93452,0.94355, &
0.95113,0.95726/
data (bm3i( 1, 5,ibeta ), ibeta=1,10)/ &
0.94988,0.95184,0.95612,0.96122,0.96628,0.97085,0.97467,0.97763, &
0.97971,0.98089/
data (bm3i( 1, 6,ibeta ), ibeta=1,10)/ &
0.98318,0.98393,0.98551,0.98728,0.98889,0.99014,0.99095,0.99124, &
0.99100,0.99020/
data (bm3i( 1, 7,ibeta ), ibeta=1,10)/ &
0.99480,0.99504,0.99551,0.99598,0.99629,0.99635,0.99611,0.99550, &
0.99450,0.99306/
data (bm3i( 1, 8,ibeta ), ibeta=1,10)/ &
0.99842,0.99848,0.99858,0.99861,0.99850,0.99819,0.99762,0.99674, &
0.99550,0.99388/
data (bm3i( 1, 9,ibeta ), ibeta=1,10)/ &
0.99951,0.99951,0.99949,0.99939,0.99915,0.99872,0.99805,0.99709, &
0.99579,0.99411/
data (bm3i( 1,10,ibeta ), ibeta=1,10)/ &
0.99984,0.99982,0.99976,0.99962,0.99934,0.99888,0.99818,0.99719, &
0.99587,0.99417/
data (bm3i( 2, 1,ibeta ), ibeta=1,10)/ &
0.72957,0.73993,0.76303,0.79178,0.82245,0.85270,0.88085,0.90578, &
0.92691,0.94415/
data (bm3i( 2, 2,ibeta ), ibeta=1,10)/ &
0.72319,0.73320,0.75547,0.78323,0.81307,0.84287,0.87107,0.89651, &
0.91852,0.93683/
data (bm3i( 2, 3,ibeta ), ibeta=1,10)/ &
0.74413,0.75205,0.76998,0.79269,0.81746,0.84258,0.86685,0.88938, &
0.90953,0.92695/
data (bm3i( 2, 4,ibeta ), ibeta=1,10)/ &
0.82588,0.83113,0.84309,0.85825,0.87456,0.89072,0.90594,0.91972, &
0.93178,0.94203/
data (bm3i( 2, 5,ibeta ), ibeta=1,10)/ &
0.91886,0.92179,0.92831,0.93624,0.94434,0.95192,0.95856,0.96409, &
0.96845,0.97164/
data (bm3i( 2, 6,ibeta ), ibeta=1,10)/ &
0.97129,0.97252,0.97515,0.97818,0.98108,0.98354,0.98542,0.98665, &
0.98721,0.98709/
data (bm3i( 2, 7,ibeta ), ibeta=1,10)/ &
0.99104,0.99145,0.99230,0.99320,0.99394,0.99439,0.99448,0.99416, &
0.99340,0.99217/
data (bm3i( 2, 8,ibeta ), ibeta=1,10)/ &
0.99730,0.99741,0.99763,0.99779,0.99782,0.99762,0.99715,0.99636, &
0.99519,0.99363/
data (bm3i( 2, 9,ibeta ), ibeta=1,10)/ &
0.99917,0.99919,0.99921,0.99915,0.99895,0.99856,0.99792,0.99698, &
0.99570,0.99404/
data (bm3i( 2,10,ibeta ), ibeta=1,10)/ &
0.99973,0.99973,0.99968,0.99955,0.99928,0.99883,0.99814,0.99716, &
0.99584,0.99415/
data (bm3i( 3, 1,ibeta ), ibeta=1,10)/ &
0.78358,0.79304,0.81445,0.84105,0.86873,0.89491,0.91805,0.93743, &
0.95300,0.96510/
data (bm3i( 3, 2,ibeta ), ibeta=1,10)/ &
0.76412,0.77404,0.79635,0.82404,0.85312,0.88101,0.90610,0.92751, &
0.94500,0.95879/
data (bm3i( 3, 3,ibeta ), ibeta=1,10)/ &
0.74239,0.75182,0.77301,0.79956,0.82809,0.85639,0.88291,0.90658, &
0.92683,0.94350/
data (bm3i( 3, 4,ibeta ), ibeta=1,10)/ &
0.78072,0.78758,0.80317,0.82293,0.84437,0.86589,0.88643,0.90526, &
0.92194,0.93625/
data (bm3i( 3, 5,ibeta ), ibeta=1,10)/ &
0.87627,0.88044,0.88981,0.90142,0.91357,0.92524,0.93585,0.94510, &
0.95285,0.95911/
data (bm3i( 3, 6,ibeta ), ibeta=1,10)/ &
0.95176,0.95371,0.95796,0.96297,0.96792,0.97233,0.97599,0.97880, &
0.98072,0.98178/
data (bm3i( 3, 7,ibeta ), ibeta=1,10)/ &
0.98453,0.98523,0.98670,0.98833,0.98980,0.99092,0.99160,0.99179, &
0.99145,0.99058/
data (bm3i( 3, 8,ibeta ), ibeta=1,10)/ &
0.99534,0.99555,0.99597,0.99637,0.99662,0.99663,0.99633,0.99569, &
0.99465,0.99318/
data (bm3i( 3, 9,ibeta ), ibeta=1,10)/ &
0.99859,0.99864,0.99872,0.99873,0.99860,0.99827,0.99768,0.99679, &
0.99555,0.99391/
data (bm3i( 3,10,ibeta ), ibeta=1,10)/ &
0.99956,0.99956,0.99953,0.99942,0.99918,0.99875,0.99807,0.99711, &
0.99580,0.99412/
data (bm3i( 4, 1,ibeta ), ibeta=1,10)/ &
0.84432,0.85223,0.86990,0.89131,0.91280,0.93223,0.94861,0.96172, &
0.97185,0.97945/
data (bm3i( 4, 2,ibeta ), ibeta=1,10)/ &
0.82299,0.83164,0.85101,0.87463,0.89857,0.92050,0.93923,0.95443, &
0.96629,0.97529/
data (bm3i( 4, 3,ibeta ), ibeta=1,10)/ &
0.77870,0.78840,0.81011,0.83690,0.86477,0.89124,0.91476,0.93460, &
0.95063,0.96316/
data (bm3i( 4, 4,ibeta ), ibeta=1,10)/ &
0.76386,0.77233,0.79147,0.81557,0.84149,0.86719,0.89126,0.91275, &
0.93116,0.94637/
data (bm3i( 4, 5,ibeta ), ibeta=1,10)/ &
0.82927,0.83488,0.84756,0.86346,0.88040,0.89704,0.91257,0.92649, &
0.93857,0.94874/
data (bm3i( 4, 6,ibeta ), ibeta=1,10)/ &
0.92184,0.92481,0.93136,0.93925,0.94724,0.95462,0.96104,0.96634, &
0.97048,0.97348/
data (bm3i( 4, 7,ibeta ), ibeta=1,10)/ &
0.97341,0.97457,0.97706,0.97991,0.98260,0.98485,0.98654,0.98760, &
0.98801,0.98777/
data (bm3i( 4, 8,ibeta ), ibeta=1,10)/ &
0.99192,0.99229,0.99305,0.99385,0.99449,0.99486,0.99487,0.99449, &
0.99367,0.99239/
data (bm3i( 4, 9,ibeta ), ibeta=1,10)/ &
0.99758,0.99768,0.99787,0.99800,0.99799,0.99777,0.99727,0.99645, &
0.99527,0.99369/
data (bm3i( 4,10,ibeta ), ibeta=1,10)/ &
0.99926,0.99928,0.99928,0.99921,0.99900,0.99860,0.99795,0.99701, &
0.99572,0.99405/
data (bm3i( 5, 1,ibeta ), ibeta=1,10)/ &
0.89577,0.90190,0.91522,0.93076,0.94575,0.95876,0.96932,0.97751, &
0.98367,0.98820/
data (bm3i( 5, 2,ibeta ), ibeta=1,10)/ &
0.87860,0.88547,0.90052,0.91828,0.93557,0.95075,0.96319,0.97292, &
0.98028,0.98572/
data (bm3i( 5, 3,ibeta ), ibeta=1,10)/ &
0.83381,0.84240,0.86141,0.88425,0.90707,0.92770,0.94510,0.95906, &
0.96986,0.97798/
data (bm3i( 5, 4,ibeta ), ibeta=1,10)/ &
0.78530,0.79463,0.81550,0.84127,0.86813,0.89367,0.91642,0.93566, &
0.95125,0.96347/
data (bm3i( 5, 5,ibeta ), ibeta=1,10)/ &
0.79614,0.80332,0.81957,0.84001,0.86190,0.88351,0.90368,0.92169, &
0.93718,0.95006/
data (bm3i( 5, 6,ibeta ), ibeta=1,10)/ &
0.88192,0.88617,0.89565,0.90728,0.91931,0.93076,0.94107,0.94997, &
0.95739,0.96333/
data (bm3i( 5, 7,ibeta ), ibeta=1,10)/ &
0.95509,0.95698,0.96105,0.96583,0.97048,0.97460,0.97796,0.98050, &
0.98218,0.98304/
data (bm3i( 5, 8,ibeta ), ibeta=1,10)/ &
0.98596,0.98660,0.98794,0.98943,0.99074,0.99172,0.99227,0.99235, &
0.99192,0.99096/
data (bm3i( 5, 9,ibeta ), ibeta=1,10)/ &
0.99581,0.99600,0.99637,0.99672,0.99691,0.99687,0.99653,0.99585, &
0.99478,0.99329/
data (bm3i( 5,10,ibeta ), ibeta=1,10)/ &
0.99873,0.99878,0.99884,0.99883,0.99869,0.99834,0.99774,0.99684, &
0.99558,0.99394/
data (bm3i( 6, 1,ibeta ), ibeta=1,10)/ &
0.93335,0.93777,0.94711,0.95764,0.96741,0.97562,0.98210,0.98701, &
0.99064,0.99327/
data (bm3i( 6, 2,ibeta ), ibeta=1,10)/ &
0.92142,0.92646,0.93723,0.94947,0.96096,0.97069,0.97842,0.98431, &
0.98868,0.99186/
data (bm3i( 6, 3,ibeta ), ibeta=1,10)/ &
0.88678,0.89351,0.90810,0.92508,0.94138,0.95549,0.96693,0.97578, &
0.98243,0.98731/
data (bm3i( 6, 4,ibeta ), ibeta=1,10)/ &
0.83249,0.84124,0.86051,0.88357,0.90655,0.92728,0.94477,0.95880, &
0.96964,0.97779/
data (bm3i( 6, 5,ibeta ), ibeta=1,10)/ &
0.79593,0.80444,0.82355,0.84725,0.87211,0.89593,0.91735,0.93566, &
0.95066,0.96255/
data (bm3i( 6, 6,ibeta ), ibeta=1,10)/ &
0.84124,0.84695,0.85980,0.87575,0.89256,0.90885,0.92383,0.93704, &
0.94830,0.95761/
data (bm3i( 6, 7,ibeta ), ibeta=1,10)/ &
0.92721,0.93011,0.93647,0.94406,0.95166,0.95862,0.96460,0.96949, &
0.97326,0.97595/
data (bm3i( 6, 8,ibeta ), ibeta=1,10)/ &
0.97573,0.97681,0.97913,0.98175,0.98421,0.98624,0.98772,0.98860, &
0.98885,0.98847/
data (bm3i( 6, 9,ibeta ), ibeta=1,10)/ &
0.99271,0.99304,0.99373,0.99444,0.99499,0.99528,0.99522,0.99477, &
0.99390,0.99258/
data (bm3i( 6,10,ibeta ), ibeta=1,10)/ &
0.99782,0.99791,0.99807,0.99817,0.99813,0.99788,0.99737,0.99653, &
0.99533,0.99374/
data (bm3i( 7, 1,ibeta ), ibeta=1,10)/ &
0.95858,0.96158,0.96780,0.97460,0.98073,0.98575,0.98963,0.99252, &
0.99463,0.99615/
data (bm3i( 7, 2,ibeta ), ibeta=1,10)/ &
0.95091,0.95438,0.96163,0.96962,0.97688,0.98286,0.98751,0.99099, &
0.99353,0.99536/
data (bm3i( 7, 3,ibeta ), ibeta=1,10)/ &
0.92751,0.93233,0.94255,0.95406,0.96473,0.97366,0.98070,0.98602, &
0.98994,0.99278/
data (bm3i( 7, 4,ibeta ), ibeta=1,10)/ &
0.88371,0.89075,0.90595,0.92351,0.94028,0.95474,0.96642,0.97544, &
0.98220,0.98715/
data (bm3i( 7, 5,ibeta ), ibeta=1,10)/ &
0.82880,0.83750,0.85671,0.87980,0.90297,0.92404,0.94195,0.95644, &
0.96772,0.97625/
data (bm3i( 7, 6,ibeta ), ibeta=1,10)/ &
0.81933,0.82655,0.84279,0.86295,0.88412,0.90449,0.92295,0.93890, &
0.95215,0.96281/
data (bm3i( 7, 7,ibeta ), ibeta=1,10)/ &
0.89099,0.89519,0.90448,0.91577,0.92732,0.93820,0.94789,0.95616, &
0.96297,0.96838/
data (bm3i( 7, 8,ibeta ), ibeta=1,10)/ &
0.95886,0.96064,0.96448,0.96894,0.97324,0.97701,0.98004,0.98228, &
0.98371,0.98435/
data (bm3i( 7, 9,ibeta ), ibeta=1,10)/ &
0.98727,0.98786,0.98908,0.99043,0.99160,0.99245,0.99288,0.99285, &
0.99234,0.99131/
data (bm3i( 7,10,ibeta ), ibeta=1,10)/ &
0.99621,0.99638,0.99671,0.99700,0.99715,0.99707,0.99670,0.99599, &
0.99489,0.99338/
data (bm3i( 8, 1,ibeta ), ibeta=1,10)/ &
0.97470,0.97666,0.98064,0.98491,0.98867,0.99169,0.99399,0.99569, &
0.99691,0.99779/
data (bm3i( 8, 2,ibeta ), ibeta=1,10)/ &
0.96996,0.97225,0.97693,0.98196,0.98643,0.99003,0.99279,0.99482, &
0.99630,0.99735/
data (bm3i( 8, 3,ibeta ), ibeta=1,10)/ &
0.95523,0.95848,0.96522,0.97260,0.97925,0.98468,0.98888,0.99200, &
0.99427,0.99590/
data (bm3i( 8, 4,ibeta ), ibeta=1,10)/ &
0.92524,0.93030,0.94098,0.95294,0.96397,0.97317,0.98038,0.98582, &
0.98981,0.99270/
data (bm3i( 8, 5,ibeta ), ibeta=1,10)/ &
0.87576,0.88323,0.89935,0.91799,0.93583,0.95126,0.96377,0.97345, &
0.98072,0.98606/
data (bm3i( 8, 6,ibeta ), ibeta=1,10)/ &
0.83078,0.83894,0.85705,0.87899,0.90126,0.92179,0.93950,0.95404, &
0.96551,0.97430/
data (bm3i( 8, 7,ibeta ), ibeta=1,10)/ &
0.85727,0.86294,0.87558,0.89111,0.90723,0.92260,0.93645,0.94841, &
0.95838,0.96643/
data (bm3i( 8, 8,ibeta ), ibeta=1,10)/ &
0.93337,0.93615,0.94220,0.94937,0.95647,0.96292,0.96840,0.97283, &
0.97619,0.97854/
data (bm3i( 8, 9,ibeta ), ibeta=1,10)/ &
0.97790,0.97891,0.98105,0.98346,0.98569,0.98751,0.98879,0.98950, &
0.98961,0.98912/
data (bm3i( 8,10,ibeta ), ibeta=1,10)/ &
0.99337,0.99367,0.99430,0.99493,0.99541,0.99562,0.99551,0.99501, &
0.99410,0.99274/
data (bm3i( 9, 1,ibeta ), ibeta=1,10)/ &
0.98470,0.98594,0.98844,0.99106,0.99334,0.99514,0.99650,0.99749, &
0.99821,0.99872/
data (bm3i( 9, 2,ibeta ), ibeta=1,10)/ &
0.98184,0.98330,0.98624,0.98934,0.99205,0.99420,0.99582,0.99701, &
0.99787,0.99848/
data (bm3i( 9, 3,ibeta ), ibeta=1,10)/ &
0.97288,0.97498,0.97927,0.98385,0.98789,0.99113,0.99360,0.99541, &
0.99673,0.99766/
data (bm3i( 9, 4,ibeta ), ibeta=1,10)/ &
0.95403,0.95741,0.96440,0.97202,0.97887,0.98444,0.98872,0.99190, &
0.99421,0.99586/
data (bm3i( 9, 5,ibeta ), ibeta=1,10)/ &
0.91845,0.92399,0.93567,0.94873,0.96076,0.97079,0.97865,0.98457, &
0.98892,0.99206/
data (bm3i( 9, 6,ibeta ), ibeta=1,10)/ &
0.86762,0.87533,0.89202,0.91148,0.93027,0.94669,0.96013,0.97062, &
0.97855,0.98441/
data (bm3i( 9, 7,ibeta ), ibeta=1,10)/ &
0.84550,0.85253,0.86816,0.88721,0.90671,0.92490,0.94083,0.95413, &
0.96481,0.97314/
data (bm3i( 9, 8,ibeta ), ibeta=1,10)/ &
0.90138,0.90544,0.91437,0.92513,0.93602,0.94615,0.95506,0.96258, &
0.96868,0.97347/
data (bm3i( 9, 9,ibeta ), ibeta=1,10)/ &
0.96248,0.96415,0.96773,0.97187,0.97583,0.97925,0.98198,0.98394, &
0.98514,0.98559/
data (bm3i( 9,10,ibeta ), ibeta=1,10)/ &
0.98837,0.98892,0.99005,0.99127,0.99232,0.99306,0.99339,0.99328, &
0.99269,0.99161/
data (bm3i(10, 1,ibeta ), ibeta=1,10)/ &
0.99080,0.99158,0.99311,0.99471,0.99607,0.99715,0.99795,0.99853, &
0.99895,0.99925/
data (bm3i(10, 2,ibeta ), ibeta=1,10)/ &
0.98910,0.99001,0.99182,0.99371,0.99533,0.99661,0.99757,0.99826, &
0.99876,0.99912/
data (bm3i(10, 3,ibeta ), ibeta=1,10)/ &
0.98374,0.98506,0.98772,0.99051,0.99294,0.99486,0.99630,0.99736, &
0.99812,0.99866/
data (bm3i(10, 4,ibeta ), ibeta=1,10)/ &
0.97238,0.97453,0.97892,0.98361,0.98773,0.99104,0.99354,0.99538, &
0.99671,0.99765/
data (bm3i(10, 5,ibeta ), ibeta=1,10)/ &
0.94961,0.95333,0.96103,0.96941,0.97693,0.98303,0.98772,0.99119, &
0.99371,0.99551/
data (bm3i(10, 6,ibeta ), ibeta=1,10)/ &
0.90943,0.91550,0.92834,0.94275,0.95608,0.96723,0.97600,0.98263, &
0.98751,0.99103/
data (bm3i(10, 7,ibeta ), ibeta=1,10)/ &
0.86454,0.87200,0.88829,0.90749,0.92630,0.94300,0.95687,0.96785, &
0.97626,0.98254/
data (bm3i(10, 8,ibeta ), ibeta=1,10)/ &
0.87498,0.88048,0.89264,0.90737,0.92240,0.93642,0.94877,0.95917, &
0.96762,0.97429/
data (bm3i(10, 9,ibeta ), ibeta=1,10)/ &
0.93946,0.94209,0.94781,0.95452,0.96111,0.96704,0.97203,0.97602, &
0.97900,0.98106/
data (bm3i(10,10,ibeta ), ibeta=1,10)/ &
0.97977,0.98071,0.98270,0.98492,0.98695,0.98858,0.98970,0.99027, &
0.99026,0.98968/
! fsb fm correction for intramodal m2 coagulation
data bm2ii / &
0.707107, 0.720583, 0.745310, 0.748056, 0.696935, &
0.604164, 0.504622, 0.416559, 0.343394, 0.283641/
! *** total correction for intramodal m2 coagulation
data bm2iitt / &
1.000000, 0.907452, 0.680931, 0.409815, 0.196425, &
0.078814, 0.028473, 0.009800, 0.003322, 0.001129/
! fsb fm correction for m2 i to j coagulation
data (bm2ij ( 1, 1,ibeta), ibeta = 1,10) / &
0.707107, 0.716828, 0.738240, 0.764827, 0.793610, &
0.822843, 0.851217, 0.877670, 0.901404, 0.921944/
data (bm2ij ( 1, 2,ibeta), ibeta = 1,10) / &
0.719180, 0.727975, 0.747638, 0.772334, 0.799234, &
0.826666, 0.853406, 0.878482, 0.901162, 0.920987/
data (bm2ij ( 1, 3,ibeta), ibeta = 1,10) / &
0.760947, 0.767874, 0.783692, 0.803890, 0.826015, &
0.848562, 0.870498, 0.891088, 0.909823, 0.926400/
data (bm2ij ( 1, 4,ibeta), ibeta = 1,10) / &
0.830926, 0.836034, 0.847708, 0.862528, 0.878521, &
0.894467, 0.909615, 0.923520, 0.935959, 0.946858/
data (bm2ij ( 1, 5,ibeta), ibeta = 1,10) / &
0.903643, 0.907035, 0.914641, 0.924017, 0.933795, &
0.943194, 0.951806, 0.959449, 0.966087, 0.971761/
data (bm2ij ( 1, 6,ibeta), ibeta = 1,10) / &
0.954216, 0.956094, 0.960211, 0.965123, 0.970068, &
0.974666, 0.978750, 0.982277, 0.985268, 0.987775/
data (bm2ij ( 1, 7,ibeta), ibeta = 1,10) / &
0.980546, 0.981433, 0.983343, 0.985568, 0.987751, &
0.989735, 0.991461, 0.992926, 0.994150, 0.995164/
data (bm2ij ( 1, 8,ibeta), ibeta = 1,10) / &
0.992142, 0.992524, 0.993338, 0.994272, 0.995174, &
0.995981, 0.996675, 0.997257, 0.997740, 0.998137/
data (bm2ij ( 1, 9,ibeta), ibeta = 1,10) / &
0.996868, 0.997026, 0.997361, 0.997742, 0.998106, &
0.998430, 0.998705, 0.998935, 0.999125, 0.999280/
data (bm2ij ( 1, 10,ibeta), ibeta = 1,10) / &
0.998737, 0.998802, 0.998939, 0.999094, 0.999241, &
0.999371, 0.999481, 0.999573, 0.999648, 0.999709/
data (bm2ij ( 2, 1,ibeta), ibeta = 1,10) / &
0.729600, 0.739948, 0.763059, 0.791817, 0.822510, &
0.852795, 0.881000, 0.905999, 0.927206, 0.944532/
data (bm2ij ( 2, 2,ibeta), ibeta = 1,10) / &
0.727025, 0.737116, 0.759615, 0.787657, 0.817740, &
0.847656, 0.875801, 0.901038, 0.922715, 0.940643/
data (bm2ij ( 2, 3,ibeta), ibeta = 1,10) / &
0.738035, 0.746779, 0.766484, 0.791340, 0.818324, &
0.845546, 0.871629, 0.895554, 0.916649, 0.934597/
data (bm2ij ( 2, 4,ibeta), ibeta = 1,10) / &
0.784185, 0.790883, 0.806132, 0.825501, 0.846545, &
0.867745, 0.888085, 0.906881, 0.923705, 0.938349/
data (bm2ij ( 2, 5,ibeta), ibeta = 1,10) / &
0.857879, 0.862591, 0.873238, 0.886539, 0.900645, &
0.914463, 0.927360, 0.939004, 0.949261, 0.958125/
data (bm2ij ( 2, 6,ibeta), ibeta = 1,10) / &
0.925441, 0.928304, 0.934645, 0.942324, 0.950181, &
0.957600, 0.964285, 0.970133, 0.975147, 0.979388/
data (bm2ij ( 2, 7,ibeta), ibeta = 1,10) / &
0.966728, 0.968176, 0.971323, 0.975027, 0.978705, &
0.982080, 0.985044, 0.987578, 0.989710, 0.991485/
data (bm2ij ( 2, 8,ibeta), ibeta = 1,10) / &
0.986335, 0.986980, 0.988362, 0.989958, 0.991511, &
0.992912, 0.994122, 0.995143, 0.995992, 0.996693/
data (bm2ij ( 2, 9,ibeta), ibeta = 1,10) / &
0.994547, 0.994817, 0.995391, 0.996046, 0.996677, &
0.997238, 0.997719, 0.998122, 0.998454, 0.998727/
data (bm2ij ( 2, 10,ibeta), ibeta = 1,10) / &
0.997817, 0.997928, 0.998163, 0.998429, 0.998683, &
0.998908, 0.999099, 0.999258, 0.999389, 0.999497/
data (bm2ij ( 3, 1,ibeta), ibeta = 1,10) / &
0.783612, 0.793055, 0.814468, 0.841073, 0.868769, &
0.894963, 0.918118, 0.937527, 0.953121, 0.965244/
data (bm2ij ( 3, 2,ibeta), ibeta = 1,10) / &
0.772083, 0.781870, 0.803911, 0.831238, 0.859802, &
0.887036, 0.911349, 0.931941, 0.948649, 0.961751/
data (bm2ij ( 3, 3,ibeta), ibeta = 1,10) / &
0.755766, 0.765509, 0.787380, 0.814630, 0.843526, &
0.871670, 0.897443, 0.919870, 0.938557, 0.953576/
data (bm2ij ( 3, 4,ibeta), ibeta = 1,10) / &
0.763816, 0.772145, 0.790997, 0.814784, 0.840434, &
0.865978, 0.890034, 0.911671, 0.930366, 0.945963/
data (bm2ij ( 3, 5,ibeta), ibeta = 1,10) / &
0.813597, 0.819809, 0.833889, 0.851618, 0.870640, &
0.889514, 0.907326, 0.923510, 0.937768, 0.950003/
data (bm2ij ( 3, 6,ibeta), ibeta = 1,10) / &
0.886317, 0.890437, 0.899643, 0.910955, 0.922730, &
0.934048, 0.944422, 0.953632, 0.961624, 0.968444/
data (bm2ij ( 3, 7,ibeta), ibeta = 1,10) / &
0.944565, 0.946855, 0.951872, 0.957854, 0.963873, &
0.969468, 0.974438, 0.978731, 0.982372, 0.985424/
data (bm2ij ( 3, 8,ibeta), ibeta = 1,10) / &
0.976358, 0.977435, 0.979759, 0.982467, 0.985125, &
0.987540, 0.989642, 0.991425, 0.992916, 0.994150/
data (bm2ij ( 3, 9,ibeta), ibeta = 1,10) / &
0.990471, 0.990932, 0.991917, 0.993048, 0.994142, &
0.995121, 0.995964, 0.996671, 0.997258, 0.997740/
data (bm2ij ( 3, 10,ibeta), ibeta = 1,10) / &
0.996199, 0.996389, 0.996794, 0.997254, 0.997694, &
0.998086, 0.998420, 0.998699, 0.998929, 0.999117/
data (bm2ij ( 4, 1,ibeta), ibeta = 1,10) / &
0.844355, 0.852251, 0.869914, 0.891330, 0.912823, &
0.932259, 0.948642, 0.961767, 0.971897, 0.979510/
data (bm2ij ( 4, 2,ibeta), ibeta = 1,10) / &
0.831550, 0.839954, 0.858754, 0.881583, 0.904592, &
0.925533, 0.943309, 0.957647, 0.968779, 0.977185/
data (bm2ij ( 4, 3,ibeta), ibeta = 1,10) / &
0.803981, 0.813288, 0.834060, 0.859400, 0.885285, &
0.909286, 0.930084, 0.947193, 0.960714, 0.971078/
data (bm2ij ( 4, 4,ibeta), ibeta = 1,10) / &
0.781787, 0.791080, 0.811931, 0.837749, 0.864768, &
0.890603, 0.913761, 0.933477, 0.949567, 0.962261/
data (bm2ij ( 4, 5,ibeta), ibeta = 1,10) / &
0.791591, 0.799355, 0.816916, 0.838961, 0.862492, &
0.885595, 0.907003, 0.925942, 0.942052, 0.955310/
data (bm2ij ( 4, 6,ibeta), ibeta = 1,10) / &
0.844933, 0.850499, 0.863022, 0.878593, 0.895038, &
0.911072, 0.925939, 0.939227, 0.950765, 0.960550/
data (bm2ij ( 4, 7,ibeta), ibeta = 1,10) / &
0.912591, 0.916022, 0.923607, 0.932777, 0.942151, &
0.951001, 0.958976, 0.965950, 0.971924, 0.976965/
data (bm2ij ( 4, 8,ibeta), ibeta = 1,10) / &
0.959859, 0.961617, 0.965433, 0.969924, 0.974382, &
0.978472, 0.982063, 0.985134, 0.987716, 0.989865/
data (bm2ij ( 4, 9,ibeta), ibeta = 1,10) / &
0.983377, 0.984162, 0.985844, 0.987788, 0.989681, &
0.991386, 0.992860, 0.994104, 0.995139, 0.995991/
data (bm2ij ( 4, 10,ibeta), ibeta = 1,10) / &
0.993343, 0.993672, 0.994370, 0.995169, 0.995937, &
0.996622, 0.997209, 0.997700, 0.998106, 0.998439/
data (bm2ij ( 5, 1,ibeta), ibeta = 1,10) / &
0.895806, 0.901918, 0.915233, 0.930783, 0.945768, &
0.958781, 0.969347, 0.977540, 0.983697, 0.988225/
data (bm2ij ( 5, 2,ibeta), ibeta = 1,10) / &
0.885634, 0.892221, 0.906629, 0.923540, 0.939918, &
0.954213, 0.965873, 0.974951, 0.981794, 0.986840/
data (bm2ij ( 5, 3,ibeta), ibeta = 1,10) / &
0.860120, 0.867858, 0.884865, 0.904996, 0.924724, &
0.942177, 0.956602, 0.967966, 0.976616, 0.983043/
data (bm2ij ( 5, 4,ibeta), ibeta = 1,10) / &
0.827462, 0.836317, 0.855885, 0.879377, 0.902897, &
0.924232, 0.942318, 0.956900, 0.968222, 0.976774/
data (bm2ij ( 5, 5,ibeta), ibeta = 1,10) / &
0.805527, 0.814279, 0.833853, 0.857892, 0.882726, &
0.906095, 0.926690, 0.943938, 0.957808, 0.968615/
data (bm2ij ( 5, 6,ibeta), ibeta = 1,10) / &
0.820143, 0.827223, 0.843166, 0.863002, 0.883905, &
0.904128, 0.922585, 0.938687, 0.952222, 0.963255/
data (bm2ij ( 5, 7,ibeta), ibeta = 1,10) / &
0.875399, 0.880208, 0.890929, 0.904065, 0.917699, &
0.930756, 0.942656, 0.953131, 0.962113, 0.969657/
data (bm2ij ( 5, 8,ibeta), ibeta = 1,10) / &
0.934782, 0.937520, 0.943515, 0.950656, 0.957840, &
0.964516, 0.970446, 0.975566, 0.979905, 0.983534/
data (bm2ij ( 5, 9,ibeta), ibeta = 1,10) / &
0.971369, 0.972679, 0.975505, 0.978797, 0.982029, &
0.984964, 0.987518, 0.989685, 0.991496, 0.992994/
data (bm2ij ( 5, 10,ibeta), ibeta = 1,10) / &
0.988329, 0.988893, 0.990099, 0.991485, 0.992825, &
0.994025, 0.995058, 0.995925, 0.996643, 0.997234/
data (bm2ij ( 6, 1,ibeta), ibeta = 1,10) / &
0.933384, 0.937784, 0.947130, 0.957655, 0.967430, &
0.975639, 0.982119, 0.987031, 0.990657, 0.993288/
data (bm2ij ( 6, 2,ibeta), ibeta = 1,10) / &
0.926445, 0.931227, 0.941426, 0.952975, 0.963754, &
0.972845, 0.980044, 0.985514, 0.989558, 0.992498/
data (bm2ij ( 6, 3,ibeta), ibeta = 1,10) / &
0.907835, 0.913621, 0.926064, 0.940308, 0.953745, &
0.965189, 0.974327, 0.981316, 0.986510, 0.990297/
data (bm2ij ( 6, 4,ibeta), ibeta = 1,10) / &
0.879088, 0.886306, 0.901945, 0.920079, 0.937460, &
0.952509, 0.964711, 0.974166, 0.981265, 0.986484/
data (bm2ij ( 6, 5,ibeta), ibeta = 1,10) / &
0.846500, 0.854862, 0.873189, 0.894891, 0.916264, &
0.935315, 0.951197, 0.963812, 0.973484, 0.980715/
data (bm2ij ( 6, 6,ibeta), ibeta = 1,10) / &
0.828137, 0.836250, 0.854310, 0.876287, 0.898710, &
0.919518, 0.937603, 0.952560, 0.964461, 0.973656/
data (bm2ij ( 6, 7,ibeta), ibeta = 1,10) / &
0.848595, 0.854886, 0.868957, 0.886262, 0.904241, &
0.921376, 0.936799, 0.950096, 0.961172, 0.970145/
data (bm2ij ( 6, 8,ibeta), ibeta = 1,10) / &
0.902919, 0.906922, 0.915760, 0.926427, 0.937312, &
0.947561, 0.956758, 0.964747, 0.971525, 0.977175/
data (bm2ij ( 6, 9,ibeta), ibeta = 1,10) / &
0.952320, 0.954434, 0.959021, 0.964418, 0.969774, &
0.974688, 0.979003, 0.982690, 0.985789, 0.988364/
data (bm2ij ( 6, 10,ibeta), ibeta = 1,10) / &
0.979689, 0.980650, 0.982712, 0.985093, 0.987413, &
0.989502, 0.991308, 0.992831, 0.994098, 0.995142/
data (bm2ij ( 7, 1,ibeta), ibeta = 1,10) / &
0.958611, 0.961598, 0.967817, 0.974620, 0.980752, &
0.985771, 0.989650, 0.992543, 0.994653, 0.996171/
data (bm2ij ( 7, 2,ibeta), ibeta = 1,10) / &
0.954225, 0.957488, 0.964305, 0.971795, 0.978576, &
0.984144, 0.988458, 0.991681, 0.994034, 0.995728/
data (bm2ij ( 7, 3,ibeta), ibeta = 1,10) / &
0.942147, 0.946158, 0.954599, 0.963967, 0.972529, &
0.979612, 0.985131, 0.989271, 0.992301, 0.994487/
data (bm2ij ( 7, 4,ibeta), ibeta = 1,10) / &
0.921821, 0.927048, 0.938140, 0.950598, 0.962118, &
0.971752, 0.979326, 0.985046, 0.989254, 0.992299/
data (bm2ij ( 7, 5,ibeta), ibeta = 1,10) / &
0.893419, 0.900158, 0.914598, 0.931070, 0.946584, &
0.959795, 0.970350, 0.978427, 0.984432, 0.988811/
data (bm2ij ( 7, 6,ibeta), ibeta = 1,10) / &
0.863302, 0.871111, 0.888103, 0.907990, 0.927305, &
0.944279, 0.958245, 0.969211, 0.977540, 0.983720/
data (bm2ij ( 7, 7,ibeta), ibeta = 1,10) / &
0.850182, 0.857560, 0.873890, 0.893568, 0.913408, &
0.931591, 0.947216, 0.960014, 0.970121, 0.977886/
data (bm2ij ( 7, 8,ibeta), ibeta = 1,10) / &
0.875837, 0.881265, 0.893310, 0.907936, 0.922910, &
0.936977, 0.949480, 0.960154, 0.968985, 0.976111/
data (bm2ij ( 7, 9,ibeta), ibeta = 1,10) / &
0.926228, 0.929445, 0.936486, 0.944868, 0.953293, &
0.961108, 0.968028, 0.973973, 0.978974, 0.983118/
data (bm2ij ( 7, 10,ibeta), ibeta = 1,10) / &
0.965533, 0.967125, 0.970558, 0.974557, 0.978484, &
0.982050, 0.985153, 0.987785, 0.989982, 0.991798/
data (bm2ij ( 8, 1,ibeta), ibeta = 1,10) / &
0.974731, 0.976674, 0.980660, 0.984926, 0.988689, &
0.991710, 0.994009, 0.995703, 0.996929, 0.997805/
data (bm2ij ( 8, 2,ibeta), ibeta = 1,10) / &
0.972062, 0.974192, 0.978571, 0.983273, 0.987432, &
0.990780, 0.993333, 0.995218, 0.996581, 0.997557/
data (bm2ij ( 8, 3,ibeta), ibeta = 1,10) / &
0.964662, 0.967300, 0.972755, 0.978659, 0.983921, &
0.988181, 0.991444, 0.993859, 0.995610, 0.996863/
data (bm2ij ( 8, 4,ibeta), ibeta = 1,10) / &
0.951782, 0.955284, 0.962581, 0.970559, 0.977737, &
0.983593, 0.988103, 0.991454, 0.993889, 0.995635/
data (bm2ij ( 8, 5,ibeta), ibeta = 1,10) / &
0.931947, 0.936723, 0.946751, 0.957843, 0.967942, &
0.976267, 0.982734, 0.987571, 0.991102, 0.993642/
data (bm2ij ( 8, 6,ibeta), ibeta = 1,10) / &
0.905410, 0.911665, 0.924950, 0.939908, 0.953798, &
0.965469, 0.974684, 0.981669, 0.986821, 0.990556/
data (bm2ij ( 8, 7,ibeta), ibeta = 1,10) / &
0.878941, 0.886132, 0.901679, 0.919688, 0.936970, &
0.951980, 0.964199, 0.973709, 0.980881, 0.986174/
data (bm2ij ( 8, 8,ibeta), ibeta = 1,10) / &
0.871653, 0.878218, 0.892652, 0.909871, 0.927034, &
0.942592, 0.955836, 0.966604, 0.975065, 0.981545/
data (bm2ij ( 8, 9,ibeta), ibeta = 1,10) / &
0.900693, 0.905239, 0.915242, 0.927232, 0.939335, &
0.950555, 0.960420, 0.968774, 0.975651, 0.981188/
data (bm2ij ( 8, 10,ibeta), ibeta = 1,10) / &
0.944922, 0.947435, 0.952894, 0.959317, 0.965689, &
0.971529, 0.976645, 0.981001, 0.984641, 0.987642/
data (bm2ij ( 9, 1,ibeta), ibeta = 1,10) / &
0.984736, 0.985963, 0.988453, 0.991078, 0.993357, &
0.995161, 0.996519, 0.997512, 0.998226, 0.998734/
data (bm2ij ( 9, 2,ibeta), ibeta = 1,10) / &
0.983141, 0.984488, 0.987227, 0.990119, 0.992636, &
0.994632, 0.996137, 0.997238, 0.998030, 0.998595/
data (bm2ij ( 9, 3,ibeta), ibeta = 1,10) / &
0.978726, 0.980401, 0.983819, 0.987450, 0.990626, &
0.993157, 0.995071, 0.996475, 0.997486, 0.998206/
data (bm2ij ( 9, 4,ibeta), ibeta = 1,10) / &
0.970986, 0.973224, 0.977818, 0.982737, 0.987072, &
0.990546, 0.993184, 0.995124, 0.996523, 0.997521/
data (bm2ij ( 9, 5,ibeta), ibeta = 1,10) / &
0.958579, 0.961700, 0.968149, 0.975116, 0.981307, &
0.986301, 0.990112, 0.992923, 0.994954, 0.996404/
data (bm2ij ( 9, 6,ibeta), ibeta = 1,10) / &
0.940111, 0.944479, 0.953572, 0.963506, 0.972436, &
0.979714, 0.985313, 0.989468, 0.992483, 0.994641/
data (bm2ij ( 9, 7,ibeta), ibeta = 1,10) / &
0.916127, 0.921878, 0.934003, 0.947506, 0.959899, &
0.970199, 0.978255, 0.984314, 0.988755, 0.991960/
data (bm2ij ( 9, 8,ibeta), ibeta = 1,10) / &
0.893848, 0.900364, 0.914368, 0.930438, 0.945700, &
0.958824, 0.969416, 0.977603, 0.983746, 0.988262/
data (bm2ij ( 9, 9,ibeta), ibeta = 1,10) / &
0.892161, 0.897863, 0.910315, 0.925021, 0.939523, &
0.952544, 0.963544, 0.972442, 0.979411, 0.984742/
data (bm2ij ( 9, 10,ibeta), ibeta = 1,10) / &
0.922260, 0.925966, 0.934047, 0.943616, 0.953152, &
0.961893, 0.969506, 0.975912, 0.981167, 0.985394/
data (bm2ij ( 10, 1,ibeta), ibeta = 1,10) / &
0.990838, 0.991598, 0.993128, 0.994723, 0.996092, &
0.997167, 0.997969, 0.998552, 0.998969, 0.999265/
data (bm2ij ( 10, 2,ibeta), ibeta = 1,10) / &
0.989892, 0.990727, 0.992411, 0.994167, 0.995678, &
0.996864, 0.997751, 0.998396, 0.998858, 0.999186/
data (bm2ij ( 10, 3,ibeta), ibeta = 1,10) / &
0.987287, 0.988327, 0.990428, 0.992629, 0.994529, &
0.996026, 0.997148, 0.997965, 0.998551, 0.998967/
data (bm2ij ( 10, 4,ibeta), ibeta = 1,10) / &
0.982740, 0.984130, 0.986952, 0.989926, 0.992508, &
0.994551, 0.996087, 0.997208, 0.998012, 0.998584/
data (bm2ij ( 10, 5,ibeta), ibeta = 1,10) / &
0.975380, 0.977330, 0.981307, 0.985529, 0.989216, &
0.992147, 0.994358, 0.995975, 0.997136, 0.997961/
data (bm2ij ( 10, 6,ibeta), ibeta = 1,10) / &
0.963911, 0.966714, 0.972465, 0.978614, 0.984022, &
0.988346, 0.991620, 0.994020, 0.995747, 0.996974/
data (bm2ij ( 10, 7,ibeta), ibeta = 1,10) / &
0.947187, 0.951161, 0.959375, 0.968258, 0.976160, &
0.982540, 0.987409, 0.991000, 0.993592, 0.995441/
data (bm2ij ( 10, 8,ibeta), ibeta = 1,10) / &
0.926045, 0.931270, 0.942218, 0.954297, 0.965273, &
0.974311, 0.981326, 0.986569, 0.990394, 0.993143/
data (bm2ij ( 10, 9,ibeta), ibeta = 1,10) / &
0.908092, 0.913891, 0.926288, 0.940393, 0.953667, &
0.964987, 0.974061, 0.981038, 0.986253, 0.990078/
data (bm2ij ( 10, 10,ibeta), ibeta = 1,10) / &
0.911143, 0.915972, 0.926455, 0.938721, 0.950701, &
0.961370, 0.970329, 0.977549, 0.983197, 0.987518/
! fsb total correction factor for m2 coagulation j from i
data (bm2ji( 1, 1,ibeta), ibeta = 1,10) / &
0.753466, 0.756888, 0.761008, 0.759432, 0.748675, &
0.726951, 0.693964, 0.650915, 0.600227, 0.545000/
data (bm2ji( 1, 2,ibeta), ibeta = 1,10) / &
0.824078, 0.828698, 0.835988, 0.838943, 0.833454, &
0.817148, 0.789149, 0.750088, 0.701887, 0.647308/
data (bm2ji( 1, 3,ibeta), ibeta = 1,10) / &
1.007389, 1.014362, 1.028151, 1.041011, 1.047939, &
1.045707, 1.032524, 1.007903, 0.972463, 0.927667/
data (bm2ji( 1, 4,ibeta), ibeta = 1,10) / &
1.246157, 1.255135, 1.274249, 1.295351, 1.313362, &
1.325187, 1.329136, 1.324491, 1.311164, 1.289459/
data (bm2ji( 1, 5,ibeta), ibeta = 1,10) / &
1.450823, 1.459551, 1.478182, 1.499143, 1.518224, &
1.533312, 1.543577, 1.548882, 1.549395, 1.545364/
data (bm2ji( 1, 6,ibeta), ibeta = 1,10) / &
1.575248, 1.581832, 1.595643, 1.610866, 1.624601, &
1.635690, 1.643913, 1.649470, 1.652688, 1.653878/
data (bm2ji( 1, 7,ibeta), ibeta = 1,10) / &
1.638426, 1.642626, 1.651293, 1.660641, 1.668926, &
1.675571, 1.680572, 1.684147, 1.686561, 1.688047/
data (bm2ji( 1, 8,ibeta), ibeta = 1,10) / &
1.669996, 1.672392, 1.677283, 1.682480, 1.687028, &
1.690651, 1.693384, 1.695372, 1.696776, 1.697734/
data (bm2ji( 1, 9,ibeta), ibeta = 1,10) / &
1.686148, 1.687419, 1.689993, 1.692704, 1.695057, &
1.696922, 1.698329, 1.699359, 1.700099, 1.700621/
data (bm2ji( 1,10,ibeta), ibeta = 1,10) / &
1.694364, 1.695010, 1.696313, 1.697676, 1.698853, &
1.699782, 1.700482, 1.700996, 1.701366, 1.701631/
data (bm2ji( 2, 1,ibeta), ibeta = 1,10) / &
0.783166, 0.779369, 0.768044, 0.747572, 0.716709, &
0.675422, 0.624981, 0.567811, 0.507057, 0.445975/
data (bm2ji( 2, 2,ibeta), ibeta = 1,10) / &
0.848390, 0.847100, 0.840874, 0.826065, 0.800296, &
0.762625, 0.713655, 0.655545, 0.591603, 0.525571/
data (bm2ji( 2, 3,ibeta), ibeta = 1,10) / &
1.039894, 1.043786, 1.049445, 1.049664, 1.039407, &
1.015322, 0.975983, 0.922180, 0.856713, 0.783634/
data (bm2ji( 2, 4,ibeta), ibeta = 1,10) / &
1.345995, 1.356064, 1.376947, 1.398304, 1.412685, &
1.414611, 1.400652, 1.369595, 1.322261, 1.260993/
data (bm2ji( 2, 5,ibeta), ibeta = 1,10) / &
1.675575, 1.689859, 1.720957, 1.756659, 1.788976, &
1.812679, 1.824773, 1.824024, 1.810412, 1.784630/
data (bm2ji( 2, 6,ibeta), ibeta = 1,10) / &
1.919835, 1.933483, 1.962973, 1.996810, 2.028377, &
2.054172, 2.072763, 2.083963, 2.088190, 2.086052/
data (bm2ji( 2, 7,ibeta), ibeta = 1,10) / &
2.064139, 2.074105, 2.095233, 2.118909, 2.140688, &
2.158661, 2.172373, 2.182087, 2.188330, 2.191650/
data (bm2ji( 2, 8,ibeta), ibeta = 1,10) / &
2.144871, 2.150990, 2.163748, 2.177731, 2.190364, &
2.200712, 2.208687, 2.214563, 2.218716, 2.221502/
data (bm2ji( 2, 9,ibeta), ibeta = 1,10) / &
2.189223, 2.192595, 2.199540, 2.207033, 2.213706, &
2.219125, 2.223297, 2.226403, 2.228660, 2.230265/
data (bm2ji( 2,10,ibeta), ibeta = 1,10) / &
2.212595, 2.214342, 2.217912, 2.221723, 2.225082, &
2.227791, 2.229869, 2.231417, 2.232551, 2.233372/
data (bm2ji( 3, 1,ibeta), ibeta = 1,10) / &
0.837870, 0.824476, 0.793119, 0.750739, 0.700950, &
0.646691, 0.590508, 0.534354, 0.479532, 0.426856/
data (bm2ji( 3, 2,ibeta), ibeta = 1,10) / &
0.896771, 0.885847, 0.859327, 0.821694, 0.775312, &
0.722402, 0.665196, 0.605731, 0.545742, 0.486687/
data (bm2ji( 3, 3,ibeta), ibeta = 1,10) / &
1.076089, 1.071727, 1.058845, 1.036171, 1.002539, &
0.957521, 0.901640, 0.836481, 0.764597, 0.689151/
data (bm2ji( 3, 4,ibeta), ibeta = 1,10) / &
1.409571, 1.415168, 1.425346, 1.432021, 1.428632, &
1.409696, 1.371485, 1.312958, 1.236092, 1.145293/
data (bm2ji( 3, 5,ibeta), ibeta = 1,10) / &
1.862757, 1.880031, 1.918394, 1.963456, 2.004070, &
2.030730, 2.036144, 2.016159, 1.970059, 1.900079/
data (bm2ji( 3, 6,ibeta), ibeta = 1,10) / &
2.289741, 2.313465, 2.366789, 2.431612, 2.495597, &
2.549838, 2.588523, 2.608665, 2.609488, 2.591662/
data (bm2ji( 3, 7,ibeta), ibeta = 1,10) / &
2.597157, 2.618731, 2.666255, 2.722597, 2.777531, &
2.825187, 2.862794, 2.889648, 2.906199, 2.913380/
data (bm2ji( 3, 8,ibeta), ibeta = 1,10) / &
2.797975, 2.813116, 2.845666, 2.882976, 2.918289, &
2.948461, 2.972524, 2.990687, 3.003664, 3.012284/
data (bm2ji( 3, 9,ibeta), ibeta = 1,10) / &
2.920832, 2.929843, 2.948848, 2.970057, 2.989632, &
3.006057, 3.019067, 3.028979, 3.036307, 3.041574/
data (bm2ji( 3,10,ibeta), ibeta = 1,10) / &
2.989627, 2.994491, 3.004620, 3.015720, 3.025789, &
3.034121, 3.040664, 3.045641, 3.049347, 3.052066/
data (bm2ji( 4, 1,ibeta), ibeta = 1,10) / &
0.893179, 0.870897, 0.820996, 0.759486, 0.695488, &
0.634582, 0.579818, 0.532143, 0.490927, 0.454618/
data (bm2ji( 4, 2,ibeta), ibeta = 1,10) / &
0.948355, 0.927427, 0.880215, 0.821146, 0.758524, &
0.697680, 0.641689, 0.591605, 0.546919, 0.506208/
data (bm2ji( 4, 3,ibeta), ibeta = 1,10) / &
1.109562, 1.093648, 1.056438, 1.007310, 0.951960, &
0.894453, 0.837364, 0.781742, 0.727415, 0.673614/
data (bm2ji( 4, 4,ibeta), ibeta = 1,10) / &
1.423321, 1.417557, 1.402442, 1.379079, 1.347687, &
1.308075, 1.259703, 1.201983, 1.134778, 1.058878/
data (bm2ji( 4, 5,ibeta), ibeta = 1,10) / &
1.933434, 1.944347, 1.968765, 1.997653, 2.023054, &
2.036554, 2.029949, 1.996982, 1.934982, 1.845473/
data (bm2ji( 4, 6,ibeta), ibeta = 1,10) / &
2.547772, 2.577105, 2.645918, 2.735407, 2.830691, &
2.917268, 2.981724, 3.013684, 3.007302, 2.961560/
data (bm2ji( 4, 7,ibeta), ibeta = 1,10) / &
3.101817, 3.139271, 3.225851, 3.336402, 3.453409, &
3.563116, 3.655406, 3.724014, 3.766113, 3.781394/
data (bm2ji( 4, 8,ibeta), ibeta = 1,10) / &
3.540920, 3.573780, 3.647439, 3.737365, 3.828468, &
3.911436, 3.981317, 4.036345, 4.076749, 4.103751/
data (bm2ji( 4, 9,ibeta), ibeta = 1,10) / &
3.856771, 3.879363, 3.928579, 3.986207, 4.042173, &
4.091411, 4.132041, 4.164052, 4.188343, 4.206118/
data (bm2ji( 4,10,ibeta), ibeta = 1,10) / &
4.053923, 4.067191, 4.095509, 4.127698, 4.158037, &
4.184055, 4.205135, 4.221592, 4.234115, 4.243463/
data (bm2ji( 5, 1,ibeta), ibeta = 1,10) / &
0.935846, 0.906814, 0.843358, 0.768710, 0.695885, &
0.631742, 0.579166, 0.538471, 0.508410, 0.486863/
data (bm2ji( 5, 2,ibeta), ibeta = 1,10) / &
0.988308, 0.959524, 0.896482, 0.821986, 0.748887, &
0.684168, 0.630908, 0.589516, 0.558676, 0.536056/
data (bm2ji( 5, 3,ibeta), ibeta = 1,10) / &
1.133795, 1.107139, 1.048168, 0.977258, 0.906341, &
0.842477, 0.789093, 0.746731, 0.713822, 0.687495/
data (bm2ji( 5, 4,ibeta), ibeta = 1,10) / &
1.405692, 1.385781, 1.340706, 1.284776, 1.227085, &
1.173532, 1.127008, 1.087509, 1.052712, 1.018960/
data (bm2ji( 5, 5,ibeta), ibeta = 1,10) / &
1.884992, 1.879859, 1.868463, 1.854995, 1.841946, &
1.829867, 1.816972, 1.799319, 1.771754, 1.729406/
data (bm2ji( 5, 6,ibeta), ibeta = 1,10) / &
2.592275, 2.612268, 2.661698, 2.731803, 2.815139, &
2.901659, 2.978389, 3.031259, 3.048045, 3.021122/
data (bm2ji( 5, 7,ibeta), ibeta = 1,10) / &
3.390321, 3.435519, 3.545615, 3.698419, 3.876958, &
4.062790, 4.236125, 4.378488, 4.475619, 4.519170/
data (bm2ji( 5, 8,ibeta), ibeta = 1,10) / &
4.161376, 4.216558, 4.346896, 4.519451, 4.711107, &
4.902416, 5.077701, 5.226048, 5.341423, 5.421764/
data (bm2ji( 5, 9,ibeta), ibeta = 1,10) / &
4.843961, 4.892035, 5.001492, 5.138515, 5.281684, &
5.416805, 5.535493, 5.634050, 5.712063, 5.770996/
data (bm2ji( 5,10,ibeta), ibeta = 1,10) / &
5.352093, 5.385119, 5.458056, 5.545311, 5.632162, &
5.710566, 5.777005, 5.830863, 5.873123, 5.905442/
data (bm2ji( 6, 1,ibeta), ibeta = 1,10) / &
0.964038, 0.930794, 0.859433, 0.777776, 0.700566, &
0.634671, 0.582396, 0.543656, 0.517284, 0.501694/
data (bm2ji( 6, 2,ibeta), ibeta = 1,10) / &
1.013416, 0.979685, 0.907197, 0.824135, 0.745552, &
0.678616, 0.625870, 0.587348, 0.561864, 0.547674/
data (bm2ji( 6, 3,ibeta), ibeta = 1,10) / &
1.145452, 1.111457, 1.038152, 0.953750, 0.873724, &
0.805955, 0.753621, 0.717052, 0.694920, 0.684910/
data (bm2ji( 6, 4,ibeta), ibeta = 1,10) / &
1.376547, 1.345004, 1.276415, 1.196704, 1.121091, &
1.058249, 1.012197, 0.983522, 0.970323, 0.968933/
data (bm2ji( 6, 5,ibeta), ibeta = 1,10) / &
1.778801, 1.755897, 1.706074, 1.649008, 1.597602, &
1.560087, 1.540365, 1.538205, 1.549738, 1.568333/
data (bm2ji( 6, 6,ibeta), ibeta = 1,10) / &
2.447603, 2.445172, 2.443762, 2.451842, 2.475877, &
2.519039, 2.580118, 2.653004, 2.727234, 2.789738/
data (bm2ji( 6, 7,ibeta), ibeta = 1,10) / &
3.368490, 3.399821, 3.481357, 3.606716, 3.772101, &
3.969416, 4.184167, 4.396163, 4.582502, 4.721838/
data (bm2ji( 6, 8,ibeta), ibeta = 1,10) / &
4.426458, 4.489861, 4.648250, 4.877510, 5.160698, &
5.477495, 5.803123, 6.111250, 6.378153, 6.586050/
data (bm2ji( 6, 9,ibeta), ibeta = 1,10) / &
5.568061, 5.644988, 5.829837, 6.081532, 6.371214, &
6.672902, 6.963737, 7.226172, 7.449199, 7.627886/
data (bm2ji( 6,10,ibeta), ibeta = 1,10) / &
6.639152, 6.707020, 6.863974, 7.065285, 7.281744, &
7.492437, 7.683587, 7.847917, 7.983296, 8.090977/
data (bm2ji( 7, 1,ibeta), ibeta = 1,10) / &
0.980853, 0.945724, 0.871244, 0.787311, 0.708818, &
0.641987, 0.588462, 0.547823, 0.518976, 0.500801/
data (bm2ji( 7, 2,ibeta), ibeta = 1,10) / &
1.026738, 0.990726, 0.914306, 0.828140, 0.747637, &
0.679351, 0.625127, 0.584662, 0.556910, 0.540749/
data (bm2ji( 7, 3,ibeta), ibeta = 1,10) / &
1.146496, 1.108808, 1.028695, 0.938291, 0.854101, &
0.783521, 0.728985, 0.690539, 0.667272, 0.657977/
data (bm2ji( 7, 4,ibeta), ibeta = 1,10) / &
1.344846, 1.306434, 1.224543, 1.132031, 1.046571, &
0.976882, 0.926488, 0.896067, 0.884808, 0.891027/
data (bm2ji( 7, 5,ibeta), ibeta = 1,10) / &
1.670227, 1.634583, 1.558421, 1.472939, 1.396496, &
1.339523, 1.307151, 1.300882, 1.319622, 1.360166/
data (bm2ji( 7, 6,ibeta), ibeta = 1,10) / &
2.224548, 2.199698, 2.148284, 2.095736, 2.059319, &
2.050496, 2.075654, 2.136382, 2.229641, 2.347958/
data (bm2ji( 7, 7,ibeta), ibeta = 1,10) / &
3.104483, 3.105947, 3.118398, 3.155809, 3.230427, &
3.350585, 3.519071, 3.731744, 3.976847, 4.235616/
data (bm2ji( 7, 8,ibeta), ibeta = 1,10) / &
4.288426, 4.331456, 4.447024, 4.633023, 4.891991, &
5.221458, 5.610060, 6.036467, 6.471113, 6.880462/
data (bm2ji( 7, 9,ibeta), ibeta = 1,10) / &
5.753934, 5.837061, 6.048530, 6.363800, 6.768061, &
7.241280, 7.755346, 8.276666, 8.771411, 9.210826/
data (bm2ji( 7,10,ibeta), ibeta = 1,10) / &
7.466219, 7.568810, 7.819032, 8.168340, 8.582973, &
9.030174, 9.478159, 9.899834, 10.275940, 10.595910/
data (bm2ji( 8, 1,ibeta), ibeta = 1,10) / &
0.990036, 0.954782, 0.880531, 0.797334, 0.719410, &
0.652220, 0.596923, 0.552910, 0.519101, 0.494529/
data (bm2ji( 8, 2,ibeta), ibeta = 1,10) / &
1.032428, 0.996125, 0.919613, 0.833853, 0.753611, &
0.684644, 0.628260, 0.583924, 0.550611, 0.527407/
data (bm2ji( 8, 3,ibeta), ibeta = 1,10) / &
1.141145, 1.102521, 1.021017, 0.929667, 0.844515, &
0.772075, 0.714086, 0.670280, 0.639824, 0.621970/
data (bm2ji( 8, 4,ibeta), ibeta = 1,10) / &
1.314164, 1.273087, 1.186318, 1.089208, 0.999476, &
0.924856, 0.867948, 0.829085, 0.807854, 0.803759/
data (bm2ji( 8, 5,ibeta), ibeta = 1,10) / &
1.580611, 1.538518, 1.449529, 1.350459, 1.260910, &
1.190526, 1.143502, 1.121328, 1.124274, 1.151974/
data (bm2ji( 8, 6,ibeta), ibeta = 1,10) / &
2.016773, 1.977721, 1.895727, 1.806974, 1.732891, &
1.685937, 1.673026, 1.697656, 1.761039, 1.862391/
data (bm2ji( 8, 7,ibeta), ibeta = 1,10) / &
2.750093, 2.723940, 2.672854, 2.628264, 2.612250, &
2.640406, 2.723211, 2.866599, 3.071893, 3.335217/
data (bm2ji( 8, 8,ibeta), ibeta = 1,10) / &
3.881905, 3.887143, 3.913667, 3.981912, 4.111099, &
4.316575, 4.608146, 4.988157, 5.449592, 5.974848/
data (bm2ji( 8, 9,ibeta), ibeta = 1,10) / &
5.438870, 5.492742, 5.640910, 5.886999, 6.241641, &
6.710609, 7.289480, 7.960725, 8.693495, 9.446644/
data (bm2ji( 8,10,ibeta), ibeta = 1,10) / &
7.521152, 7.624621, 7.892039, 8.300444, 8.839787, &
9.493227, 10.231770, 11.015642, 11.799990, 12.542260/
data (bm2ji( 9, 1,ibeta), ibeta = 1,10) / &
0.994285, 0.960012, 0.887939, 0.807040, 0.730578, &
0.663410, 0.606466, 0.559137, 0.520426, 0.489429/
data (bm2ji( 9, 2,ibeta), ibeta = 1,10) / &
1.033505, 0.998153, 0.923772, 0.840261, 0.761383, &
0.692242, 0.633873, 0.585709, 0.546777, 0.516215/
data (bm2ji( 9, 3,ibeta), ibeta = 1,10) / &
1.132774, 1.094907, 1.015161, 0.925627, 0.841293, &
0.767888, 0.706741, 0.657439, 0.619135, 0.591119/
data (bm2ji( 9, 4,ibeta), ibeta = 1,10) / &
1.286308, 1.245273, 1.158809, 1.061889, 0.971208, &
0.893476, 0.830599, 0.782561, 0.748870, 0.729198/
data (bm2ji( 9, 5,ibeta), ibeta = 1,10) / &
1.511105, 1.467141, 1.374520, 1.271162, 1.175871, &
1.096887, 1.037243, 0.997820, 0.978924, 0.980962/
data (bm2ji( 9, 6,ibeta), ibeta = 1,10) / &
1.857468, 1.812177, 1.717002, 1.612197, 1.519171, &
1.448660, 1.405871, 1.393541, 1.413549, 1.467532/
data (bm2ji( 9, 7,ibeta), ibeta = 1,10) / &
2.430619, 2.388452, 2.301326, 2.210241, 2.139724, &
2.104571, 2.114085, 2.174696, 2.291294, 2.467500/
data (bm2ji( 9, 8,ibeta), ibeta = 1,10) / &
3.385332, 3.357690, 3.306611, 3.269804, 3.274462, &
3.340862, 3.484609, 3.717740, 4.048748, 4.481588/
data (bm2ji( 9, 9,ibeta), ibeta = 1,10) / &
4.850497, 4.858280, 4.896008, 4.991467, 5.171511, &
5.459421, 5.873700, 6.426128, 7.119061, 7.942603/
data (bm2ji( 9,10,ibeta), ibeta = 1,10) / &
6.957098, 7.020164, 7.197272, 7.499331, 7.946554, &
8.555048, 9.330503, 10.263610, 11.327454, 12.478332/
data (bm2ji(10, 1,ibeta), ibeta = 1,10) / &
0.994567, 0.961842, 0.892854, 0.814874, 0.740198, &
0.673303, 0.615105, 0.565139, 0.522558, 0.486556/
data (bm2ji(10, 2,ibeta), ibeta = 1,10) / &
1.031058, 0.997292, 0.926082, 0.845571, 0.768501, &
0.699549, 0.639710, 0.588538, 0.545197, 0.508894/
data (bm2ji(10, 3,ibeta), ibeta = 1,10) / &
1.122535, 1.086287, 1.009790, 0.923292, 0.840626, &
0.766982, 0.703562, 0.650004, 0.605525, 0.569411/
data (bm2ji(10, 4,ibeta), ibeta = 1,10) / &
1.261142, 1.221555, 1.137979, 1.043576, 0.953745, &
0.874456, 0.807292, 0.752109, 0.708326, 0.675477/
data (bm2ji(10, 5,ibeta), ibeta = 1,10) / &
1.456711, 1.413432, 1.322096, 1.219264, 1.122319, &
1.038381, 0.969743, 0.916811, 0.879544, 0.858099/
data (bm2ji(10, 6,ibeta), ibeta = 1,10) / &
1.741792, 1.695157, 1.596897, 1.487124, 1.385734, &
1.301670, 1.238638, 1.198284, 1.181809, 1.190689/
data (bm2ji(10, 7,ibeta), ibeta = 1,10) / &
2.190197, 2.141721, 2.040226, 1.929245, 1.832051, &
1.760702, 1.721723, 1.719436, 1.757705, 1.840677/
data (bm2ji(10, 8,ibeta), ibeta = 1,10) / &
2.940764, 2.895085, 2.801873, 2.707112, 2.638603, &
2.613764, 2.644686, 2.741255, 2.912790, 3.168519/
data (bm2ji(10, 9,ibeta), ibeta = 1,10) / &
4.186191, 4.155844, 4.101953, 4.069102, 4.089886, &
4.189530, 4.389145, 4.707528, 5.161567, 5.765283/
data (bm2ji(10,10,ibeta), ibeta = 1,10) / &
6.119526, 6.127611, 6.171174, 6.286528, 6.508738, &
6.869521, 7.396912, 8.113749, 9.034683, 10.162190/
! *** end of data statements.
! *** start calculations:
constii = abs( half * ( two ) ** two3rds - one )
sqrttwo = sqrt(two)
dlgsqt2 = one / log( sqrttwo )
esat01 = exp( 0.125d0 * xxlsgat * xxlsgat )
esac01 = exp( 0.125d0 * xxlsgac * xxlsgac )
esat04 = esat01 ** 4
esac04 = esac01 ** 4
esat05 = esat04 * esat01
esac05 = esac04 * esac01
esat08 = esat04 * esat04
esac08 = esac04 * esac04
esat09 = esat08 * esat01
esac09 = esac08 * esac01
esat16 = esat08 * esat08
esac16 = esac08 * esac08
esat20 = esat16 * esat04
esac20 = esac16 * esac04
esat24 = esat20 * esat04
esac24 = esac20 * esac04
esat25 = esat20 * esat05
esac25 = esac20 * esac05
esat36 = esat20 * esat16
esac36 = esac20 * esac16
esat49 = esat24 * esat25
esat64 = esat20 * esat20 * esat24
esac64 = esac20 * esac20 * esac24
esat100 = esat64 * esat36
dgat2 = dgatk * dgatk
dgat3 = dgatk * dgatk * dgatk
dgac2 = dgacc * dgacc
dgac3 = dgacc * dgacc * dgacc
sqdgat = sqrt( dgatk )
sqdgac = sqrt( dgacc )
sqdgat5 = dgat2 * sqdgat
sqdgac5 = dgac2 * sqdgac
sqdgat7 = dgat3 * sqdgat
xm2at = dgat2 * esat16
xm3at = dgat3 * esat36
xm2ac = dgac2 * esac16
xm3ac = dgac3 * esac36
! *** for the free molecular regime: page h.3 of whitby et al. (1991)
r = sqdgac / sqdgat
r2 = r * r
r3 = r2 * r
rx4 = r2 * r2
r5 = r3 * r2
r6 = r3 * r3
rx8 = rx4 * rx4
ri1 = one / r
ri2 = one / r2
ri3 = one / r3
ri4 = ri2 * ri2
kngat = two * lamda / dgatk
kngac = two * lamda / dgacc
! *** calculate ratio of geometric mean diameters
rat = dgacc / dgatk
! *** trap subscripts for bm0 and bm0i, between 1 and 10
! see page h.5 of whitby et al. (1991)
n2n = max( 1, min( 10, &
nint( 4.0 * ( sgatk - 0.75d0 ) ) ) )
n2a = max( 1, min( 10, &
nint( 4.0 * ( sgacc - 0.75d0 ) ) ) )
n1 = max( 1, min( 10, &
1 + nint( dlgsqt2 * log( rat ) ) ) )
! *** intermodal coagulation
! *** set up for zeroeth moment
! *** near-continuum form: equation h.10a of whitby et al. (1991)
coagnc0 = knc * ( &
two + a * ( kngat * ( esat04 + r2 * esat16 * esac04 ) &
+ kngac * ( esac04 + ri2 * esac16 * esat04 ) ) &
+ ( r2 + ri2 ) * esat04 * esac04 )
! *** free-molecular form: equation h.7a of whitby et al. (1991)
coagfm0 = kfmatac * sqdgat * bm0ij(n1,n2n,n2a) * ( &
esat01 + r * esac01 + two * r2 * esat01 * esac04 &
+ rx4 * esat09 * esac16 + ri3 * esat16 * esac09 &
+ two * ri1 * esat04 + esac01 )
! *** loss to accumulation mode
! *** harmonic mean
coagatac0 = coagnc0 * coagfm0 / ( coagnc0 + coagfm0 )
qn12 = coagatac0
! *** set up for second moment
! the second moment equations are new and begin with equations a1
! through a4 of binkowski and shankar (1995). after some algebraic
! rearrangement and application of the extended mean value theorem
! of integral calculus, equations are obtained that can be solved
! analytically with correction factors as has been done by
! whitby et al. (1991)
! *** the term ( dp1 + dp2 ) ** (2/3) in equations a3 and a4 of
! binkowski and shankar (1995) is approximated by
! (dgat ** 3 + dgac **3 ) ** 2/3
! *** near-continuum form
i1nc = knc * dgat2 * ( &
two * esat16 &
+ r2 * esat04 * esac04 &
+ ri2 * esat36 * esac04 &
+ a * kngat * ( &
esat04 &
+ ri2 * esat16 * esac04 &
+ ri4 * esat36 * esac16 &
+ r2 * esac04 ) )
! *** free-molecular form
i1fm = kfmatac * sqdgat5 * bm2ij(n1,n2n,n2a) * ( &
esat25 &
+ two * r2 * esat09 * esac04 &
+ rx4 * esat01 * esac16 &
+ ri3 * esat64 * esac09 &
+ two * ri1 * esat36 * esac01 &
+ r * esat16 * esac01 )
! *** loss to accumulation mode
! *** harmonic mean
i1 = ( i1fm * i1nc ) / ( i1fm + i1nc )
coagatac2 = i1
qs12 = coagatac2
! *** gain by accumulation mode
coagacat2 = ( ( one + r6 ) ** two3rds - rx4 ) * i1
qs21 = coagacat2 * bm2ji(n1,n2n,n2a)
! *** set up for third moment
! *** near-continuum form: equation h.10b of whitby et al. (1991)
coagnc3 = knc * dgat3 * ( &
two * esat36 &
+ a * kngat * ( esat16 + r2 * esat04 * esac04 ) &
+ a * kngac * ( esat36 * esac04 + ri2 * esat64 * esac16 ) &
+ r2 * esat16 * esac04 + ri2 * esat64 * esac04 )
! *** free_molecular form: equation h.7b of whitby et al. (1991)
coagfm3 = kfmatac * sqdgat7 * bm3i( n1, n2n, n2a ) * ( &
esat49 &
+ r * esat36 * esac01 &
+ two * r2 * esat25 * esac04 &
+ rx4 * esat09 * esac16 &
+ ri3 * esat100 * esac09 &
+ two * ri1 * esat64 * esac01 )
! *** gain by accumulation mode = loss from aitken mode
! *** harmonic mean
coagatac3 = coagnc3 * coagfm3 / ( coagnc3 + coagfm3 )
qv12 = coagatac3
! *** intramodal coagulation
! *** zeroeth moment
! *** aitken mode
! *** near-continuum form: equation h.12a of whitby et al. (1991)
coagnc_at = knc * (one + esat08 + a * kngat * (esat20 + esat04))
! *** free-molecular form: equation h.11a of whitby et al. (1991)
coagfm_at = kfmat * sqdgat * bm0(n2n) * &
( esat01 + esat25 + two * esat05 )
! *** harmonic mean
coagatat0 = coagfm_at * coagnc_at / ( coagfm_at + coagnc_at )
qn11 = coagatat0
! *** accumulation mode
! *** near-continuum form: equation h.12a of whitby et al. (1991)
coagnc_ac = knc * (one + esac08 + a * kngac * (esac20 + esac04))
! *** free-molecular form: equation h.11a of whitby et al. (1991)
coagfm_ac = kfmac * sqdgac * bm0(n2a) * &
( esac01 + esac25 + two * esac05 )
! *** harmonic mean
coagacac0 = coagfm_ac * coagnc_ac / ( coagfm_ac + coagnc_ac )
qn22 = coagacac0
! *** set up for second moment
! the second moment equations are new and begin with 3.11a on page
! 45 of whitby et al. (1991). after some algebraic rearrangement and
! application of the extended mean value theorem of integral calculus
! equations are obtained that can be solved analytically with
! correction factors as has been done by whitby et al. (1991)
! *** aitken mode
! *** near-continuum
i1nc_at = knc * dgat2 * ( &
two * esat16 &
+ esat04 * esat04 &
+ esat36 * esat04 &
+ a * kngat * ( &
two * esat04 &
+ esat16 * esat04 &
+ esat36 * esat16 ) )
! *** free- molecular form
i1fm_at = kfmat * sqdgat5 * bm2ii(n2n) * ( &
esat25 &
+ two * esat09 * esat04 &
+ esat01 * esat16 &
+ esat64 * esat09 &
+ two * esat36 * esat01 &
+ esat16 * esat01 )
i1_at = ( i1nc_at * i1fm_at ) / ( i1nc_at + i1fm_at )
coagatat2 = constii * i1_at
qs11 = coagatat2 * bm2iitt(n2n)
! *** accumulation mode
! *** near-continuum
i1nc_ac = knc * dgac2 * ( &
two * esac16 &
+ esac04 * esac04 &
+ esac36 * esac04 &
+ a * kngac * ( &
two * esac04 &
+ esac16 * esac04 &
+ esac36 * esac16 ) )
! *** free- molecular form
i1fm_ac = kfmac * sqdgac5 * bm2ii(n2a) * ( &
esac25 &
+ two * esac09 * esac04 &
+ esac01 * esac16 &
+ esac64 * esac09 &
+ two * esac36 * esac01 &
+ esac16 * esac01 )
i1_ac = ( i1nc_ac * i1fm_ac ) / ( i1nc_ac + i1fm_ac )
coagacac2 = constii * i1_ac
qs22 = coagacac2 * bm2iitt(n2a)
return
end subroutine getcoags
!----------------------------------------------------------------------
!----------------------------------------------------------------------
end module modal_aero_coag