Line data    Source code

       1             : ! modal_aero_newnuc.F90
       2             : 
       3             : 
       4             : !----------------------------------------------------------------------
       5             : !BOP
       6             : !
       7             : ! !MODULE: modal_aero_newnuc --- modal aerosol new-particle nucleation
       8             : !
       9             : ! !INTERFACE:
      10             : module modal_aero_newnuc
      11             : 
      12             : ! !USES:
      13             :    use shr_kind_mod,  only:  r8 => shr_kind_r8
      14             :    use shr_kind_mod,  only:  r4 => shr_kind_r4
      15             :    use mo_constants,  only:  pi
      16             :    use chem_mods,     only:  gas_pcnst
      17             : 
      18             :   implicit none
      19             :   private
      20             :   save
      21             : 
      22             : ! !PUBLIC MEMBER FUNCTIONS:
      23             :   public modal_aero_newnuc_sub, modal_aero_newnuc_init
      24             : 
      25             : ! !PUBLIC DATA MEMBERS:
      26             :   integer, parameter  :: pcnstxx = gas_pcnst
      27             :   integer  :: l_h2so4_sv, l_nh3_sv, lnumait_sv, lnh4ait_sv, lso4ait_sv
      28             : 
      29             : ! min h2so4 vapor for nuc calcs = 4.0e-16 mol/mol-air ~= 1.0e4 molecules/cm3, 
      30             :   real(r8), parameter :: qh2so4_cutoff = 4.0e-16_r8
      31             : 
      32             :   real(r8) :: dens_so4a_host
      33             :   real(r8) :: mw_nh4a_host, mw_so4a_host
      34             : 
      35             : ! !DESCRIPTION: This module implements ...
      36             : !
      37             : ! !REVISION HISTORY:
      38             : !
      39             : !   R.Easter 2007.09.14:  Adapted from MIRAGE2 code
      40             : !
      41             : !EOP
      42             : !----------------------------------------------------------------------
      43             : !BOC
      44             : 
      45             : ! list private module data here
      46             : 
      47             : !EOC
      48             : !----------------------------------------------------------------------
      49             : 
      50             : 
      51             :   contains
      52             : 
      53             : !----------------------------------------------------------------------
      54             : !----------------------------------------------------------------------
      55             : !BOP
      56             : ! !ROUTINE:  modal_aero_newnuc_sub --- ...
      57             : !
      58             : ! !INTERFACE:
      59     1489176 :    subroutine modal_aero_newnuc_sub(                             &
      60             :                         lchnk,    ncol,     nstep,               &
      61             :                         loffset,  deltat,                        &
      62             :                         t,        pmid,     pdel,                &
      63             :                         zm,       pblh,                          &
      64     1489176 :                         qv,       cld,                           &
      65     1489176 :                         q,                                       &
      66     1489176 :                         del_h2so4_gasprod,  del_h2so4_aeruptk    )
      67             : 
      68             : 
      69             : ! !USES:
      70             :    use modal_aero_data
      71             :    use cam_abortutils,    only: endrun
      72             :    use cam_history,       only: outfld, fieldname_len
      73             :    use chem_mods,         only: adv_mass
      74             :    use constituents,      only: pcnst, cnst_name
      75             :    use physconst,         only: gravit, mwdry, r_universal
      76             :    use ppgrid,            only: pcols, pver
      77             :    use spmd_utils,        only: iam, masterproc
      78             :    use wv_saturation,     only: qsat
      79             :    use ref_pres,          only: top_lev=>clim_modal_aero_top_lev
      80             : 
      81             :    implicit none
      82             : 
      83             : ! !PARAMETERS:
      84             :    integer, intent(in)  :: lchnk            ! chunk identifier
      85             :    integer, intent(in)  :: ncol             ! number of columns in chunk
      86             :    integer, intent(in)  :: nstep            ! model step
      87             :    integer, intent(in)  :: loffset          ! offset applied to modal aero "pointers"
      88             :    real(r8), intent(in) :: deltat           ! model timestep (s)
      89             : 
      90             :    real(r8), intent(in) :: t(pcols,pver)    ! temperature (K)
      91             :    real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
      92             :    real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
      93             :    real(r8), intent(in) :: zm(pcols,pver)   ! midpoint height above surface (m)
      94             :    real(r8), intent(in) :: pblh(pcols)      ! pbl height (m)
      95             :    real(r8), intent(in) :: qv(pcols,pver)   ! specific humidity (kg/kg)
      96             :    real(r8), intent(in) :: cld(ncol,pver)   ! stratiform cloud fraction
      97             :                                             ! *** NOTE ncol dimension
      98             :    real(r8), intent(inout) :: q(ncol,pver,pcnstxx) 
      99             :                                             ! tracer mixing ratio (TMR) array
     100             :                                             ! *** MUST BE mol/mol-air or #/mol-air
     101             :                                             ! *** NOTE ncol & pcnstxx dimensions
     102             :    real(r8), intent(in) :: del_h2so4_gasprod(ncol,pver) 
     103             :                                             ! h2so4 gas-phase production
     104             :                                             ! change over deltat (mol/mol)
     105             :    real(r8), intent(in) :: del_h2so4_aeruptk(ncol,pver) 
     106             :                                             ! h2so4 gas-phase loss to
     107             :                                             ! aerosol over deltat (mol/mol)
     108             : 
     109             : ! !DESCRIPTION: 
     110             : !   computes changes due to aerosol nucleation (new particle formation)
     111             : !       treats both nucleation and subsequent growth of new particles
     112             : !           to aitken mode size
     113             : !   uses the following parameterizations
     114             : !       vehkamaki et al. (2002) parameterization for binary
     115             : !           homogeneous nucleation (h2so4-h2o) plus
     116             : !       kerminen and kulmala (2002) parameterization for
     117             : !           new particle loss during growth to aitken size
     118             : !
     119             : ! !REVISION HISTORY:
     120             : !   R.Easter 2007.09.14:  Adapted from MIRAGE2 code and CMAQ V4.6 code
     121             : !
     122             : !EOP
     123             : !----------------------------------------------------------------------
     124             : !BOC
     125             : 
     126             : !   local variables
     127             :         integer :: i, itmp, k, l, lmz, lun, m, mait
     128             :         integer :: lnumait, lso4ait, lnh4ait
     129             :         integer :: l_h2so4, l_nh3
     130             :         integer :: ldiagveh02
     131             :         integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
     132             :         integer, parameter :: newnuc_method_flagaa = 11
     133             : !       integer, parameter :: newnuc_method_flagaa = 12
     134             :         !  1=merikanto et al (2007) ternary   2=vehkamaki et al (2002) binary
     135             :         ! 11=merikanto ternary + first-order boundary layer
     136             :         ! 12=merikanto ternary + second-order boundary layer
     137             : 
     138             :         real(r8) :: adjust_factor
     139             :         real(r8) :: aircon
     140             :         real(r8) :: cldx 
     141             :         real(r8) :: dens_nh4so4a
     142             :         real(r8) :: dmdt_ait, dmdt_aitsv1, dmdt_aitsv2, dmdt_aitsv3
     143             :         real(r8) :: dndt_ait, dndt_aitsv1, dndt_aitsv2, dndt_aitsv3
     144             :         real(r8) :: dndt(pcols,pver) ! nucleation rate (#/m3/s)
     145             :         real(r8) :: dnh4dt_ait, dso4dt_ait
     146             :         real(r8) :: dpnuc
     147             :         real(r8) :: dplom_mode(1), dphim_mode(1)
     148             :         real(r8) :: ev_sat(pcols,pver)
     149             :         real(r8) :: mass1p
     150             :         real(r8) :: mass1p_aithi, mass1p_aitlo 
     151             :         real(r8) :: pdel_fac
     152             :         real(r8) :: qh2so4_cur, qh2so4_avg, qh2so4_del
     153             :         real(r8) :: qnh3_cur, qnh3_del, qnh4a_del
     154             :         real(r8) :: qnuma_del
     155             :         real(r8) :: qso4a_del
     156             :         real(r8) :: qv_sat(pcols,pver)
     157             :         real(r8) :: qvswtr
     158             :         real(r8) :: relhum, relhumav, relhumnn
     159             :         real(r8) :: tmpa, tmpb, tmpc
     160             :         real(r8) :: tmp_q1, tmp_q2, tmp_q3
     161             :         real(r8) :: tmp_frso4, tmp_uptkrate
     162             : 
     163             :         integer, parameter :: nsrflx = 1     ! last dimension of qsrflx
     164             :         real(r8) :: qsrflx(pcols,pcnst,nsrflx)
     165             :                               ! process-specific column tracer tendencies
     166             :                               ! 1 = nucleation (for aerocom)
     167     2978352 :         real(r8) :: dqdt(ncol,pver,pcnstxx)  ! TMR tendency array -- NOTE dims
     168             :         logical  :: dotend(pcnst)            ! flag for doing tendency
     169             :         logical  :: do_nh3                   ! flag for doing nh3/nh4
     170             : 
     171             : 
     172             :         character(len=1) :: tmpch1, tmpch2, tmpch3
     173             :         character(len=fieldname_len+3) :: fieldname
     174             : 
     175             : 
     176             : ! begin
     177     1489176 :         lun = 6
     178             : 
     179             : !--------------------------------------------------------------------------------
     180             : !!$   if (ldiag1 > 0) then
     181             : !!$   do i = 1, ncol
     182             : !!$   if (lonndx(i) /= 37) cycle
     183             : !!$   if (latndx(i) /= 23) cycle
     184             : !!$   if (nstep > 3)       cycle
     185             : !!$   write( lun, '(/a,i7,3i5,f10.2)' )   &
     186             : !!$         '*** modal_aero_newnuc_sub -- nstep, iam, lat, lon =',   &
     187             : !!$         nstep, iam, latndx(i), lonndx(i)
     188             : !!$   end do
     189             : !!$   if (nstep > 3) call endrun( '*** modal_aero_newnuc_sub -- testing halt after step 3' )
     190             : !!$!  if (ncol /= -999888777) return
     191             : !!$   end if
     192             : !--------------------------------------------------------------------------------
     193             : 
     194             : !-----------------------------------------------------------------------
     195     1489176 :         l_h2so4 = l_h2so4_sv - loffset
     196     1489176 :         l_nh3   = l_nh3_sv   - loffset
     197     1489176 :         lnumait = lnumait_sv - loffset
     198     1489176 :         lnh4ait = lnh4ait_sv - loffset
     199     1489176 :         lso4ait = lso4ait_sv - loffset
     200             : 
     201             : !   skip if no aitken mode OR if no h2so4 species
     202     1489176 :         if ((l_h2so4 <= 0) .or. (lso4ait <= 0) .or. (lnumait <= 0)) return
     203             : 
     204     1489176 :         dotend(:) = .false.
     205 71735685840 :         dqdt(1:ncol,:,:) = 0.0_r8
     206  1022475168 :         qsrflx(1:ncol,:,:) = 0.0_r8
     207  2314006344 :         dndt(1:ncol,:) = 0.0_r8
     208             : 
     209             : !   set dotend
     210     1489176 :         mait = modeptr_aitken
     211     1489176 :         dotend(lnumait) = .true.
     212     1489176 :         dotend(lso4ait) = .true.
     213     1489176 :         dotend(l_h2so4) = .true.
     214             : 
     215     1489176 :         lnh4ait = lptr_nh4_a_amode(mait) - loffset
     216             :         if ((l_nh3   > 0) .and. (l_nh3   <= pcnst) .and. &
     217     1489176 :             (lnh4ait > 0) .and. (lnh4ait <= pcnst)) then
     218           0 :             do_nh3 = .true.
     219           0 :             dotend(lnh4ait) = .true.
     220           0 :             dotend(l_nh3) = .true.
     221             :         else
     222             :             do_nh3 = .false.
     223             :         end if
     224             : 
     225             : 
     226             : !   dry-diameter limits for "grown" new particles
     227           0 :         dplom_mode(1) = exp( 0.67_r8*log(dgnumlo_amode(mait))   &
     228     1489176 :                            + 0.33_r8*log(dgnum_amode(mait)) )
     229     1489176 :         dphim_mode(1) = dgnumhi_amode(mait)
     230             : 
     231             : !   mass1p_... = mass (kg) of so4 & nh4 in a single particle of diameter ...
     232             : !                (assuming same dry density for so4 & nh4)
     233             : !       mass1p_aitlo - dp = dplom_mode(1)
     234             : !       mass1p_aithi - dp = dphim_mode(1)
     235     1489176 :         tmpa = dens_so4a_host*pi/6.0_r8
     236     1489176 :         mass1p_aitlo = tmpa*(dplom_mode(1)**3)
     237     1489176 :         mass1p_aithi = tmpa*(dphim_mode(1)**3)
     238             : 
     239             : !   compute qv_sat = saturation specific humidity
     240   139982544 :         do k = 1, pver
     241   139982544 :            call qsat(t(1:ncol,k), pmid(1:ncol,k), ev_sat(1:ncol,k), qv_sat(1:ncol,k), ncol)
     242             :         end do
     243             : !
     244             : !   loop over levels and columns to calc the renaming
     245             : !
     246   139982544 : main_k: do k = top_lev, pver
     247  2314006344 : main_i: do i = 1, ncol
     248             : 
     249             : !   skip if completely cloudy, 
     250             : !   because all h2so4 vapor should be cloud-borne
     251  2174023800 :         if (cld(i,k) >= 0.99_r8) cycle main_i
     252             : 
     253             : !   qh2so4_cur = current qh2so4, after aeruptk
     254  2045827763 :         qh2so4_cur = q(i,k,l_h2so4)
     255             : !   skip if h2so4 vapor < qh2so4_cutoff
     256  2045827763 :         if (qh2so4_cur <= qh2so4_cutoff) cycle main_i
     257             : 
     258  1033492327 :         tmpa = max( 0.0_r8, del_h2so4_gasprod(i,k) )
     259  1033492327 :         tmp_q3 = qh2so4_cur
     260             : !   tmp_q2 = qh2so4 before aeruptk
     261             : !   (note tmp_q3, tmp_q2 both >= 0.0)
     262  1033492327 :         tmp_q2 = tmp_q3 + max( 0.0_r8, -del_h2so4_aeruptk(i,k) )
     263             : 
     264             : !   *** temporary -- in order to get more nucleation
     265             : !       qh2so4_cur = qh2so4_cur*1.0e1
     266             : !       tmp_q3 = tmp_q3*1.0e1
     267             : !       tmp_q2 = tmp_q2*1.0e1
     268             : !       tmpa   = tmpa  *1.0e1
     269             : 
     270             : !   tmpb = log( tmp_q2/tmp_q3 ) BUT with some checks added
     271             : !       tmp_uptkrate = tmpb/deltat
     272  1033492327 :         if (tmp_q2 <= tmp_q3) then
     273             :            tmpb = 0.0_r8
     274             :         else
     275  1033492327 :            tmpc = tmp_q2 * exp( -20.0_r8 )
     276  1033492327 :            if (tmp_q3 <= tmpc) then
     277             :               tmp_q3 = tmpc
     278             :               tmpb = 20.0_r8
     279             :            else
     280  1033492327 :               tmpb = log( tmp_q2/tmp_q3 )
     281             :            end if
     282             :         end if
     283             : !   d[ln(qh2so4)]/dt (1/s) from uptake (condensation) to aerosol
     284  1033492327 :         tmp_uptkrate = tmpb/deltat
     285             : 
     286             : !   qh2so4_avg = estimated average qh2so4
     287             : !   when production & loss are done simultaneously
     288  1033492327 :         if (tmpb <= 0.1_r8) then
     289   120761731 :            qh2so4_avg = tmp_q3*(1.0_r8 + 0.5_r8*tmpb) - 0.5_r8*tmpa
     290             :         else
     291   912730596 :            tmpc = tmpa/tmpb
     292   912730596 :            qh2so4_avg = (tmp_q3 - tmpc)*((exp(tmpb)-1.0_r8)/tmpb) + tmpc
     293             :         end if
     294  1033492327 :         if (qh2so4_avg <= qh2so4_cutoff) cycle main_i
     295             : 
     296             : 
     297   834837703 :         if ( do_nh3 ) then
     298           0 :             qnh3_cur = max( 0.0_r8, q(i,k,l_nh3) )
     299             :         else
     300   834837703 :             qnh3_cur = 0.0_r8
     301             :         end if
     302             : 
     303             : 
     304             : !   relhumav = grid average RH
     305   834837703 :         qvswtr = qv_sat(i,k)
     306   834837703 :         qvswtr = max( qvswtr, 1.0e-20_r8 )
     307   834837703 :         relhumav = qv(i,k) / qvswtr
     308   834837703 :         relhumav = max( 0.0_r8, min( 1.0_r8, relhumav ) )
     309             : !   relhum = non-cloudy area RH
     310   834837703 :         cldx = max( 0.0_r8, cld(i,k) )
     311   834837703 :         relhum = (relhumav - cldx) / (1.0_r8 - cldx)
     312   834837703 :         relhum = max( 0.0_r8, min( 1.0_r8, relhum ) )
     313             : !   limit RH to between 0.1% and 99%
     314   834837703 :         relhumnn = relhum
     315   834837703 :         relhumnn = max( 0.01_r8, min( 0.99_r8, relhumnn ) )
     316             : 
     317             : !   aircon = air concentration (mol-air/m3)
     318   834837703 :         aircon = 1.0e3_r8*pmid(i,k)/(r_universal*t(i,k))
     319             : 
     320             : 
     321             : !   call ... routine to get nucleation rates
     322   834837703 :         ldiagveh02 = -1
     323             : !!$     if (ldiag2 > 0) then
     324             : !!$     if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
     325             : !!$     if ((k >= 24) .or. (mod(k,4) == 0)) then
     326             : !!$         ldiagveh02 = +1
     327             : !!$            write(lun,'(/a,i8,3i4,f8.2,1p,4e10.2)')   &
     328             : !!$             'veh02 call - nstep,lat,lon,k; tk,rh,p,cair',   &
     329             : !!$             nstep, latndx(i), lonndx(i), k,   &
     330             : !!$             t(i,k), relhumnn, pmid(k,k), aircon
     331             : !!$     end if
     332             : !!$     end if
     333             : !!$     end if
     334             :         call mer07_veh02_nuc_mosaic_1box(   &
     335             :            newnuc_method_flagaa,   &
     336             :            deltat, t(i,k), relhumnn, pmid(i,k),   &
     337             :            zm(i,k), pblh(i),   &
     338             :            qh2so4_cur, qh2so4_avg, qnh3_cur, tmp_uptkrate,   &
     339             :            mw_so4a_host,   &
     340             :            1, 1, dplom_mode, dphim_mode,   &
     341             :            itmp, qnuma_del, qso4a_del, qnh4a_del,   &
     342   834837703 :            qh2so4_del, qnh3_del, dens_nh4so4a, ldiagveh02 )
     343             : !          qh2so4_del, qnh3_del, dens_nh4so4a )
     344             : !----------------------------------------------------------------------
     345             : !       subr mer07_veh02_nuc_mosaic_1box(   &
     346             : !          newnuc_method_flagaa,   &
     347             : !          dtnuc, temp_in, rh_in, press_in,   &
     348             : !          qh2so4_cur, qh2so4_avg, qnh3_cur, h2so4_uptkrate,   &
     349             : !          nsize, maxd_asize, dplom_sect, dphim_sect,   &
     350             : !          isize_nuc, qnuma_del, qso4a_del, qnh4a_del,   &
     351             : !          qh2so4_del, qnh3_del, dens_nh4so4a )
     352             : !
     353             : !! subr arguments (in)
     354             : !        real(r8), intent(in) :: dtnuc             ! nucleation time step (s)
     355             : !        real(r8), intent(in) :: temp_in           ! temperature, in k
     356             : !        real(r8), intent(in) :: rh_in             ! relative humidity, as fraction
     357             : !        real(r8), intent(in) :: press_in          ! air pressure (pa)
     358             : !
     359             : !        real(r8), intent(in) :: qh2so4_cur, qh2so4_avg
     360             : !                                                  ! gas h2so4 mixing ratios (mol/mol-air)
     361             : !        real(r8), intent(in) :: qnh3_cur          ! gas nh3 mixing ratios (mol/mol-air)
     362             : !             ! qxxx_cur = current value (after gas chem and condensation)
     363             : !             ! qxxx_avg = estimated average value (for simultaneous source/sink calcs)
     364             : !        real(r8), intent(in) :: h2so4_uptkrate    ! h2so4 uptake rate to aerosol (1/s)
     365             : 
     366             : !
     367             : !        integer, intent(in) :: nsize                    ! number of aerosol size bins
     368             : !        integer, intent(in) :: maxd_asize               ! dimension for dplom_sect, ...
     369             : !        real(r8), intent(in) :: dplom_sect(maxd_asize)  ! dry diameter at lower bnd of bin (m)
     370             : !        real(r8), intent(in) :: dphim_sect(maxd_asize)  ! dry diameter at upper bnd of bin (m)
     371             : !
     372             : !! subr arguments (out)
     373             : !        integer, intent(out) :: isize_nuc         ! size bin into which new particles go
     374             : !        real(r8), intent(out) :: qnuma_del        ! change to aerosol number mixing ratio (#/mol-air)
     375             : !        real(r8), intent(out) :: qso4a_del        ! change to aerosol so4 mixing ratio (mol/mol-air)
     376             : !        real(r8), intent(out) :: qnh4a_del        ! change to aerosol nh4 mixing ratio (mol/mol-air)
     377             : !        real(r8), intent(out) :: qh2so4_del       ! change to gas h2so4 mixing ratio (mol/mol-air)
     378             : !        real(r8), intent(out) :: qnh3_del         ! change to gas nh3 mixing ratio (mol/mol-air)
     379             : !                                                  ! aerosol changes are > 0; gas changes are < 0
     380             : !        real(r8), intent(out) :: dens_nh4so4a     ! dry-density of the new nh4-so4 aerosol mass (kg/m3)
     381             : !----------------------------------------------------------------------
     382             : 
     383             : 
     384             : !   convert qnuma_del from (#/mol-air) to (#/kmol-air)
     385   834837703 :         qnuma_del = qnuma_del*1.0e3_r8
     386             : !   number nuc rate (#/kmol-air/s) from number nuc amt
     387   834837703 :         dndt_ait = qnuma_del/deltat
     388             : !   fraction of mass nuc going to so4
     389   834837703 :         tmpa = qso4a_del*mw_so4a_host
     390   834837703 :         tmpb = tmpa + qnh4a_del*mw_nh4a_host
     391   834837703 :         tmp_frso4 = max( tmpa, 1.0e-35_r8 )/max( tmpb, 1.0e-35_r8 )
     392             : !   mass nuc rate (kg/kmol-air/s or g/mol...) hhfrom mass nuc amts
     393   834837703 :         dmdt_ait = max( 0.0_r8, (tmpb/deltat) ) 
     394             : 
     395   834837703 :         dndt_aitsv1 = dndt_ait
     396   834837703 :         dmdt_aitsv1 = dmdt_ait
     397   834837703 :         dndt_aitsv2 = 0.0_r8
     398   834837703 :         dmdt_aitsv2 = 0.0_r8
     399   834837703 :         dndt_aitsv3 = 0.0_r8
     400   834837703 :         dmdt_aitsv3 = 0.0_r8
     401   834837703 :         tmpch1 = ' '
     402   834837703 :         tmpch2 = ' '
     403             : 
     404   834837703 :         if (dndt_ait < 1.0e2_r8) then
     405             : !   ignore newnuc if number rate < 100 #/kmol-air/s ~= 0.3 #/mg-air/d
     406   619405355 :             dndt_ait = 0.0_r8
     407   619405355 :             dmdt_ait = 0.0_r8
     408   619405355 :             tmpch1 = 'A'
     409             : 
     410             :         else
     411   215432348 :             dndt_aitsv2 = dndt_ait
     412   215432348 :             dmdt_aitsv2 = dmdt_ait
     413   215432348 :             tmpch1 = 'B'
     414             : 
     415             : !   mirage2 code checked for complete h2so4 depletion here,
     416             : !   but this is now done in mer07_veh02_nuc_mosaic_1box
     417   215432348 :             mass1p = dmdt_ait/dndt_ait
     418   215432348 :             dndt_aitsv3 = dndt_ait
     419   215432348 :             dmdt_aitsv3 = dmdt_ait
     420             : 
     421             : !   apply particle size constraints
     422   215432348 :             if (mass1p < mass1p_aitlo) then
     423             : !   reduce dndt to increase new particle size
     424           0 :                 dndt_ait = dmdt_ait/mass1p_aitlo
     425           0 :                 tmpch1 = 'C'
     426   215432348 :             else if (mass1p > mass1p_aithi) then
     427             : !   reduce dmdt to decrease new particle size
     428           0 :                 dmdt_ait = dndt_ait*mass1p_aithi
     429           0 :                 tmpch1 = 'E'
     430             :             end if
     431             :         end if
     432             : 
     433             : ! *** apply adjustment factor to avoid unrealistically high
     434             : !     aitken number concentrations in mid and upper troposphere
     435             : !       adjust_factor = 0.5
     436             : !       dndt_ait = dndt_ait * adjust_factor
     437             : !       dmdt_ait = dmdt_ait * adjust_factor
     438             : 
     439             : !   set tendencies
     440   834837703 :         pdel_fac = pdel(i,k)/gravit
     441             : 
     442             : !   dso4dt_ait, dnh4dt_ait are (kmol/kmol-air/s)
     443   834837703 :         dso4dt_ait = dmdt_ait*tmp_frso4/mw_so4a_host
     444   834837703 :         dnh4dt_ait = dmdt_ait*(1.0_r8 - tmp_frso4)/mw_nh4a_host
     445             : 
     446   834837703 :         dqdt(i,k,l_h2so4) = -dso4dt_ait*(1.0_r8-cldx)
     447   834837703 :         qsrflx(i,l_h2so4,1) = qsrflx(i,l_h2so4,1) + dqdt(i,k,l_h2so4)*pdel_fac
     448   834837703 :         q(i,k,l_h2so4) = q(i,k,l_h2so4) + dqdt(i,k,l_h2so4)*deltat
     449             : 
     450   834837703 :         dqdt(i,k,lso4ait) = dso4dt_ait*(1.0_r8-cldx)
     451   834837703 :         qsrflx(i,lso4ait,1) = qsrflx(i,lso4ait,1) + dqdt(i,k,lso4ait)*pdel_fac
     452   834837703 :         q(i,k,lso4ait) = q(i,k,lso4ait) + dqdt(i,k,lso4ait)*deltat
     453   834837703 :         if (lnumait > 0) then
     454   834837703 :             dqdt(i,k,lnumait) = dndt_ait*(1.0_r8-cldx)
     455             : !   dndt is (#/m3/s), dqdt(:,:,lnumait) is (#/kmol-air/s), aircon is (mol-air/m3)
     456   834837703 :             dndt(i,k) = dqdt(i,k,lnumait)*aircon*1.0e-3_r8
     457             :             qsrflx(i,lnumait,1) = qsrflx(i,lnumait,1)   &
     458   834837703 :                                 + dqdt(i,k,lnumait)*pdel_fac
     459   834837703 :             q(i,k,lnumait) = q(i,k,lnumait) + dqdt(i,k,lnumait)*deltat
     460             :         end if
     461             : 
     462   973331071 :         if (( do_nh3 ) .and. (dnh4dt_ait > 0.0_r8)) then
     463           0 :             dqdt(i,k,l_nh3) = -dnh4dt_ait*(1.0_r8-cldx)
     464           0 :             qsrflx(i,l_nh3,1) = qsrflx(i,l_nh3,1) + dqdt(i,k,l_nh3)*pdel_fac
     465           0 :             q(i,k,l_nh3) = q(i,k,l_nh3) + dqdt(i,k,l_nh3)*deltat
     466             : 
     467           0 :             dqdt(i,k,lnh4ait) = dnh4dt_ait*(1.0_r8-cldx)
     468           0 :             qsrflx(i,lnh4ait,1) = qsrflx(i,lnh4ait,1) + dqdt(i,k,lnh4ait)*pdel_fac
     469           0 :             q(i,k,lnh4ait) = q(i,k,lnh4ait) + dqdt(i,k,lnh4ait)*deltat
     470             :         end if
     471             : 
     472             : !!   temporary diagnostic
     473             : !        if (ldiag3 > 0) then
     474             : !        if ((dndt_ait /= 0.0_r8) .or. (dmdt_ait /= 0.0_r8)) then
     475             : !           write(lun,'(3a,1x,i7,3i5,1p,5e12.4)')   &
     476             : !              'newnucxx', tmpch1, tmpch2, nstep, lchnk, i, k,   &
     477             : !              dndt_ait, dmdt_ait, cldx
     478             : !!          call endrun( 'modal_aero_newnuc_sub' )
     479             : !        end if
     480             : !        end if
     481             : 
     482             : 
     483             : !   diagnostic output start ----------------------------------------
     484             : !!$     if (ldiag4 > 0) then
     485             : !!$     if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
     486             : !!$     if ((k >= 24) .or. (mod(k,4) == 0)) then
     487             : !!$        write(lun,97010) nstep, latndx(i), lonndx(i), k, t(i,k), aircon
     488             : !!$        write(lun,97020) 'pmid, pdel                   ',   &
     489             : !!$                pmid(i,k), pdel(i,k)
     490             : !!$        write(lun,97030) 'qv,qvsw, cld, rh_av, rh_clr  ',   &
     491             : !!$                qv(i,k), qvswtr, cldx, relhumav, relhum
     492             : !!$        write(lun,97020) 'h2so4_cur, _pre, _av, nh3_cur',   &
     493             : !!$             qh2so4_cur, tmp_q2, qh2so4_avg, qnh3_cur
     494             : !!$        write(lun,97020) 'del_h2so4_gasprod, _aeruptk  ',   &
     495             : !!$             del_h2so4_gasprod(i,k), del_h2so4_aeruptk(i,k),   &
     496             : !!$             tmp_uptkrate*3600.0_r8
     497             : !!$        write(lun,97020) ' '
     498             : !!$        write(lun,97050) 'tmpch1, tmpch2               ', tmpch1, tmpch2
     499             : !!$        write(lun,97020) 'dndt_, dmdt_aitsv1           ',   &
     500             : !!$                              dndt_aitsv1, dmdt_aitsv1
     501             : !!$        write(lun,97020) 'dndt_, dmdt_aitsv2           ',   &
     502             : !!$                              dndt_aitsv2, dmdt_aitsv2
     503             : !!$        write(lun,97020) 'dndt_, dmdt_aitsv3           ',   &
     504             : !!$                              dndt_aitsv3, dmdt_aitsv3
     505             : !!$        write(lun,97020) 'dndt_, dmdt_ait              ',   &
     506             : !!$                              dndt_ait, dmdt_ait
     507             : !!$        write(lun,97020) 'dso4dt_, dnh4dt_ait          ',   &
     508             : !!$                              dso4dt_ait, dnh4dt_ait
     509             : !!$        write(lun,97020) 'qso4a_del, qh2so4_del        ',   &
     510             : !!$                              qso4a_del, qh2so4_del
     511             : !!$        write(lun,97020) 'qnh4a_del, qnh3_del          ',   &
     512             : !!$                              qnh4a_del, qnh3_del
     513             : !!$        write(lun,97020) 'dqdt(h2so4), (nh3)           ',   &
     514             : !!$              dqdt(i,k,l_h2so4), dqdt(i,k,l_nh3) 
     515             : !!$        write(lun,97020) 'dqdt(so4a), (nh4a), (numa)   ',   &
     516             : !!$              dqdt(i,k,lso4ait), dqdt(i,k,lnh4ait), dqdt(i,k,lnumait)
     517             : !!$ 
     518             : !!$     dpnuc = 0.0_r8
     519             : !!$     if (dndt_aitsv1 > 1.0e-5_r8) dpnuc = (6.0_r8*dmdt_aitsv1/   &
     520             : !!$                     (pi*dens_so4a_host*dndt_aitsv1))**0.3333333_r8
     521             : !!$        if (dpnuc > 0.0_r8) then
     522             : !!$        write(lun,97020) 'dpnuc,      dp_aitlo, _aithi ',   &
     523             : !!$                      dpnuc, dplom_mode(1), dphim_mode(1)
     524             : !!$        write(lun,97020) 'mass1p, mass1p_aitlo, _aithi ',   &
     525             : !!$                      mass1p, mass1p_aitlo, mass1p_aithi
     526             : !!$        end if
     527             : !!$ 
     528             : !!$ 97010  format( / 'NEWNUC nstep,lat,lon,k,tk,cair', i8, 3i4, f8.2, 1pe12.4 )
     529             : !!$ 97020  format( a, 1p, 6e12.4 )
     530             : !!$ 97030  format( a, 1p, 2e12.4, 0p, 5f10.6 )
     531             : !!$ 97040  format( 29x, 1p, 6e12.4 )
     532             : !!$ 97050  format( a, 2(3x,a) )
     533             : !!$        end if
     534             : !!$        end if
     535             : !!$        end if
     536             : !   diagnostic output end   ------------------------------------------
     537             : 
     538             : 
     539             :         end do main_i
     540             :         end do main_k
     541             : 
     542             : 
     543             : !   do history file column-tendency fields
     544    46164456 :         do l = loffset+1, pcnst
     545    44675280 :             lmz = l - loffset
     546    44675280 :             if ( .not. dotend(lmz) ) cycle
     547             : 
     548    74597328 :             do i = 1, ncol
     549    74597328 :                 qsrflx(i,lmz,1) = qsrflx(i,lmz,1)*(adv_mass(lmz)/mwdry)
     550             :             end do
     551     4467528 :             fieldname = trim(cnst_name(l)) // '_sfnnuc1'
     552    46164456 :             call outfld( fieldname, qsrflx(:,lmz,1), pcols, lchnk )
     553             : 
     554             : !           if (( masterproc ) .and. (nstep < 1)) &
     555             : !               write(lun,'(2(a,2x),1p,e11.3)') &
     556             : !               'modal_aero_newnuc_sub outfld', fieldname, adv_mass(lmz)
     557             :         end do ! l = ...
     558             : 
     559             : 
     560             :         return
     561             : !EOC
     562     1489176 :         end subroutine modal_aero_newnuc_sub
     563             : 
     564             : 
     565             : 
     566             : !----------------------------------------------------------------------
     567             : !-----------------------------------------------------------------------
     568   834837703 :         subroutine mer07_veh02_nuc_mosaic_1box(   &
     569             :            newnuc_method_flagaa, dtnuc, temp_in, rh_in, press_in,   &
     570             :            zm_in, pblh_in,   &
     571             :            qh2so4_cur, qh2so4_avg, qnh3_cur, h2so4_uptkrate,   &
     572             :            mw_so4a_host,   &
     573   834837703 :            nsize, maxd_asize, dplom_sect, dphim_sect,   &
     574             :            isize_nuc, qnuma_del, qso4a_del, qnh4a_del,   &
     575             :            qh2so4_del, qnh3_del, dens_nh4so4a, ldiagaa )
     576             : !          qh2so4_del, qnh3_del, dens_nh4so4a )
     577     1489176 :           use mo_constants, only: rgas, &               ! Gas constant (J/K/kmol)
     578             :                                   avogad => avogadro    ! Avogadro's number (1/kmol)
     579             :           use physconst,    only: mw_so4a => mwso4, &   ! Molecular weight of sulfate
     580             :                                   mw_nh4a => mwnh4      ! Molecular weight of ammonium
     581             : !.......................................................................
     582             : !
     583             : ! calculates new particle production from homogeneous nucleation
     584             : !    over timestep dtnuc, using nucleation rates from either
     585             : !    merikanto et al. (2007) h2so4-nh3-h2o ternary parameterization
     586             : !    vehkamaki et al. (2002) h2so4-h2o binary parameterization
     587             : !
     588             : ! the new particles are "grown" to the lower-bound size of the host code's 
     589             : !    smallest size bin.  (this "growth" is somewhat ad hoc, and would not be
     590             : !    necessary if the host code's size bins extend down to ~1 nm.)
     591             : !
     592             : !    if the h2so4 and nh3 mass mixing ratios (mixrats) of the grown new 
     593             : !    particles exceed the current gas mixrats, the new particle production
     594             : !    is reduced so that the new particle mass mixrats match the gas mixrats.
     595             : !
     596             : !    the correction of kerminen and kulmala (2002) is applied to account
     597             : !    for loss of the new particles by coagulation as they are
     598             : !    growing to the "host code mininum size"
     599             : !
     600             : ! revision history
     601             : !    coded by rc easter, pnnl, xx-apr-2007
     602             : !
     603             : ! key routines called: subr ternary_nuc_napari
     604             : !
     605             : ! references:
     606             : !    merikanto, j., i. napari, h. vehkamaki, t. anttila,
     607             : !     and m. kulmala, 2007, new parameterization of
     608             : !     sulfuric acid-ammonia-water ternary nucleation
     609             : !     rates at tropospheric conditions,
     610             : !       j. geophys. res., 112, d15207, doi:10.1029/2006jd0027977
     611             : !
     612             : !    vehkamäki, h., m. kulmala, i. napari, k.e.j. lehtinen,
     613             : !       c. timmreck, m. noppel and a. laaksonen, 2002,
     614             : !       an improved parameterization for sulfuric acid-water nucleation
     615             : !       rates for tropospheric and stratospheric conditions,
     616             : !       j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
     617             : !
     618             : !    kerminen, v., and m. kulmala, 2002,
     619             : !       analytical formulae connecting the "real" and the "apparent"
     620             : !       nucleation rate and the nuclei number concentration
     621             : !       for atmospheric nucleation events
     622             : !
     623             : !.......................................................................
     624             :       implicit none
     625             : 
     626             : ! subr arguments (in)
     627             :         real(r8), intent(in) :: dtnuc             ! nucleation time step (s)
     628             :         real(r8), intent(in) :: temp_in           ! temperature, in k
     629             :         real(r8), intent(in) :: rh_in             ! relative humidity, as fraction
     630             :         real(r8), intent(in) :: press_in          ! air pressure (pa)
     631             :         real(r8), intent(in) :: zm_in             ! layer midpoint height (m)
     632             :         real(r8), intent(in) :: pblh_in           ! pbl height (m)
     633             : 
     634             :         real(r8), intent(in) :: qh2so4_cur, qh2so4_avg
     635             :                                                   ! gas h2so4 mixing ratios (mol/mol-air)
     636             :         real(r8), intent(in) :: qnh3_cur          ! gas nh3 mixing ratios (mol/mol-air)
     637             :              ! qxxx_cur = current value (after gas chem and condensation)
     638             :              ! qxxx_avg = estimated average value (for simultaneous source/sink calcs)
     639             :         real(r8), intent(in) :: h2so4_uptkrate    ! h2so4 uptake rate to aerosol (1/s)
     640             :         real(r8), intent(in) :: mw_so4a_host      ! mw of so4 aerosol in host code (g/mol)
     641             : 
     642             :         integer, intent(in) :: newnuc_method_flagaa     ! 1=merikanto et al (2007) ternary
     643             :                                                         ! 2=vehkamaki et al (2002) binary
     644             :         integer, intent(in) :: nsize                    ! number of aerosol size bins
     645             :         integer, intent(in) :: maxd_asize               ! dimension for dplom_sect, ...
     646             :         real(r8), intent(in) :: dplom_sect(maxd_asize)  ! dry diameter at lower bnd of bin (m)
     647             :         real(r8), intent(in) :: dphim_sect(maxd_asize)  ! dry diameter at upper bnd of bin (m)
     648             :         integer, intent(in) :: ldiagaa
     649             : 
     650             : ! subr arguments (out)
     651             :         integer, intent(out) :: isize_nuc         ! size bin into which new particles go
     652             :         real(r8), intent(out) :: qnuma_del        ! change to aerosol number mixing ratio (#/mol-air)
     653             :         real(r8), intent(out) :: qso4a_del        ! change to aerosol so4 mixing ratio (mol/mol-air)
     654             :         real(r8), intent(out) :: qnh4a_del        ! change to aerosol nh4 mixing ratio (mol/mol-air)
     655             :         real(r8), intent(out) :: qh2so4_del       ! change to gas h2so4 mixing ratio (mol/mol-air)
     656             :         real(r8), intent(out) :: qnh3_del         ! change to gas nh3 mixing ratio (mol/mol-air)
     657             :                                                   ! aerosol changes are > 0; gas changes are < 0
     658             :         real(r8), intent(out) :: dens_nh4so4a     ! dry-density of the new nh4-so4 aerosol mass (kg/m3)
     659             : 
     660             : ! subr arguments (out) passed via common block  
     661             : !    these are used to duplicate the outputs of yang zhang's original test driver
     662             : !    they are not really needed in wrf-chem
     663             :         real(r8) :: ratenuclt        ! j = ternary nucleation rate from napari param. (cm-3 s-1)
     664             :         real(r8) :: rateloge         ! ln (j)
     665             :         real(r8) :: cnum_h2so4       ! number of h2so4 molecules in the critical nucleus
     666             :         real(r8) :: cnum_nh3         ! number of nh3   molecules in the critical nucleus
     667             :         real(r8) :: cnum_tot         ! total number of molecules in the critical nucleus
     668             :         real(r8) :: radius_cluster   ! the radius of cluster (nm)
     669             : 
     670             : 
     671             : ! local variables
     672             :         integer :: i
     673             :         integer :: igrow
     674             :         integer, save :: icase = 0, icase_reldiffmax = 0
     675             : !       integer, parameter :: ldiagaa = -1
     676             :         integer :: lun
     677             :         integer :: newnuc_method_flagaa2
     678             : 
     679             :         real(r8), parameter :: onethird = 1.0_r8/3.0_r8
     680             : 
     681             :         real(r8), parameter :: accom_coef_h2so4 = 0.65_r8   ! accomodation coef for h2so4 conden
     682             : 
     683             : ! dry densities (kg/m3) molecular weights of aerosol 
     684             : ! ammsulf, ammbisulf, and sulfacid (from mosaic  dens_electrolyte values)
     685             : !       real(r8), parameter :: dens_ammsulf   = 1.769e3
     686             : !       real(r8), parameter :: dens_ammbisulf = 1.78e3
     687             : !       real(r8), parameter :: dens_sulfacid  = 1.841e3
     688             : ! use following to match cam3 modal_aero densities
     689             :         real(r8), parameter :: dens_ammsulf   = 1.770e3_r8
     690             :         real(r8), parameter :: dens_ammbisulf = 1.770e3_r8
     691             :         real(r8), parameter :: dens_sulfacid  = 1.770e3_r8
     692             : 
     693             : ! molecular weights (g/mol) of aerosol ammsulf, ammbisulf, and sulfacid
     694             : !    for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
     695             : !    because we don't keep track of aerosol hion mass
     696             :         real(r8), parameter :: mw_ammsulf   = 132.0_r8
     697             :         real(r8), parameter :: mw_ammbisulf = 114.0_r8
     698             :         real(r8), parameter :: mw_sulfacid  =  96.0_r8
     699             : 
     700             :         real(r8), save :: reldiffmax = 0.0_r8
     701             : 
     702             :         real(r8) cair                     ! dry-air molar density (mol/m3)
     703             :         real(r8) cs_prime_kk              ! kk2002 "cs_prime" parameter (1/m2)
     704             :         real(r8) cs_kk                    ! kk2002 "cs" parameter (1/s)
     705             :         real(r8) dens_part                ! "grown" single-particle dry density (kg/m3)
     706             :         real(r8) dfin_kk, dnuc_kk         ! kk2002 final/initial new particle wet diameter (nm)
     707             :         real(r8) dpdry_clus               ! critical cluster diameter (m)
     708             :         real(r8) dpdry_part               ! "grown" single-particle dry diameter (m)
     709             :         real(r8) tmpa, tmpb, tmpc, tmpe, tmpq
     710             :         real(r8) tmpa1, tmpb1
     711             :         real(r8) tmp_m1, tmp_m2, tmp_m3, tmp_n1, tmp_n2, tmp_n3
     712             :         real(r8) tmp_spd                  ! h2so4 vapor molecular speed (m/s)
     713             :         real(r8) factor_kk
     714             :         real(r8) fogas, foso4a, fonh4a, fonuma
     715             :         real(r8) freduce                  ! reduction factor applied to nucleation rate
     716             :                                           ! due to limited availability of h2so4 & nh3 gases
     717             :         real(r8) freducea, freduceb
     718             :         real(r8) gamma_kk                 ! kk2002 "gamma" parameter (nm2*m2/h)
     719             :         real(r8) gr_kk                    ! kk2002 "gr" parameter (nm/h)
     720             :         real(r8) kgaero_per_moleso4a      ! (kg dry aerosol)/(mol aerosol so4)
     721             :         real(r8) mass_part                ! "grown" single-particle dry mass (kg)
     722             :         real(r8) molenh4a_per_moleso4a    ! (mol aerosol nh4)/(mol aerosol so4)
     723             :         real(r8) nh3ppt, nh3ppt_bb        ! actual and bounded nh3 (ppt)
     724             :         real(r8) nu_kk                    ! kk2002 "nu" parameter (nm)
     725             :         real(r8) qmolnh4a_del_max         ! max production of aerosol nh4 over dtnuc (mol/mol-air)
     726             :         real(r8) qmolso4a_del_max         ! max production of aerosol so4 over dtnuc (mol/mol-air)
     727             :         real(r8) ratenuclt_bb             ! nucleation rate (#/m3/s)
     728             :         real(r8) ratenuclt_kk             ! nucleation rate after kk2002 adjustment (#/m3/s)
     729             :         real(r8) rh_bb                    ! bounded value of rh_in
     730             :         real(r8) so4vol_in                ! concentration of h2so4 for nucl. calc., molecules cm-3
     731             :         real(r8) so4vol_bb                ! bounded value of so4vol_in
     732             :         real(r8) temp_bb                  ! bounded value of temp_in
     733             :         real(r8) voldry_clus              ! critical-cluster dry volume (m3)
     734             :         real(r8) voldry_part              ! "grown" single-particle dry volume (m3)
     735             :         real(r8) wetvol_dryvol            ! grown particle (wet-volume)/(dry-volume)
     736             :         real(r8) wet_volfrac_so4a         ! grown particle (dry-volume-from-so4)/(wet-volume)
     737             : 
     738             : 
     739             : 
     740             : !
     741             : ! if h2so4 vapor < qh2so4_cutoff
     742             : ! exit with new particle formation = 0
     743             : !
     744   834837703 :         isize_nuc = 1
     745   834837703 :         qnuma_del = 0.0_r8
     746   834837703 :         qso4a_del = 0.0_r8
     747   834837703 :         qnh4a_del = 0.0_r8
     748   834837703 :         qh2so4_del = 0.0_r8
     749   834837703 :         qnh3_del = 0.0_r8
     750             : !       if (qh2so4_avg .le. qh2so4_cutoff) return   ! this no longer needed
     751             : !       if (qh2so4_cur .le. qh2so4_cutoff) return   ! this no longer needed
     752             : 
     753             :         if ((newnuc_method_flagaa /=  1) .and. &
     754             :             (newnuc_method_flagaa /=  2) .and. &
     755   834837703 :             (newnuc_method_flagaa /= 11) .and. &
     756             :             (newnuc_method_flagaa /= 12)) return
     757             : 
     758             : 
     759             : !
     760             : ! make call to parameterization routine
     761             : !
     762             : 
     763             : ! calc h2so4 in molecules/cm3 and nh3 in ppt
     764   834837703 :         cair = press_in/(temp_in*rgas)
     765   834837703 :         so4vol_in  = qh2so4_avg * cair * avogad * 1.0e-6_r8
     766   834837703 :         nh3ppt    = qnh3_cur * 1.0e12_r8
     767   834837703 :         ratenuclt = 1.0e-38_r8
     768   834837703 :         rateloge = log( ratenuclt )
     769             : 
     770   834837703 :         if ( (newnuc_method_flagaa /=  2) .and. &
     771             :              (nh3ppt >= 0.1_r8) ) then
     772             : ! make call to merikanto ternary parameterization routine
     773             : ! (when nh3ppt < 0.1, use binary param instead)
     774             : 
     775           0 :             if (so4vol_in >= 5.0e4_r8) then
     776           0 :                temp_bb = max( 235.0_r8, min( 295.0_r8, temp_in ) )
     777           0 :                rh_bb = max( 0.05_r8, min( 0.95_r8, rh_in ) )
     778           0 :                so4vol_bb = max( 5.0e4_r8, min( 1.0e9_r8, so4vol_in ) )
     779           0 :                nh3ppt_bb = max( 0.1_r8, min( 1.0e3_r8, nh3ppt ) )
     780             :                call ternary_nuc_merik2007(   &
     781             :                   temp_bb, rh_bb, so4vol_bb, nh3ppt_bb,   &
     782             :                   rateloge,   &
     783           0 :                   cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
     784             :             end if
     785           0 :             newnuc_method_flagaa2 = 1
     786             : 
     787             :         else
     788             : ! make call to vehkamaki binary parameterization routine
     789             : 
     790   834837703 :             if (so4vol_in >= 1.0e4_r8) then
     791   739780990 :                temp_bb = max( 230.15_r8, min( 305.15_r8, temp_in ) )
     792   739780990 :                rh_bb = max( 1.0e-4_r8, min( 1.0_r8, rh_in ) )
     793   739780990 :                so4vol_bb = max( 1.0e4_r8, min( 1.0e11_r8, so4vol_in ) )
     794             :                call binary_nuc_vehk2002(   &
     795             :                   temp_bb, rh_bb, so4vol_bb,   &
     796             :                   ratenuclt, rateloge,   &
     797   739780990 :                   cnum_h2so4, cnum_tot, radius_cluster )
     798             :             end if
     799   834837703 :             cnum_nh3 = 0.0_r8
     800   834837703 :             newnuc_method_flagaa2 = 2
     801             : 
     802             :         end if
     803             : 
     804             : 
     805             : ! do boundary layer nuc
     806   834837703 :         if ((newnuc_method_flagaa == 11) .or.   &
     807             :             (newnuc_method_flagaa == 12)) then
     808   834837703 :            if ( zm_in <= max(pblh_in,100.0_r8) ) then
     809    22327328 :               so4vol_bb = so4vol_in
     810             :               call pbl_nuc_wang2008( so4vol_bb,   &
     811             :                  newnuc_method_flagaa, newnuc_method_flagaa2,   &
     812             :                  ratenuclt, rateloge,   &
     813    22327328 :                  cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
     814             :            end if
     815             :         end if
     816             : 
     817             : 
     818             : ! if nucleation rate is less than 1e-6 #/m3/s ~= 0.1 #/cm3/day,
     819             : ! exit with new particle formation = 0
     820   834837703 :         if (rateloge  .le. -13.82_r8) return
     821             : !       if (ratenuclt .le. 1.0e-6) return
     822   278834791 :         ratenuclt = exp( rateloge )
     823   278834791 :         ratenuclt_bb = ratenuclt*1.0e6_r8
     824             : 
     825             : 
     826             : ! wet/dry volume ratio - use simple kohler approx for ammsulf/ammbisulf
     827   278834791 :         tmpa = max( 0.10_r8, min( 0.95_r8, rh_in ) )
     828   278834791 :         wetvol_dryvol = 1.0_r8 - 0.56_r8/log(tmpa)
     829             : 
     830             : 
     831             : ! determine size bin into which the new particles go
     832             : ! (probably it will always be bin #1, but ...)
     833             :         voldry_clus = ( max(cnum_h2so4,1.0_r8)*mw_so4a + cnum_nh3*mw_nh4a ) /   &
     834   278834791 :                       (1.0e3_r8*dens_sulfacid*avogad)
     835             : ! correction when host code sulfate is really ammonium bisulfate/sulfate
     836   278834791 :         voldry_clus = voldry_clus * (mw_so4a_host/mw_so4a)
     837   278834791 :         dpdry_clus = (voldry_clus*6.0_r8/pi)**onethird
     838             : 
     839   278834791 :         isize_nuc = 1
     840   278834791 :         dpdry_part = dplom_sect(1)
     841   278834791 :         if (dpdry_clus <= dplom_sect(1)) then
     842             :            igrow = 1   ! need to clusters to larger size
     843           0 :         else if (dpdry_clus >= dphim_sect(nsize)) then
     844           0 :            igrow = 0
     845           0 :            isize_nuc = nsize
     846             :            dpdry_part = dphim_sect(nsize)
     847             :         else
     848           0 :            igrow = 0
     849           0 :            do i = 1, nsize
     850           0 :               if (dpdry_clus < dphim_sect(i)) then
     851           0 :                  isize_nuc = i
     852           0 :                  dpdry_part = dpdry_clus
     853           0 :                  dpdry_part = min( dpdry_part, dphim_sect(i) )
     854           0 :                  dpdry_part = max( dpdry_part, dplom_sect(i) )
     855           0 :                  exit
     856             :               end if
     857             :            end do
     858             :         end if
     859   278834791 :         voldry_part = (pi/6.0_r8)*(dpdry_part**3)
     860             : 
     861             : 
     862             : !
     863             : ! determine composition and density of the "grown particles"
     864             : ! the grown particles are assumed to be liquid
     865             : !    (since critical clusters contain water)
     866             : !    so any (nh4/so4) molar ratio between 0 and 2 is allowed
     867             : ! assume that the grown particles will have 
     868             : !    (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
     869             : !
     870   278834791 :         if (igrow .le. 0) then
     871             : ! no "growing" so pure sulfuric acid
     872             :            tmp_n1 = 0.0_r8
     873             :            tmp_n2 = 0.0_r8
     874             :            tmp_n3 = 1.0_r8
     875   278834791 :         else if (qnh3_cur .ge. qh2so4_cur) then
     876             : ! combination of ammonium sulfate and ammonium bisulfate
     877             : ! tmp_n1 & tmp_n2 = mole fractions of the ammsulf & ammbisulf
     878           0 :            tmp_n1 = (qnh3_cur/qh2so4_cur) - 1.0_r8
     879           0 :            tmp_n1 = max( 0.0_r8, min( 1.0_r8, tmp_n1 ) )
     880           0 :            tmp_n2 = 1.0_r8 - tmp_n1
     881           0 :            tmp_n3 = 0.0_r8
     882             :         else
     883             : ! combination of ammonium bisulfate and sulfuric acid
     884             : ! tmp_n2 & tmp_n3 = mole fractions of the ammbisulf & sulfacid
     885   278834791 :            tmp_n1 = 0.0_r8
     886   278834791 :            tmp_n2 = (qnh3_cur/qh2so4_cur)
     887   278834791 :            tmp_n2 = max( 0.0_r8, min( 1.0_r8, tmp_n2 ) )
     888   278834791 :            tmp_n3 = 1.0_r8 - tmp_n2
     889             :         end if
     890             : 
     891   278834791 :         tmp_m1 = tmp_n1*mw_ammsulf
     892   278834791 :         tmp_m2 = tmp_n2*mw_ammbisulf
     893   278834791 :         tmp_m3 = tmp_n3*mw_sulfacid
     894             :         dens_part = (tmp_m1 + tmp_m2 + tmp_m3)/   &
     895             :            ((tmp_m1/dens_ammsulf) + (tmp_m2/dens_ammbisulf)   &
     896   278834791 :                                   + (tmp_m3/dens_sulfacid))
     897   278834791 :         dens_nh4so4a = dens_part
     898   278834791 :         mass_part  = voldry_part*dens_part 
     899             : ! (mol aerosol nh4)/(mol aerosol so4)
     900   278834791 :         molenh4a_per_moleso4a = 2.0_r8*tmp_n1 + tmp_n2  
     901             : ! (kg dry aerosol)/(mol aerosol so4)
     902   278834791 :         kgaero_per_moleso4a = 1.0e-3_r8*(tmp_m1 + tmp_m2 + tmp_m3)  
     903             : ! correction when host code sulfate is really ammonium bisulfate/sulfate
     904   278834791 :         kgaero_per_moleso4a = kgaero_per_moleso4a * (mw_so4a_host/mw_so4a)
     905             : 
     906             : ! fraction of wet volume due to so4a
     907   278834791 :         tmpb = 1.0_r8 + molenh4a_per_moleso4a*17.0_r8/98.0_r8
     908   278834791 :         wet_volfrac_so4a = 1.0_r8 / ( wetvol_dryvol * tmpb )
     909             : 
     910             : 
     911             : !
     912             : ! calc kerminen & kulmala (2002) correction
     913             : !
     914   278834791 :         if (igrow <=  0) then
     915           0 :             factor_kk = 1.0_r8
     916             : 
     917             :         else
     918             : ! "gr" parameter (nm/h) = condensation growth rate of new particles
     919             : ! use kk2002 eqn 21 for h2so4 uptake, and correct for nh3 & h2o uptake
     920   278834791 :             tmp_spd = 14.7_r8*sqrt(temp_in)   ! h2so4 molecular speed (m/s)
     921             :             gr_kk = 3.0e-9_r8*tmp_spd*mw_sulfacid*so4vol_in/   &
     922   278834791 :                     (dens_part*wet_volfrac_so4a)
     923             : 
     924             : ! "gamma" parameter (nm2/m2/h)
     925             : ! use kk2002 eqn 22
     926             : !
     927             : ! dfin_kk = wet diam (nm) of grown particle having dry dia = dpdry_part (m)
     928   278834791 :             dfin_kk = 1.0e9_r8 * dpdry_part * (wetvol_dryvol**onethird)
     929             : ! dnuc_kk = wet diam (nm) of cluster
     930   278834791 :             dnuc_kk = 2.0_r8*radius_cluster
     931   278834791 :             dnuc_kk = max( dnuc_kk, 1.0_r8 )
     932             : ! neglect (dmean/150)**0.048 factor, 
     933             : ! which should be very close to 1.0 because of small exponent
     934             :             gamma_kk = 0.23_r8 * (dnuc_kk)**0.2_r8   &
     935             :                      * (dfin_kk/3.0_r8)**0.075_r8   &
     936             :                      * (dens_part*1.0e-3_r8)**(-0.33_r8)   &
     937   278834791 :                      * (temp_in/293.0_r8)**(-0.75_r8)
     938             : 
     939             : ! "cs_prime parameter" (1/m2) 
     940             : ! instead kk2002 eqn 3, use
     941             : !     cs_prime ~= tmpa / (4*pi*tmpb * h2so4_accom_coef)
     942             : ! where
     943             : !     tmpa = -d(ln(h2so4))/dt by conden to particles   (1/h units)
     944             : !     tmpb = h2so4 vapor diffusivity (m2/h units)
     945             : ! this approx is generally within a few percent of the cs_prime
     946             : !     calculated directly from eqn 2, 
     947             : !     which is acceptable, given overall uncertainties
     948             : ! tmpa = -d(ln(h2so4))/dt by conden to particles   (1/h units)
     949   278834791 :             tmpa = h2so4_uptkrate * 3600.0_r8
     950   278834791 :             tmpa1 = tmpa
     951   278834791 :             tmpa = max( tmpa, 0.0_r8 )
     952             : ! tmpb = h2so4 gas diffusivity (m2/s, then m2/h)
     953   278834791 :             tmpb = 6.7037e-6_r8 * (temp_in**0.75_r8) / cair
     954   278834791 :             tmpb1 = tmpb         ! m2/s
     955   278834791 :             tmpb = tmpb*3600.0_r8   ! m2/h
     956   278834791 :             cs_prime_kk = tmpa/(4.0_r8*pi*tmpb*accom_coef_h2so4)
     957   278834791 :             cs_kk = cs_prime_kk*4.0_r8*pi*tmpb1
     958             : 
     959             : ! "nu" parameter (nm) -- kk2002 eqn 11
     960   278834791 :             nu_kk = gamma_kk*cs_prime_kk/gr_kk
     961             : ! nucleation rate adjustment factor (--) -- kk2002 eqn 13
     962   278834791 :             factor_kk = exp( (nu_kk/dfin_kk) - (nu_kk/dnuc_kk) )
     963             : 
     964             :         end if
     965   278834791 :         ratenuclt_kk = ratenuclt_bb*factor_kk
     966             : 
     967             : 
     968             : ! max production of aerosol dry mass (kg-aero/m3-air)
     969   278834791 :         tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc*mass_part) )
     970             : ! max production of aerosol so4 (mol-so4a/mol-air)
     971   278834791 :         tmpe = tmpa/(kgaero_per_moleso4a*cair)
     972             : ! max production of aerosol so4 (mol/mol-air)
     973             : ! based on ratenuclt_kk and mass_part
     974   278834791 :         qmolso4a_del_max = tmpe
     975             : 
     976             : ! check if max production exceeds available h2so4 vapor
     977   278834791 :         freducea = 1.0_r8
     978   278834791 :         if (qmolso4a_del_max .gt. qh2so4_cur) then
     979     7300673 :            freducea = qh2so4_cur/qmolso4a_del_max
     980             :         end if
     981             : 
     982             : ! check if max production exceeds available nh3 vapor
     983   278834791 :         freduceb = 1.0_r8
     984   278834791 :         if (molenh4a_per_moleso4a .ge. 1.0e-10_r8) then
     985             : ! max production of aerosol nh4 (ppm) based on ratenuclt_kk and mass_part
     986           0 :            qmolnh4a_del_max = qmolso4a_del_max*molenh4a_per_moleso4a
     987           0 :            if (qmolnh4a_del_max .gt. qnh3_cur) then
     988           0 :               freduceb = qnh3_cur/qmolnh4a_del_max
     989             :            end if
     990             :         end if
     991   278834791 :         freduce = min( freducea, freduceb )
     992             : 
     993             : ! if adjusted nucleation rate is less than 1e-12 #/m3/s ~= 0.1 #/cm3/day,
     994             : ! exit with new particle formation = 0
     995   278834791 :         if (freduce*ratenuclt_kk .le. 1.0e-12_r8) return
     996             : 
     997             : 
     998             : ! note:  suppose that at this point, freduce < 1.0 (no gas-available 
     999             : !    constraints) and molenh4a_per_moleso4a < 2.0
    1000             : ! if the gas-available constraints is do to h2so4 availability,
    1001             : !    then it would be possible to condense "additional" nh3 and have
    1002             : !    (nh3/h2so4 gas molar ratio) < (nh4/so4 aerosol molar ratio) <= 2 
    1003             : ! one could do some additional calculations of 
    1004             : !    dens_part & molenh4a_per_moleso4a to realize this
    1005             : ! however, the particle "growing" is a crude approximate way to get
    1006             : !    the new particles to the host code's minimum particle size,
    1007             : ! are such refinements worth the effort?
    1008             : 
    1009             : 
    1010             : ! changes to h2so4 & nh3 gas (in mol/mol-air), limited by amounts available
    1011   271270680 :         tmpa = 0.9999_r8
    1012   271270680 :         qh2so4_del = min( tmpa*qh2so4_cur, freduce*qmolso4a_del_max )
    1013   271270680 :         qnh3_del   = min( tmpa*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
    1014   271270680 :         qh2so4_del = -qh2so4_del
    1015   271270680 :         qnh3_del   = -qnh3_del
    1016             : 
    1017             : ! changes to so4 & nh4 aerosol (in mol/mol-air)
    1018   271270680 :         qso4a_del = -qh2so4_del
    1019   271270680 :         qnh4a_del =   -qnh3_del
    1020             : ! change to aerosol number (in #/mol-air)
    1021   271270680 :         qnuma_del = 1.0e-3_r8*(qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
    1022             : 
    1023             : ! do the following (tmpa, tmpb, tmpc) calculations as a check
    1024             : ! max production of aerosol number (#/mol-air)
    1025   271270680 :         tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc/cair) )
    1026             : ! adjusted production of aerosol number (#/mol-air)
    1027   271270680 :         tmpb = tmpa*freduce
    1028             : ! relative difference from qnuma_del
    1029   271270680 :         tmpc = (tmpb - qnuma_del)/max(tmpb, qnuma_del, 1.0e-35_r8)
    1030             : 
    1031             : 
    1032             : !
    1033             : ! diagnostic output to fort.41
    1034             : ! (this should be commented-out or deleted in the wrf-chem version)
    1035             : !
    1036   271270680 :         if (ldiagaa <= 0) return
    1037             : 
    1038           0 :         icase = icase + 1
    1039           0 :         if (abs(tmpc) .gt. abs(reldiffmax)) then
    1040           0 :            reldiffmax = tmpc
    1041           0 :            icase_reldiffmax = icase
    1042             :         end if
    1043             : !       do lun = 41, 51, 10
    1044           0 :         do lun = 6, 6
    1045             : !          write(lun,'(/)')
    1046             :            write(lun,'(a,2i9,1p,e10.2)')   &
    1047           0 :                'vehkam bin-nuc icase, icase_rdmax =',   &
    1048           0 :                icase, icase_reldiffmax, reldiffmax
    1049           0 :            if (freduceb .lt. freducea) then
    1050           0 :               if (abs(freducea-freduceb) .gt.   &
    1051             :                    3.0e-7_r8*max(freduceb,freducea)) write(lun,'(a,1p,2e15.7)')   &
    1052           0 :                  'freducea, b =', freducea, freduceb
    1053             :            end if
    1054             :         end do
    1055             : 
    1056             : ! output factors so that output matches that of ternucl03
    1057             : !       fogas  = 1.0e6                     ! convert mol/mol-air to ppm
    1058             : !       foso4a = 1.0e9*mw_so4a/mw_air      ! convert mol-so4a/mol-air to ug/kg-air
    1059             : !       fonh4a = 1.0e9*mw_nh4a/mw_air      ! convert mol-nh4a/mol-air to ug/kg-air
    1060             : !       fonuma = 1.0e3/mw_air              ! convert #/mol-air to #/kg-air
    1061             :         fogas  = 1.0_r8
    1062             :         foso4a = 1.0_r8
    1063             :         fonh4a = 1.0_r8
    1064             :         fonuma = 1.0_r8
    1065             : 
    1066             : !       do lun = 41, 51, 10
    1067           0 :         do lun = 6, 6
    1068             : 
    1069           0 :         write(lun,'(a,2i5)') 'newnuc_method_flagaa/aa2',   &
    1070           0 :            newnuc_method_flagaa, newnuc_method_flagaa2
    1071             : 
    1072           0 :         write(lun,9210)
    1073           0 :         write(lun,9201) temp_in, rh_in,   &
    1074           0 :            ratenuclt, 2.0_r8*radius_cluster*1.0e-7_r8, dpdry_part*1.0e2_r8,   &
    1075           0 :            voldry_part*1.0e6_r8, float(igrow)
    1076           0 :         write(lun,9215)
    1077             :         write(lun,9201)   &
    1078           0 :            qh2so4_avg*fogas, 0.0_r8,  &
    1079           0 :            qh2so4_cur*fogas, qnh3_cur*fogas,  &
    1080           0 :            qh2so4_del*fogas, qnh3_del*fogas,  &
    1081           0 :            qso4a_del*foso4a, qnh4a_del*fonh4a
    1082             : 
    1083           0 :         write(lun,9220)
    1084             :         write(lun,9201)   &
    1085           0 :            dtnuc, dens_nh4so4a*1.0e-3_r8,   &
    1086           0 :            (qnh3_cur/qh2so4_cur), molenh4a_per_moleso4a,   &
    1087           0 :            qnuma_del*fonuma, tmpb*fonuma, tmpc, freduce
    1088             : 
    1089             :         end do
    1090             : 
    1091             : !       lun = 51
    1092           0 :         lun = 6
    1093           0 :         write(lun,9230)
    1094             :         write(lun,9201)   &
    1095           0 :            press_in, cair*1.0e-6_r8, so4vol_in,   &
    1096           0 :            wet_volfrac_so4a, wetvol_dryvol, dens_part*1.0e-3_r8
    1097             : 
    1098           0 :         if (igrow > 0) then
    1099           0 :         write(lun,9240)
    1100             :         write(lun,9201)   &
    1101           0 :            tmp_spd, gr_kk, dnuc_kk, dfin_kk,   &
    1102           0 :            gamma_kk, tmpa1, tmpb1, cs_kk
    1103             : 
    1104           0 :         write(lun,9250)
    1105             :         write(lun,9201)   &
    1106           0 :            cs_prime_kk, nu_kk, factor_kk, ratenuclt,   &
    1107           0 :            ratenuclt_kk*1.0e-6_r8
    1108             :         end if
    1109             : 
    1110             : 9201    format ( 1p, 40e10.2  )
    1111             : 9210    format (   &
    1112             :         '      temp        rh',   &
    1113             :         '   ratenuc  dia_clus ddry_part',   &
    1114             :         ' vdry_part     igrow' )
    1115             : 9215    format (   &
    1116             :         '  h2so4avg  h2so4pre',   &
    1117             :         '  h2so4cur   nh3_cur',   &
    1118             :         '  h2so4del   nh3_del',   &
    1119             :         '  so4a_del  nh4a_del' )
    1120             : 9220    format (    &
    1121             :         '     dtnuc    dens_a   nh/so g   nh/so a',   &
    1122             :         '  numa_del  numa_dl2   reldiff   freduce' )
    1123             : 9230    format (   &
    1124             :         '  press_in      cair so4_volin',   &
    1125             :         ' wet_volfr wetv_dryv dens_part' )
    1126             : 9240    format (   &
    1127             :         '   tmp_spd     gr_kk   dnuc_kk   dfin_kk',   &
    1128             :         '  gamma_kk     tmpa1     tmpb1     cs_kk' )
    1129             : 9250    format (   &
    1130             :         ' cs_pri_kk     nu_kk factor_kk ratenuclt',   &
    1131             :         ' ratenu_kk' )
    1132             : 
    1133             : 
    1134             :         return
    1135             :         end subroutine mer07_veh02_nuc_mosaic_1box
    1136             : 
    1137             : 
    1138             : 
    1139             : !-----------------------------------------------------------------------
    1140             : !-----------------------------------------------------------------------
    1141    22327328 :         subroutine pbl_nuc_wang2008( so4vol,   &
    1142             :             newnuc_method_flagaa, newnuc_method_flagaa2,   &
    1143             :             ratenucl, rateloge,   &
    1144             :             cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
    1145             : !
    1146             : ! calculates boundary nucleation nucleation rate
    1147             : ! using the first or second-order parameterization in  
    1148             : !     wang, m., and j.e. penner, 2008,
    1149             : !        aerosol indirect forcing in a global model with particle nucleation,
    1150             : !        atmos. chem. phys. discuss., 8, 13943-13998
    1151             : !
    1152             :         implicit none
    1153             : 
    1154             : ! subr arguments (in)
    1155             :         real(r8), intent(in) :: so4vol            ! concentration of h2so4 (molecules cm-3)
    1156             :         integer, intent(in)  :: newnuc_method_flagaa  
    1157             :                                 ! [11,12] value selects [first,second]-order parameterization
    1158             : 
    1159             : ! subr arguments (inout)
    1160             :         integer, intent(inout)  :: newnuc_method_flagaa2
    1161             :         real(r8), intent(inout) :: ratenucl         ! binary nucleation rate, j (# cm-3 s-1)
    1162             :         real(r8), intent(inout) :: rateloge         ! log( ratenucl )
    1163             : 
    1164             :         real(r8), intent(inout) :: cnum_tot         ! total number of molecules
    1165             :                                                     ! in the critical nucleus
    1166             :         real(r8), intent(inout) :: cnum_h2so4       ! number of h2so4 molecules
    1167             :         real(r8), intent(inout) :: cnum_nh3         ! number of nh3 molecules
    1168             :         real(r8), intent(inout) :: radius_cluster   ! the radius of cluster (nm)
    1169             : 
    1170             : 
    1171             : ! local variables
    1172             :         real(r8) :: tmp_diam, tmp_mass, tmp_volu
    1173             :         real(r8) :: tmp_rateloge, tmp_ratenucl
    1174             : 
    1175             : ! executable
    1176             : 
    1177             : 
    1178             : ! nucleation rate
    1179    22327328 :         if (newnuc_method_flagaa == 11) then
    1180    22327328 :            tmp_ratenucl = 1.0e-6_r8 * so4vol
    1181           0 :         else if (newnuc_method_flagaa == 12) then
    1182           0 :            tmp_ratenucl = 1.0e-12_r8 * (so4vol**2)
    1183             :         else
    1184             :            return
    1185             :         end if
    1186    22327328 :         tmp_rateloge = log( tmp_ratenucl )
    1187             : 
    1188             : ! exit if pbl nuc rate is lower than (incoming) ternary/binary rate
    1189    22327328 :         if (tmp_rateloge <= rateloge) return
    1190             : 
    1191    22265208 :         rateloge = tmp_rateloge
    1192    22265208 :         ratenucl = tmp_ratenucl
    1193    22265208 :         newnuc_method_flagaa2 = newnuc_method_flagaa
    1194             : 
    1195             : ! following wang 2002, assume fresh nuclei are 1 nm diameter
    1196             : !    subsequent code will "grow" them to aitken mode size
    1197    22265208 :         radius_cluster = 0.5_r8
    1198             : 
    1199             : ! assume fresh nuclei are pure h2so4
    1200             : !    since aitken size >> initial size, the initial composition 
    1201             : !    has very little impact on the results
    1202    22265208 :         tmp_diam = radius_cluster * 2.0e-7_r8   ! diameter in cm
    1203    22265208 :         tmp_volu = (tmp_diam**3) * (pi/6.0_r8)  ! volume in cm^3
    1204    22265208 :         tmp_mass = tmp_volu * 1.8_r8            ! mass in g
    1205    22265208 :         cnum_h2so4 = (tmp_mass / 98.0_r8) * 6.023e23_r8   ! no. of h2so4 molec assuming pure h2so4
    1206    22265208 :         cnum_tot = cnum_h2so4
    1207    22265208 :         cnum_nh3 = 0.0_r8
    1208             : 
    1209             : 
    1210    22265208 :         return
    1211             :         end subroutine pbl_nuc_wang2008
    1212             : 
    1213             : 
    1214             : 
    1215             : !-----------------------------------------------------------------------
    1216             : !-----------------------------------------------------------------------
    1217   739780990 :         subroutine binary_nuc_vehk2002( temp, rh, so4vol,   &
    1218             :             ratenucl, rateloge,   &
    1219             :             cnum_h2so4, cnum_tot, radius_cluster )
    1220             : !
    1221             : ! calculates binary nucleation rate and critical cluster size
    1222             : ! using the parameterization in  
    1223             : !     vehkamäki, h., m. kulmala, i. napari, k.e.j. lehtinen,
    1224             : !        c. timmreck, m. noppel and a. laaksonen, 2002,
    1225             : !        an improved parameterization for sulfuric acid-water nucleation
    1226             : !        rates for tropospheric and stratospheric conditions,
    1227             : !        j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
    1228             : !
    1229             :         implicit none
    1230             : 
    1231             : ! subr arguments (in)
    1232             :         real(r8), intent(in) :: temp              ! temperature (k)  
    1233             :         real(r8), intent(in) :: rh                ! relative humidity (0-1)
    1234             :         real(r8), intent(in) :: so4vol            ! concentration of h2so4 (molecules cm-3)
    1235             : 
    1236             : ! subr arguments (out)
    1237             :         real(r8), intent(out) :: ratenucl         ! binary nucleation rate, j (# cm-3 s-1)
    1238             :         real(r8), intent(out) :: rateloge         ! log( ratenucl )
    1239             : 
    1240             :         real(r8), intent(out) :: cnum_h2so4       ! number of h2so4 molecules
    1241             :                                                   ! in the critical nucleus
    1242             :         real(r8), intent(out) :: cnum_tot         ! total number of molecules
    1243             :                                                   ! in the critical nucleus
    1244             :         real(r8), intent(out) :: radius_cluster   ! the radius of cluster (nm)
    1245             : 
    1246             : 
    1247             : ! local variables
    1248             :         real(r8) :: crit_x
    1249             :         real(r8) :: acoe, bcoe, ccoe, dcoe, ecoe, fcoe, gcoe, hcoe, icoe, jcoe
    1250             :         real(r8) :: tmpa, tmpb
    1251             : 
    1252             : ! executable
    1253             : 
    1254             : 
    1255             : ! calc sulfuric acid mole fraction in critical cluster
    1256             :         crit_x = 0.740997_r8 - 0.00266379_r8 * temp   &
    1257             :                - 0.00349998_r8 * log (so4vol)   &
    1258             :                + 0.0000504022_r8 * temp * log (so4vol)   &
    1259             :                + 0.00201048_r8 * log (rh)   &
    1260             :                - 0.000183289_r8 * temp * log (rh)   &
    1261             :                + 0.00157407_r8 * (log (rh)) ** 2.0_r8   &
    1262             :                - 0.0000179059_r8 * temp * (log (rh)) ** 2.0_r8   &
    1263             :                + 0.000184403_r8 * (log (rh)) ** 3.0_r8   &
    1264   739780990 :                - 1.50345e-6_r8 * temp * (log (rh)) ** 3.0_r8
    1265             : 
    1266             : 
    1267             : ! calc nucleation rate
    1268             :         acoe    = 0.14309_r8+2.21956_r8*temp   &
    1269             :                 - 0.0273911_r8 * temp**2.0_r8   &
    1270   739780990 :                 + 0.0000722811_r8 * temp**3.0_r8 + 5.91822_r8/crit_x
    1271             : 
    1272             :         bcoe    = 0.117489_r8 + 0.462532_r8 *temp   &
    1273             :                 - 0.0118059_r8 * temp**2.0_r8   &
    1274   739780990 :                 + 0.0000404196_r8 * temp**3.0_r8 + 15.7963_r8/crit_x
    1275             : 
    1276             :         ccoe    = -0.215554_r8-0.0810269_r8 * temp   &
    1277             :                 + 0.00143581_r8 * temp**2.0_r8   &
    1278             :                 - 4.7758e-6_r8 * temp**3.0_r8   &
    1279   739780990 :                 - 2.91297_r8/crit_x
    1280             : 
    1281             :         dcoe    = -3.58856_r8+0.049508_r8 * temp   &
    1282             :                 - 0.00021382_r8 * temp**2.0_r8   &
    1283             :                 + 3.10801e-7_r8 * temp**3.0_r8   &
    1284   739780990 :                 - 0.0293333_r8/crit_x
    1285             : 
    1286             :         ecoe    = 1.14598_r8 - 0.600796_r8 * temp   &
    1287             :                 + 0.00864245_r8 * temp**2.0_r8   &
    1288             :                 - 0.0000228947_r8 * temp**3.0_r8   &
    1289   739780990 :                 - 8.44985_r8/crit_x
    1290             : 
    1291             :         fcoe    = 2.15855_r8 + 0.0808121_r8 * temp   &
    1292             :                 -0.000407382_r8 * temp**2.0_r8   &
    1293             :                 -4.01957e-7_r8 * temp**3.0_r8   &
    1294   739780990 :                 + 0.721326_r8/crit_x
    1295             : 
    1296             :         gcoe    = 1.6241_r8 - 0.0160106_r8 * temp   &
    1297             :                 + 0.0000377124_r8 * temp**2.0_r8   &
    1298             :                 + 3.21794e-8_r8 * temp**3.0_r8   &
    1299   739780990 :                 - 0.0113255_r8/crit_x
    1300             : 
    1301             :         hcoe    = 9.71682_r8 - 0.115048_r8 * temp   &
    1302             :                 + 0.000157098_r8 * temp**2.0_r8   &
    1303             :                 + 4.00914e-7_r8 * temp**3.0_r8   &
    1304   739780990 :                 + 0.71186_r8/crit_x
    1305             : 
    1306             :         icoe    = -1.05611_r8 + 0.00903378_r8 * temp   &
    1307             :                 - 0.0000198417_r8 * temp**2.0_r8   &
    1308             :                 + 2.46048e-8_r8  * temp**3.0_r8   &
    1309   739780990 :                 - 0.0579087_r8/crit_x
    1310             : 
    1311             :         jcoe    = -0.148712_r8 + 0.00283508_r8 * temp   &
    1312             :                 - 9.24619e-6_r8  * temp**2.0_r8   &
    1313             :                 + 5.00427e-9_r8 * temp**3.0_r8   &
    1314   739780990 :                 - 0.0127081_r8/crit_x
    1315             : 
    1316             :         tmpa     =     (   &
    1317             :                   acoe   &
    1318             :                 + bcoe * log (rh)   &
    1319             :                 + ccoe * ( log (rh))**2.0_r8   &
    1320             :                 + dcoe * ( log (rh))**3.0_r8   &
    1321             :                 + ecoe * log (so4vol)   &
    1322             :                 + fcoe * (log (rh)) * (log (so4vol))   &
    1323             :                 + gcoe * ((log (rh) ) **2.0_r8)   &
    1324             :                        * (log (so4vol))   &
    1325             :                 + hcoe * (log (so4vol)) **2.0_r8   &
    1326             :                 + icoe * log (rh)   &
    1327             :                        * ((log (so4vol)) **2.0_r8)   &
    1328             :                 + jcoe * (log (so4vol)) **3.0_r8   &
    1329   739780990 :                 )
    1330   739780990 :         rateloge = tmpa
    1331   739780990 :         tmpa = min( tmpa, log(1.0e38_r8) )
    1332   739780990 :         ratenucl = exp ( tmpa )
    1333             : !       write(*,*) 'tmpa, ratenucl =', tmpa, ratenucl
    1334             : 
    1335             : 
    1336             : 
    1337             : ! calc number of molecules in critical cluster
    1338             :         acoe    = -0.00295413_r8 - 0.0976834_r8*temp   &
    1339             :                 + 0.00102485_r8 * temp**2.0_r8   &
    1340   739780990 :                 - 2.18646e-6_r8 * temp**3.0_r8 - 0.101717_r8/crit_x
    1341             : 
    1342             :         bcoe    = -0.00205064_r8 - 0.00758504_r8*temp   &
    1343             :                 + 0.000192654_r8 * temp**2.0_r8   &
    1344   739780990 :                 - 6.7043e-7_r8 * temp**3.0_r8 - 0.255774_r8/crit_x
    1345             : 
    1346             :         ccoe    = +0.00322308_r8 + 0.000852637_r8 * temp   &
    1347             :                 - 0.0000154757_r8 * temp**2.0_r8   &
    1348             :                 + 5.66661e-8_r8 * temp**3.0_r8   &
    1349   739780990 :                 + 0.0338444_r8/crit_x
    1350             : 
    1351             :         dcoe    = +0.0474323_r8 - 0.000625104_r8 * temp   &
    1352             :                 + 2.65066e-6_r8 * temp**2.0_r8   &
    1353             :                 - 3.67471e-9_r8 * temp**3.0_r8   &
    1354   739780990 :                 - 0.000267251_r8/crit_x
    1355             : 
    1356             :         ecoe    = -0.0125211_r8 + 0.00580655_r8 * temp   &
    1357             :                 - 0.000101674_r8 * temp**2.0_r8   &
    1358             :                 + 2.88195e-7_r8 * temp**3.0_r8   &
    1359   739780990 :                 + 0.0942243_r8/crit_x
    1360             : 
    1361             :         fcoe    = -0.038546_r8 - 0.000672316_r8 * temp   &
    1362             :                 + 2.60288e-6_r8 * temp**2.0_r8   &
    1363             :                 + 1.19416e-8_r8 * temp**3.0_r8   &
    1364   739780990 :                 - 0.00851515_r8/crit_x
    1365             : 
    1366             :         gcoe    = -0.0183749_r8 + 0.000172072_r8 * temp   &
    1367             :                 - 3.71766e-7_r8 * temp**2.0_r8   &
    1368             :                 - 5.14875e-10_r8 * temp**3.0_r8   &
    1369   739780990 :                 + 0.00026866_r8/crit_x
    1370             : 
    1371             :         hcoe    = -0.0619974_r8 + 0.000906958_r8 * temp   &
    1372             :                 - 9.11728e-7_r8 * temp**2.0_r8   &
    1373             :                 - 5.36796e-9_r8 * temp**3.0_r8   &
    1374   739780990 :                 - 0.00774234_r8/crit_x
    1375             : 
    1376             :         icoe    = +0.0121827_r8 - 0.00010665_r8 * temp   &
    1377             :                 + 2.5346e-7_r8 * temp**2.0_r8   &
    1378             :                 - 3.63519e-10_r8 * temp**3.0_r8   &
    1379   739780990 :                 + 0.000610065_r8/crit_x
    1380             : 
    1381             :         jcoe    = +0.000320184_r8 - 0.0000174762_r8 * temp   &
    1382             :                 + 6.06504e-8_r8 * temp**2.0_r8   &
    1383             :                 - 1.4177e-11_r8 * temp**3.0_r8   &
    1384   739780990 :                 + 0.000135751_r8/crit_x
    1385             : 
    1386             :         cnum_tot = exp (   &
    1387             :                   acoe   &
    1388             :                 + bcoe * log (rh)   &
    1389             :                 + ccoe * ( log (rh))**2.0_r8   &
    1390             :                 + dcoe * ( log (rh))**3.0_r8   &
    1391             :                 + ecoe * log (so4vol)   &
    1392             :                 + fcoe * (log (rh)) * (log (so4vol))   &
    1393             :                 + gcoe * ((log (rh) ) **2.0_r8)   &
    1394             :                        * (log (so4vol))   &
    1395             :                 + hcoe * (log (so4vol)) **2.0_r8   &
    1396             :                 + icoe * log (rh)   &
    1397             :                        * ((log (so4vol)) **2.0_r8)   &
    1398             :                 + jcoe * (log (so4vol)) **3.0_r8   &
    1399   739780990 :                 )
    1400             : 
    1401   739780990 :         cnum_h2so4 = cnum_tot * crit_x
    1402             : 
    1403             : !   calc radius (nm) of critical cluster
    1404             :         radius_cluster = exp( -1.6524245_r8 + 0.42316402_r8*crit_x   &
    1405   739780990 :                               + 0.3346648_r8*log(cnum_tot) )
    1406             :       
    1407             : 
    1408   739780990 :       return
    1409             :       end subroutine binary_nuc_vehk2002
    1410             : 
    1411             : 
    1412             : 
    1413             : !----------------------------------------------------------------------
    1414             : !----------------------------------------------------------------------
    1415        1536 : subroutine modal_aero_newnuc_init
    1416             : 
    1417             : !-----------------------------------------------------------------------
    1418             : !
    1419             : ! Purpose:
    1420             : !    set do_adjust and do_aitken flags
    1421             : !    create history fields for column tendencies associated with
    1422             : !       modal_aero_calcsize
    1423             : !
    1424             : ! Author: R. Easter
    1425             : !
    1426             : !-----------------------------------------------------------------------
    1427             : 
    1428             : use modal_aero_data
    1429             : use modal_aero_rename
    1430             : 
    1431             : use cam_abortutils,   only:  endrun
    1432             : use cam_history,      only:  addfld, add_default, fieldname_len, horiz_only
    1433             : use constituents,     only:  pcnst, cnst_get_ind, cnst_name
    1434             : use spmd_utils,       only:  masterproc
    1435             : use phys_control,     only: phys_getopts
    1436             : 
    1437             : 
    1438             : implicit none
    1439             : 
    1440             : !-----------------------------------------------------------------------
    1441             : ! arguments
    1442             : 
    1443             : !-----------------------------------------------------------------------
    1444             : ! local
    1445             :    integer  :: l_h2so4, l_nh3
    1446             :    integer  :: lnumait, lnh4ait, lso4ait
    1447             :    integer  :: l, l1, l2
    1448             :    integer  :: m, mait
    1449             : 
    1450             :    character(len=fieldname_len)   :: tmpname
    1451             :    character(len=fieldname_len+3) :: fieldname
    1452             :    character(128)                 :: long_name
    1453             :    character(8)                   :: unit
    1454             : 
    1455             :    logical                        :: dotend(pcnst)
    1456             :    logical                        :: history_aerosol      ! Output the MAM aerosol tendencies
    1457             : 
    1458             :    !-----------------------------------------------------------------------     
    1459             :    
    1460        1536 :         call phys_getopts( history_aerosol_out        = history_aerosol   )
    1461             : 
    1462             : 
    1463             : !   set these indices
    1464             : !   skip if no h2so4 species
    1465             : !   skip if no aitken mode so4 or num species
    1466        1536 :         l_h2so4_sv = 0
    1467        1536 :         l_nh3_sv = 0
    1468        1536 :         lnumait_sv = 0
    1469        1536 :         lnh4ait_sv = 0
    1470        1536 :         lso4ait_sv = 0
    1471             : 
    1472        1536 :         call cnst_get_ind( 'H2SO4', l_h2so4, .false. )
    1473        1536 :         call cnst_get_ind( 'NH3', l_nh3, .false. )
    1474             : 
    1475        1536 :         mait = modeptr_aitken
    1476        1536 :         if (mait > 0) then
    1477        1536 :             lnumait = numptr_amode(mait)
    1478        1536 :             lso4ait = lptr_so4_a_amode(mait)
    1479        1536 :             lnh4ait = lptr_nh4_a_amode(mait)
    1480             :         end if
    1481        1536 :         if ((l_h2so4  <= 0) .or. (l_h2so4 > pcnst)) then
    1482             :             write(*,'(/a/)')   &
    1483           0 :                 '*** modal_aero_newnuc bypass -- l_h2so4 <= 0'
    1484           0 :             return
    1485        1536 :         else if ((lso4ait <= 0) .or. (lso4ait > pcnst)) then
    1486             :             write(*,'(/a/)')   &
    1487           0 :                 '*** modal_aero_newnuc bypass -- lso4ait <= 0'
    1488           0 :             return
    1489        1536 :         else if ((lnumait <= 0) .or. (lnumait > pcnst)) then
    1490             :             write(*,'(/a/)')   &
    1491           0 :                 '*** modal_aero_newnuc bypass -- lnumait <= 0'
    1492           0 :             return
    1493        1536 :         else if ((mait <= 0) .or. (mait > ntot_amode)) then
    1494             :             write(*,'(/a/)')   &
    1495           0 :                 '*** modal_aero_newnuc bypass -- modeptr_aitken <= 0'
    1496           0 :             return
    1497             :         end if
    1498             : 
    1499        1536 :         l_h2so4_sv = l_h2so4
    1500        1536 :         l_nh3_sv   = l_nh3
    1501        1536 :         lnumait_sv = lnumait
    1502        1536 :         lnh4ait_sv = lnh4ait
    1503        1536 :         lso4ait_sv = lso4ait
    1504             : 
    1505             : !   set these constants
    1506             : !      mw_so4a_host is molec-wght of sulfate aerosol in host code
    1507             : !         96 when nh3/nh4 are simulated
    1508             : !         something else when nh3/nh4 are not simulated
    1509        1536 :         l = lptr_so4_a_amode(mait) ; l2 = -1
    1510        1536 :         if (l <= 0) call endrun( 'modal_aero_newnuch_init error a001 finding aitken so4' )
    1511        7680 :         do l1 = 1, nspec_amode(mait) 
    1512        7680 :            if (lmassptr_amode(l1,mait) == l) then
    1513        1536 :               l2 = l1
    1514        1536 :               mw_so4a_host = specmw_amode(l1,mait)
    1515        1536 :               dens_so4a_host = specdens_amode(l1,mait)
    1516             :            end if
    1517             :         end do
    1518        1536 :         if (l2 <= 0) call endrun( 'modal_aero_newnuch_init error a002 finding aitken so4' )
    1519             : 
    1520        1536 :         l = lptr_nh4_a_amode(mait) ; l2 = -1
    1521        1536 :         if (l > 0) then
    1522           0 :            do l1 = 1, nspec_amode(mait) 
    1523           0 :               if (lmassptr_amode(l1,mait) == l) then
    1524           0 :                  l2 = l1
    1525           0 :                  mw_nh4a_host = specmw_amode(l1,mait)
    1526             :               end if
    1527             :            end do
    1528           0 :            if (l2 <= 0) call endrun( 'modal_aero_newnuch_init error a002 finding aitken nh4' )
    1529             :         else
    1530        1536 :            mw_nh4a_host = mw_so4a_host
    1531             :         end if
    1532             : 
    1533             : !
    1534             : !   create history file column-tendency fields
    1535             : !
    1536        1536 :         dotend(:) = .false.
    1537        1536 :         dotend(lnumait) = .true.
    1538        1536 :         dotend(lso4ait) = .true.
    1539        1536 :         dotend(l_h2so4) = .true.
    1540             :         if ((l_nh3   > 0) .and. (l_nh3   <= pcnst) .and. &
    1541        1536 :             (lnh4ait > 0) .and. (lnh4ait <= pcnst)) then
    1542           0 :             dotend(lnh4ait) = .true.
    1543           0 :             dotend(l_nh3) = .true.
    1544             :         end if
    1545             : 
    1546       64512 :         do l = 1, pcnst
    1547       62976 :             if ( .not. dotend(l) ) cycle
    1548        4608 :             tmpname = cnst_name(l)
    1549        4608 :             unit = 'kg/m2/s'
    1550       23040 :             do m = 1, ntot_amode
    1551       23040 :                 if (l == numptr_amode(m)) unit = '#/m2/s'
    1552             :             end do
    1553        4608 :             fieldname = trim(tmpname) // '_sfnnuc1'
    1554        4608 :             long_name = trim(tmpname) // ' modal_aero new particle nucleation column tendency'
    1555        4608 :             call addfld( fieldname, horiz_only, 'A', unit, long_name )
    1556        4608 :             if ( history_aerosol ) then 
    1557           0 :                call add_default( fieldname, 1, ' ' )
    1558             :             endif
    1559        4608 :             if ( masterproc ) write(*,'(3(a,2x))') &
    1560        1542 :                 'modal_aero_newnuc_init addfld', fieldname, unit
    1561             :         end do ! l = ...
    1562             : 
    1563             : 
    1564             :       return
    1565        1536 :       end subroutine modal_aero_newnuc_init
    1566             : 
    1567             : 
    1568             : 
    1569             : !----------------------------------------------------------------------
    1570             : !----------------------------------------------------------------------
    1571           0 : subroutine ternary_nuc_merik2007( t, rh, c2, c3, j_log, ntot, nacid, namm, r )
    1572             : !subroutine ternary_fit(          t, rh, c2, c3, j_log, ntot, nacid, namm, r )
    1573             : ! *************************** ternary_fit.f90 ********************************
    1574             : ! joonas merikanto, 2006
    1575             : !
    1576             : ! fortran 90 subroutine that calculates the parameterized composition 
    1577             : ! and nucleation rate of critical clusters in h2o-h2so4-nh3 vapor
    1578             : !
    1579             : ! warning: the fit should not be used outside its limits of validity
    1580             : ! (limits indicated below)
    1581             : !
    1582             : ! in:
    1583             : ! t:     temperature (k), limits 235-295 k
    1584             : ! rh:    relative humidity as fraction (eg. 0.5=50%) limits 0.05-0.95
    1585             : ! c2:    sulfuric acid concentration (molecules/cm3) limits 5x10^4 - 10^9 molecules/cm3
    1586             : ! c3:    ammonia mixing ratio (ppt) limits 0.1 - 1000 ppt
    1587             : !
    1588             : ! out:
    1589             : ! j_log: logarithm of nucleation rate (1/(s cm3))
    1590             : ! ntot:  total number of molecules in the critical cluster
    1591             : ! nacid: number of sulfuric acid molecules in the critical cluster
    1592             : ! namm:  number of ammonia molecules in the critical cluster
    1593             : ! r:     radius of the critical cluster (nm)
    1594             : !  ****************************************************************************
    1595             : implicit none
    1596             : 
    1597             : real(r8), intent(in) :: t, rh, c2, c3
    1598             : real(r8), intent(out) :: j_log, ntot, nacid, namm, r
    1599             : real(r8) :: j, t_onset
    1600             : 
    1601             : t_onset=143.6002929064716_r8 + 1.0178856665693992_r8*rh + &
    1602             :    10.196398812974294_r8*log(c2) - &
    1603             :    0.1849879416839113_r8*log(c2)**2 - 17.161783213150173_r8*log(c3) + &
    1604             :    (109.92469248546053_r8*log(c3))/log(c2) + &
    1605           0 :    0.7734119613144357_r8*log(c2)*log(c3) - 0.15576469879527022_r8*log(c3)**2
    1606             : 
    1607           0 : if(t_onset.gt.t) then 
    1608             : 
    1609             :    j_log=-12.861848898625231_r8 + 4.905527742256349_r8*c3 - 358.2337705052991_r8*rh -& 
    1610             :    0.05463019231872484_r8*c3*t + 4.8630382337426985_r8*rh*t + &
    1611             :    0.00020258394697064567_r8*c3*t**2 - 0.02175548069741675_r8*rh*t**2 - &
    1612             :    2.502406532869512e-7_r8*c3*t**3 + 0.00003212869941055865_r8*rh*t**3 - &
    1613             :    4.39129415725234e6_r8/log(c2)**2 + (56383.93843154586_r8*t)/log(c2)**2 -& 
    1614             :    (239.835990963361_r8*t**2)/log(c2)**2 + &
    1615             :    (0.33765136625580167_r8*t**3)/log(c2)**2 - &
    1616             :    (629.7882041830943_r8*rh)/(c3**3*log(c2)) + &
    1617             :    (7.772806552631709_r8*rh*t)/(c3**3*log(c2)) - &
    1618             :    (0.031974053936299256_r8*rh*t**2)/(c3**3*log(c2)) + &
    1619             :    (0.00004383764128775082_r8*rh*t**3)/(c3**3*log(c2)) + &
    1620             :    1200.472096232311_r8*log(c2) - 17.37107890065621_r8*t*log(c2) + &
    1621             :    0.08170681335921742_r8*t**2*log(c2) - 0.00012534476159729881_r8*t**3*log(c2) - &
    1622             :    14.833042158178936_r8*log(c2)**2 + 0.2932631303555295_r8*t*log(c2)**2 - &
    1623             :    0.0016497524241142845_r8*t**2*log(c2)**2 + &
    1624             :    2.844074805239367e-6_r8*t**3*log(c2)**2 - 231375.56676032578_r8*log(c3) - &
    1625             :    100.21645273730675_r8*rh*log(c3) + 2919.2852552424706_r8*t*log(c3) + &
    1626             :    0.977886555834732_r8*rh*t*log(c3) - 12.286497122264588_r8*t**2*log(c3) - &
    1627             :    0.0030511783284506377_r8*rh*t**2*log(c3) + &
    1628             :    0.017249301826661612_r8*t**3*log(c3) + 2.967320346100855e-6_r8*rh*t**3*log(c3) + &
    1629             :    (2.360931724951942e6_r8*log(c3))/log(c2) - &
    1630             :    (29752.130254319443_r8*t*log(c3))/log(c2) + &
    1631             :    (125.04965118142027_r8*t**2*log(c3))/log(c2) - &
    1632             :    (0.1752996881934318_r8*t**3*log(c3))/log(c2) + &
    1633             :    5599.912337254629_r8*log(c2)*log(c3) - 70.70896612937771_r8*t*log(c2)*log(c3) + &
    1634             :    0.2978801613269466_r8*t**2*log(c2)*log(c3) - &
    1635             :    0.00041866525019504_r8*t**3*log(c2)*log(c3) + 75061.15281456841_r8*log(c3)**2 - &
    1636             :    931.8802278173565_r8*t*log(c3)**2 + 3.863266220840964_r8*t**2*log(c3)**2 - &
    1637             :    0.005349472062284983_r8*t**3*log(c3)**2 - &
    1638             :    (732006.8180571689_r8*log(c3)**2)/log(c2) + &
    1639             :    (9100.06398573816_r8*t*log(c3)**2)/log(c2) - &
    1640             :    (37.771091915932004_r8*t**2*log(c3)**2)/log(c2) + &
    1641             :    (0.05235455395566905_r8*t**3*log(c3)**2)/log(c2) - &
    1642             :    1911.0303773001353_r8*log(c2)*log(c3)**2 + &
    1643             :    23.6903969622286_r8*t*log(c2)*log(c3)**2 - &
    1644             :    0.09807872005428583_r8*t**2*log(c2)*log(c3)**2 + &
    1645             :    0.00013564560238552576_r8*t**3*log(c2)*log(c3)**2 - &
    1646             :    3180.5610833308_r8*log(c3)**3 + 39.08268568672095_r8*t*log(c3)**3 - &
    1647             :    0.16048521066690752_r8*t**2*log(c3)**3 + &
    1648             :    0.00022031380023793877_r8*t**3*log(c3)**3 + &
    1649             :    (40751.075322248245_r8*log(c3)**3)/log(c2) - &
    1650             :    (501.66977622013934_r8*t*log(c3)**3)/log(c2) + &
    1651             :    (2.063469732254135_r8*t**2*log(c3)**3)/log(c2) - &
    1652             :    (0.002836873785758324_r8*t**3*log(c3)**3)/log(c2) + &
    1653             :    2.792313345723013_r8*log(c2)**2*log(c3)**3 - &
    1654             :    0.03422552111802899_r8*t*log(c2)**2*log(c3)**3 + &
    1655             :    0.00014019195277521142_r8*t**2*log(c2)**2*log(c3)**3 - &
    1656             :    1.9201227328396297e-7_r8*t**3*log(c2)**2*log(c3)**3 - &
    1657             :    980.923146020468_r8*log(rh) + 10.054155220444462_r8*t*log(rh) - &
    1658             :    0.03306644502023841_r8*t**2*log(rh) + 0.000034274041225891804_r8*t**3*log(rh) + &
    1659             :    (16597.75554295064_r8*log(rh))/log(c2) - &
    1660             :    (175.2365504237746_r8*t*log(rh))/log(c2) + &
    1661             :    (0.6033215603167458_r8*t**2*log(rh))/log(c2) - &
    1662             :    (0.0006731787599587544_r8*t**3*log(rh))/log(c2) - &
    1663             :    89.38961120336789_r8*log(c3)*log(rh) + 1.153344219304926_r8*t*log(c3)*log(rh) - &
    1664             :    0.004954549700267233_r8*t**2*log(c3)*log(rh) + &
    1665             :    7.096309866238719e-6_r8*t**3*log(c3)*log(rh) + &
    1666             :    3.1712136610383244_r8*log(c3)**3*log(rh) - &
    1667             :    0.037822330602328806_r8*t*log(c3)**3*log(rh) + &
    1668             :    0.0001500555743561457_r8*t**2*log(c3)**3*log(rh) - &
    1669           0 :    1.9828365865570703e-7_r8*t**3*log(c3)**3*log(rh)
    1670             : 
    1671           0 :    j=exp(j_log)
    1672             : 
    1673             :    ntot=57.40091052369212_r8 - 0.2996341884645408_r8*t + &
    1674             :    0.0007395477768531926_r8*t**2 - &
    1675             :    5.090604835032423_r8*log(c2) + 0.011016634044531128_r8*t*log(c2) + &
    1676             :    0.06750032251225707_r8*log(c2)**2 - 0.8102831333223962_r8*log(c3) + &
    1677             :    0.015905081275952426_r8*t*log(c3) - 0.2044174683159531_r8*log(c2)*log(c3) + &
    1678             :    0.08918159167625832_r8*log(c3)**2 - 0.0004969033586666147_r8*t*log(c3)**2 + &
    1679             :    0.005704394549007816_r8*log(c3)**3 + 3.4098703903474368_r8*log(j) - &
    1680             :    0.014916956508210809_r8*t*log(j) + 0.08459090011666293_r8*log(c3)*log(j) - &
    1681           0 :    0.00014800625143907616_r8*t*log(c3)*log(j) + 0.00503804694656905_r8*log(j)**2
    1682             :  
    1683             :    r=3.2888553966535506e-10_r8 - 3.374171768439839e-12_r8*t + &
    1684             :    1.8347359507774313e-14_r8*t**2 + 2.5419844298881856e-12_r8*log(c2) - &
    1685             :    9.498107643050827e-14_r8*t*log(c2) + 7.446266520834559e-13_r8*log(c2)**2 + &
    1686             :    2.4303397746137294e-11_r8*log(c3) + 1.589324325956633e-14_r8*t*log(c3) - &
    1687             :    2.034596219775266e-12_r8*log(c2)*log(c3) - 5.59303954457172e-13_r8*log(c3)**2 - &
    1688             :    4.889507104645867e-16_r8*t*log(c3)**2 + 1.3847024107506764e-13_r8*log(c3)**3 + &
    1689             :    4.141077193427042e-15_r8*log(j) - 2.6813110884009767e-14_r8*t*log(j) + &
    1690             :    1.2879071621313094e-12_r8*log(c3)*log(j) - &
    1691           0 :    3.80352446061867e-15_r8*t*log(c3)*log(j) - 1.8790172502456827e-14_r8*log(j)**2
    1692             :  
    1693             :    nacid=-4.7154180661803595_r8 + 0.13436423483953885_r8*t - & 
    1694             :    0.00047184686478816176_r8*t**2 - & 
    1695             :    2.564010713640308_r8*log(c2) + 0.011353312899114723_r8*t*log(c2) + &
    1696             :    0.0010801941974317014_r8*log(c2)**2 + 0.5171368624197119_r8*log(c3) - &
    1697             :    0.0027882479896204665_r8*t*log(c3) + 0.8066971907026886_r8*log(c3)**2 - & 
    1698             :    0.0031849094214409335_r8*t*log(c3)**2 - 0.09951184152927882_r8*log(c3)**3 + &
    1699             :    0.00040072788891745513_r8*t*log(c3)**3 + 1.3276469271073974_r8*log(j) - &
    1700             :    0.006167654171986281_r8*t*log(j) - 0.11061390967822708_r8*log(c3)*log(j) + &
    1701           0 :    0.0004367575329273496_r8*t*log(c3)*log(j) + 0.000916366357266258_r8*log(j)**2
    1702             :  
    1703             :    namm=71.20073903979772_r8 - 0.8409600103431923_r8*t + &
    1704             :    0.0024803006590334922_r8*t**2 + &
    1705             :    2.7798606841602607_r8*log(c2) - 0.01475023348171676_r8*t*log(c2) + &
    1706             :    0.012264508212031405_r8*log(c2)**2 - 2.009926050440182_r8*log(c3) + &
    1707             :    0.008689123511431527_r8*t*log(c3) - 0.009141180198955415_r8*log(c2)*log(c3) + &
    1708             :    0.1374122553905617_r8*log(c3)**2 - 0.0006253227821679215_r8*t*log(c3)**2 + &
    1709             :    0.00009377332742098946_r8*log(c3)**3 + 0.5202974341687757_r8*log(j) - &
    1710             :    0.002419872323052805_r8*t*log(j) + 0.07916392322884074_r8*log(c3)*log(j) - &
    1711           0 :    0.0003021586030317366_r8*t*log(c3)*log(j) + 0.0046977006608603395_r8*log(j)**2
    1712             : 
    1713             : else
    1714             : ! nucleation rate less that 5e-6, setting j_log arbitrary small
    1715           0 :    j_log=-300._r8
    1716             : end if
    1717             : 
    1718           0 : return
    1719             : 
    1720        1536 : end  subroutine ternary_nuc_merik2007
    1721             : 
    1722             : !----------------------------------------------------------------------
    1723             : end module modal_aero_newnuc
    1724             : 
    1725             : 
    1726             :