chem/module_mosaic_astem.F

   1   module module_mosaic_astem
   2
   3
   4   use module_mosaic_support, only: mosaic_warn_mess, mosaic_err_mess
   5   use module_data_mosaic_kind, only: r8
   6
   7   implicit none
   8
   9
  10   contains
  11
  12
  13 ! feb 22. new flux_mix
  14
  15 !***********************************************************************
  16 ! ASTEM: Adaptive Step Time-Split Euler Method
  17 !
  18 ! author: Rahul A. Zaveri
  19 ! update: jan 2007
  20 !-----------------------------------------------------------------------
  21
  22   subroutine ASTEM(   mcall_print_aer,                                               &!intent-ins
  23        dtchem,        sigmag_a,  aH2O,   T_K,          RH_pc,     P_atm,             &
  24        kappa_nonelectro,                                                             &
  25        jaerosolstate, flux_s,            flux_l,       volatile_s,iprint_input,      &!intent -inout
  26        phi_volatile_s,phi_volatile_l,    jphase,       aer,       kg,       gas,     &
  27        gas_avg,       gas_netprod_otrproc,                                           &
  28        jhyst_leg,     electrolyte,       epercent,     activity,  mc,       sat_soa, &
  29        num_a,         Dp_dry_a,          Dp_wet_a,     dp_core_a, mass_dry_a,        &
  30        mass_soluble_a,vol_dry_a,         dens_dry_a,   water_a,   water_a_hyst,      &
  31        water_a_up,    aH2O_a,            total_species,tot_cl_in, ma,       gam,     &
  32        log_gamZ,      gam_ratio,         Keq_ll,       Keq_gl,    Keq_sg,   Kp_nh4cl,&
  33        Kp_nh4no3,     sigma_water,       Keq_sl,       MDRH_T,    molality0,         &
  34        uptkrate_h2so4,                   mosaic_vars_aa,                             &
  35        area_dry_a,    area_wet_a,        mass_wet_a,vol_wet_a,                       &!intent-out
  36        dens_wet_a,    ri_shell_a,        ri_avg_a,     ri_core_a                     )
  37
  38   use module_data_mosaic_aero, only: r8, nbin_a_max,                               &
  39        ngas_aerchtot, ngas_volatile, nelectrolyte,                                 &
  40        Ncation, naer, mYES, no_aerosol, Nanion, nrxn_aer_gl, nrxn_aer_ll,          &
  41        nrxn_aer_sg, nrxn_aer_sl, naercomp, nsalt, MDRH_T_NUM, jsulf_poor_NUM,      &
  42        jsulf_rich_NUM, nbin_a, zc, za,                                             &
  43        a_zsr, b_zsr, mw_electrolyte, partial_molar_vol, mw_aer_mac, dens_aer_mac,  &
  44        MW_a, MW_c, dens_comp_a, mw_comp_a, ref_index_a, rtol_mesa, jsalt_index,    &
  45        jsulf_poor, jsulf_rich, ih2so4_g,                                           &
  46        iso4_a, jtotal, msoa_flag1,                                                 & !for debug only remove it later BALLI
  47        mosaic_vars_aa_type
  48
  49   use module_mosaic_ext, only: aerosol_phase_state,calc_dry_n_wet_aerosol_props,   &
  50        aerosolmtc, dumpxx
  51
  52   use module_mosaic_soa_vbs, only: mosaic_soa_vbs_intr
  53
  54 ! use module_print_aer,  only: print_aer
  55
  56
  57
  58   !Subroutine Arguments
  59   !Intent-ins
  60   integer, intent(in) :: mcall_print_aer
  61
  62   real(r8), intent(in) :: dtchem
  63   real(r8), intent(in) :: aH2O
  64   real(r8), intent(in) :: T_K, RH_pc, P_atm
  65   real(r8), intent(in), dimension(nbin_a_max) :: sigmag_a
  66   real(r8), intent(in), dimension(ngas_aerchtot) :: gas_netprod_otrproc
  67   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
  68
  69   !intent-inouts
  70   integer, intent(inout) :: iprint_input
  71
  72   integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate, jphase
  73   integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
  74
  75   real(r8), intent(inout) :: tot_cl_in
  76   real(r8), intent(inout) :: Kp_nh4cl, Kp_nh4no3
  77   real(r8), intent(inout) :: sigma_water
  78
  79   real(r8), intent(inout), dimension(nbin_a_max)     :: num_a, Dp_dry_a, Dp_wet_a
  80   real(r8), intent(inout), dimension(nbin_a_max)     :: dp_core_a
  81   real(r8), intent(inout), dimension(nbin_a_max)     :: mass_dry_a, mass_soluble_a
  82   real(r8), intent(inout), dimension(nbin_a_max)     :: vol_dry_a, dens_dry_a
  83   real(r8), intent(inout), dimension(nbin_a_max)     :: water_a
  84   real(r8), intent(inout), dimension(nbin_a_max)     :: water_a_hyst,water_a_up
  85   real(r8), intent(inout), dimension(nbin_a_max)     :: aH2O_a
  86   real(r8), intent(inout), dimension(nbin_a_max)     :: gam_ratio
  87   real(r8), intent(inout), dimension(ngas_aerchtot)  :: gas
  88   real(r8), intent(inout), dimension(ngas_aerchtot)  :: gas_avg  ! average gas conc. over dtchem time step (nmol/m3)
  89   real(r8), intent(inout)                            :: uptkrate_h2so4  ! rate of h2so4 uptake by aerosols (1/s)
  90
  91   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
  92
  93   ! gas_netprod_otrproc = gas net production rate from other processes
  94   !    such as gas-phase chemistry and emissions (nmol/m3/s)
  95   ! this allows the condensation (gasaerexch) routine to apply production and condensation loss
  96   !    together, which is more accurate numerically
  97   ! NOTE - currently for mosaic, only the value for h2so4 can be non-zero
  98   real(r8), intent(inout), dimension(ngas_volatile)  :: sat_soa
  99   real(r8), intent(inout), dimension(ngas_volatile)  :: total_species
 100   real(r8), intent(inout), dimension(nrxn_aer_ll)    :: Keq_ll
 101   real(r8), intent(inout), dimension(nrxn_aer_gl)    :: Keq_gl
 102   real(r8), intent(inout), dimension(nrxn_aer_sg)    :: Keq_sg
 103   real(r8), intent(inout), dimension(nrxn_aer_sl)    :: Keq_sl
 104   real(r8), intent(inout), dimension(MDRH_T_NUM)     :: MDRH_T
 105   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)   :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
 106
 107   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max)  :: flux_s,flux_l
 108   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max)  :: volatile_s
 109   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max)  :: phi_volatile_s
 110   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max)  :: phi_volatile_l
 111   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max)  :: kg
 112   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)   :: activity
 113   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)   :: gam
 114   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
 115   real(r8), intent(inout), dimension(Ncation,nbin_a_max)        :: mc
 116   real(r8), intent(inout), dimension(Nanion,nbin_a_max)         :: ma
 117
 118   real(r8), intent(inout), dimension(naer,3,nbin_a_max)         :: aer
 119   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 120   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
 121
 122
 123   !Intent-out
 124   real(r8), intent(out), dimension(nbin_a_max) :: area_dry_a
 125   real(r8), intent(out), dimension(nbin_a_max) :: area_wet_a,mass_wet_a
 126   real(r8), intent(out), dimension(nbin_a_max) :: vol_wet_a
 127   real(r8), intent(out), dimension(nbin_a_max) :: dens_wet_a
 128
 129   complex,  intent(out), dimension(nbin_a_max) :: ri_shell_a,ri_avg_a,ri_core_a
 130
 131   !Local variables
 132   integer :: ibin, iv, itmpa
 133
 134   integer, dimension(nsalt) :: jsalt_present
 135   integer, dimension(ngas_volatile,3,nbin_a_max) :: integrate
 136
 137
 138   real(r8) :: Keq_nh4cl
 139   real(r8) :: swdown_cell
 140
 141   real(r8), dimension(nsalt) :: phi_salt_old
 142   real(r8), dimension(Ncation) :: nc_Mc,xeq_c
 143   real(r8), dimension(Nanion)  :: xeq_a,na_Ma
 144   real(r8), dimension(nbin_a_max) :: sigma_soln
 145   real(r8), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
 146   real(r8), dimension(nbin_a_max) :: delta_hcl_max
 147   real(r8), dimension(nbin_a_max) :: growth_factor
 148   real(r8), dimension(nbin_a_max) :: MDRH
 149
 150   real(r8), dimension(ngas_volatile) ::sfc_a
 151   real(r8), dimension(ngas_volatile,nbin_a_max) :: Heff
 152   real(r8), dimension(ngas_aerchtot,nbin_a_max) :: kel
 153
 154
 155   call dumpxx( 'cc', dtchem, t_k, p_atm, ah2o, &
 156          jaerosolstate, jphase, jhyst_leg, &
 157          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 158          mosaic_vars_aa )
 159
 160   mosaic_vars_aa%niter_MESA     = 0.0_r8
 161   mosaic_vars_aa%niter_MESA_max = 0
 162   mosaic_vars_aa%jMESA_fail     = 0
 163   mosaic_vars_aa%jMESA_call     = 0
 164   mosaic_vars_aa%iter_MESA(:)   = 0
 165
 166   phi_salt_old(:)  = 0.0_r8
 167   integrate(:,:,:) = 0.0_r8 !BALLI- Ask Dick about this initialization
 168   heff(:,:)        = 0.0_r8 !BALLI- Ask Dick about this initialization
 169
 170   gas_avg(:) = gas(:)  ! RCE:  set avg. gas conc. = initial conc.
 171
 172   ! update ASTEM call counter
 173   mosaic_vars_aa%jASTEM_call  = mosaic_vars_aa%jASTEM_call + 1
 174
 175   ! reset input print flag
 176   iprint_input = mYES
 177
 178   ! compute aerosol phase state before starting integration
 179   do ibin = 1, nbin_a
 180       area_dry_a(ibin) = 0.0_r8 !BSINGH - Ask Dick about it. The code blows up in print_aer
 181       area_wet_a(ibin) = 0.0_r8 !BSINGH - Ask Dick about it. The code blows up in print_aer
 182       mass_wet_a(ibin) = 0.0_r8 !BSINGH - Ask Dick about it. The code blows up in print_aer
 183      if(jaerosolstate(ibin) .ne. no_aerosol)then
 184         call aerosol_phase_state( ibin, jaerosolstate, jphase,  &
 185              aer, jhyst_leg, electrolyte, epercent, kel, activity, mc, num_a, mass_wet_a, &
 186              mass_dry_a, mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst,  &
 187              water_a_up, aH2O_a, aH2O, ma, gam, & !BALLI
 188              log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c,              & ! RAZ deleted a_zsr
 189              mw_electrolyte, partial_molar_vol, sigma_soln, T_K, RH_pc, mw_aer_mac,    &
 190              dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a, MW_c, growth_factor, MDRH, &
 191              MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index, jsulf_poor,     &
 192              jsulf_rich, phi_salt_old,                                      &
 193              kappa_nonelectro, mosaic_vars_aa )
 194
 195         call calc_dry_n_wet_aerosol_props(                          &
 196            ibin, jaerosolstate, aer, electrolyte, water_a, num_a,         &  ! input
 197            dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &  ! input
 198            Dp_dry_a, Dp_wet_a, dp_core_a,                                 &  ! output
 199            area_dry_a, area_wet_a, mass_dry_a, mass_wet_a,                &  ! output
 200            vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a,                  &  ! output
 201            ri_shell_a, ri_core_a, ri_avg_a                                )  ! output
 202      endif
 203   enddo
 204   call check_astem_negative( 1, mosaic_vars_aa%xnerr_astem_negative, &
 205                                 mosaic_vars_aa%fix_astem_negative, aer, gas )
 206
 207   ! BOX
 208   if (mcall_print_aer == 2) then
 209      !call print_aer(0,jaerosolstate,isteps_ASTEM,iter_MESA,aer,gas,electrolyte,  &
 210      !     mc,num_a,Dp_dry_a,Dp_wet_a,area_dry_a,area_wet_a,mass_wet_a,mass_dry_a,&
 211      !     water_a)     ! UNCOMMENT THIS LINE
 212   endif                 ! UNCOMMENT THIS LINE
 213
 214
 215   ! compute new gas-aerosol mass transfer coefficients
 216   call aerosolmtc( jaerosolstate, num_a, Dp_wet_a, sigmag_a, P_atm, T_K, kg )
 217
 218   uptkrate_h2so4 = sum( kg(ih2so4_g,1:nbin_a) )
 219
 220   call dumpxx( 'ee', dtchem, t_k, p_atm, ah2o, &
 221          jaerosolstate, jphase, jhyst_leg, &
 222          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 223          mosaic_vars_aa )
 224
 225   ! condense h2so4, msa, and nh3 only
 226   call ASTEM_non_volatiles( dtchem, jaerosolstate, jphase, aer,  &
 227        kg, gas, gas_avg, gas_netprod_otrproc,                                        &
 228        jhyst_leg, electrolyte, epercent, kel, activity, mc, delta_nh3_max,              &
 229        delta_hno3_max, delta_hcl_max, num_a, mass_wet_a, mass_dry_a, mass_soluble_a,   &
 230        vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, total_species,  &
 231        tot_cl_in,                                                                 &
 232        aH2O, ma, gam, log_gamZ, zc, za,          &
 233        gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, a_zsr, mw_electrolyte, partial_molar_vol,  &
 234        sigma_soln, T_K, RH_pc, mw_aer_mac, dens_aer_mac, sigma_water, Keq_ll, Keq_sl,    &
 235        MW_a, MW_c, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present,     &
 236        jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old,                     &
 237        kappa_nonelectro, mosaic_vars_aa )       ! analytical solution
 238
 239   call dumpxx( 'gg', dtchem, t_k, p_atm, ah2o, &
 240          jaerosolstate, jphase, jhyst_leg, &
 241          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 242          mosaic_vars_aa )
 243
 244   call check_astem_negative( 2, mosaic_vars_aa%xnerr_astem_negative, &
 245                                 mosaic_vars_aa%fix_astem_negative, aer, gas )
 246
 247   ! condense inorganic semi-volatile gases hno3, hcl, nh3, and co2
 248   call ASTEM_semi_volatiles( iprint_input, dtchem, jaerosolstate,                     &
 249        sfc_a, flux_s, flux_l, Heff, volatile_s, phi_volatile_s,   &
 250        jphase, aer, kg, gas, jhyst_leg, electrolyte, epercent, kel, activity, mc,          &
 251        delta_nh3_max, delta_hno3_max, delta_hcl_max,                                 &
 252        num_a, mass_dry_a, mass_wet_a, mass_soluble_a,   &
 253        vol_dry_a, vol_wet_a, water_a, total_species, tot_cl_in,                       &
 254        aH2O_a, aH2O, ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte,    &
 255        Keq_ll, Keq_gl, Keq_sg, Kp_nh4cl, Kp_nh4no3, Keq_nh4cl, MW_c, MW_a, mw_aer_mac,    &
 256        dens_aer_mac, Keq_sl, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa,         &
 257        jsalt_present, jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old,    &
 258        integrate, phi_volatile_l,                                                      &
 259        kappa_nonelectro, mosaic_vars_aa )       ! semi-implicit + explicit euler
 260
 261   call dumpxx( 'ii', dtchem, t_k, p_atm, ah2o, &
 262          jaerosolstate, jphase, jhyst_leg, &
 263          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 264          mosaic_vars_aa )
 265
 266   if (mosaic_vars_aa%f_mos_fail > 0 ) then
 267      return
 268   endif
 269   call check_astem_negative( 3, mosaic_vars_aa%xnerr_astem_negative, &
 270                                 mosaic_vars_aa%fix_astem_negative, aer, gas )
 271
 272   ! condense secondary organic gases (8 sorgam species)
 273
 274   if (msoa_flag1 == 1) then
 275      itmpa = 1
 276      call ASTEM_secondary_organics(dtchem,jaerosolstate,sfc_a,Heff,phi_volatile_l,  &
 277         integrate,aer,kg,gas,sat_soa,total_species) ! semi-implicit euler
 278
 279   else if (msoa_flag1 == 1000) then
 280      itmpa = 2
 281      swdown_cell = mosaic_vars_aa%swdown
 282      call mosaic_soa_vbs_intr( &
 283         dtchem, p_atm, t_k, swdown_cell, &
 284         jaerosolstate, &
 285         aer, gas, water_a, area_wet_a, dp_wet_a, &
 286         kg, sat_soa, total_species, &
 287         ma, mc, mosaic_vars_aa )
 288
 289   else
 290      itmpa = 0
 291   end if
 292
 293   if (itmpa > 0) then
 294   call check_astem_negative( 4, mosaic_vars_aa%xnerr_astem_negative, &
 295                                 mosaic_vars_aa%fix_astem_negative, aer, gas )
 296   end if
 297
 298   do iv = 1, ngas_aerchtot
 299      if (iv == ih2so4_g) cycle
 300      ! RCE:  avg. gas conc. = 0.5*( initial conc. + current conc. )
 301      gas_avg(iv) = 0.5_r8*(gas_avg(iv) + gas(iv))
 302   end do
 303
 304   call dumpxx( 'kk', dtchem, t_k, p_atm, ah2o, &
 305          jaerosolstate, jphase, jhyst_leg, &
 306          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 307          mosaic_vars_aa )
 308
 309   return
 310 end subroutine ASTEM
 311
 312
 313
 314 !-----------------------------------------------------------------------
 315   subroutine check_astem_negative( n, xnerr_astem_negative, &
 316                                       fix_astem_negative, aer, gas )
 317 !
 318 ! checks for negative values in gas and aer arrays
 319 ! when a negative value is found
 320 !    xnerr_astem_negative is incremented
 321 !    gas/aer value is set to 0.0
 322 !
 323   use module_data_mosaic_aero, only: naer, nbin_a, nbin_a_max, ngas_aerchtot
 324
 325   integer,  intent(in)    :: n, fix_astem_negative
 326   real(r8), intent(inout) :: xnerr_astem_negative(5,4)
 327   real(r8), intent(inout), dimension(ngas_aerchtot)     :: gas
 328   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 329
 330   character(len = 100) :: tmp_str
 331   integer  :: iaer, ibin, igas, j, m
 332   real(r8) :: tmpa
 333
 334   if ( n<1 .or. n>4 ) then
 335      write(tmp_str,'(/a,i10/)') '*** check_astem_negative fatal error, n =', n
 336      call mosaic_err_mess(tmp_str)
 337   end if
 338
 339   do igas = 1, ngas_aerchtot
 340      tmpa = gas(igas)
 341      if (tmpa >= 0.0_r8) then
 342         cycle
 343      else if (tmpa <= -1.0e-5_r8 ) then
 344         m = 1
 345      else if (tmpa <= -1.0e-10_r8) then
 346         m = 2
 347      else if (tmpa <= -1.0e-20_r8) then
 348         m = 3
 349      else if (tmpa <= -1.0e-30_r8) then
 350         m = 4
 351      else
 352         m = 5
 353      end if
 354      xnerr_astem_negative(m,n) = xnerr_astem_negative(m,n) + 1.0_r8
 355      if (fix_astem_negative > 0) gas(igas) = 0.0_r8
 356   end do
 357
 358   do ibin = 1, nbin_a
 359      do j = 1, 3
 360         do iaer = 1, naer
 361            tmpa = aer(iaer,j,ibin)
 362            if (tmpa >= 0.0_r8) then
 363               cycle
 364            else if (tmpa <= -1.0e-5_r8 ) then
 365               m = 1
 366            else if (tmpa <= -1.0e-10_r8) then
 367               m = 2
 368            else if (tmpa <= -1.0e-20_r8) then
 369               m = 3
 370            else if (tmpa <= -1.0e-30_r8) then
 371               m = 4
 372            else
 373               m = 5
 374            end if
 375            xnerr_astem_negative(m,n) = xnerr_astem_negative(m,n) + 1.0_r8
 376            if (fix_astem_negative > 0) aer(iaer,j,ibin) = 0.0_r8
 377         end do
 378      enddo
 379   enddo
 380
 381   end subroutine check_astem_negative
 382
 383
 384
 385 !***********************************************************************
 386 ! part of ASTEM: integrates semi-volatile inorganic gases
 387 !
 388 ! author: Rahul A. Zaveri
 389 ! update: feb 2015
 390 !-----------------------------------------------------------------------
 391 subroutine ASTEM_semi_volatiles( iprint_input,  dtchem, jaerosolstate,                &
 392      sfc_a, flux_s, flux_l, Heff, volatile_s, phi_volatile_s,     &
 393      jphase, aer, kg, gas, jhyst_leg, electrolyte, epercent, kel, activity, mc,            &
 394      delta_nh3_max, delta_hno3_max, delta_hcl_max,                                &
 395      num_a, mass_dry_a, mass_wet_a, mass_soluble_a,     &
 396      vol_dry_a, vol_wet_a, water_a, total_species, tot_cl_in,                         &
 397      aH2O_a, aH2O, ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, Keq_ll,  &
 398      Keq_gl, Keq_sg, Kp_nh4cl, Kp_nh4no3, Keq_nh4cl, MW_c, MW_a, mw_aer_mac,             &
 399      dens_aer_mac, Keq_sl, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa,           &
 400      jsalt_present, jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old,      &
 401      integrate, phi_volatile_l,                                                        &
 402      kappa_nonelectro, mosaic_vars_aa )
 403
 404   use module_data_mosaic_aero, only: nbin_a_max, nbin_a,                         &
 405        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
 406        Ncation, naer, mYES, mNO, ngas_ioa, jsolid, jliquid, all_solid, all_liquid,        &
 407        mixed, no_aerosol, jtotal, jhyst_lo, Nanion, nrxn_aer_gl, nrxn_aer_ll,           &
 408        nrxn_aer_sg, nrxn_aer_sl, nsalt, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM,    &
 409        jnh4cl, jnh4no3,                                                            &!TBD
 410        iso4_a, inh3_g, ihno3_g, ihcl_g,                                            & ! RAZ 2/2/2015: bugfix
 411        mosaic_vars_aa_type
 412
 413   use module_mosaic_ext, only: do_full_deliquescence,form_electrolytes
 414
 415
 416
 417   ! subr arguments
 418   integer, intent(inout) :: iprint_input
 419   integer, intent(in), dimension(nsalt) :: jsalt_index
 420   integer, intent(inout), dimension(nsalt) :: jsalt_present
 421   integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase
 422   integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
 423   integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
 424   integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich
 425   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
 426
 427   real(r8), intent(in)  :: dtchem
 428   real(r8), intent(in) :: aH2O,rtol_mesa
 429   real(r8), intent(inout) :: Kp_nh4cl,Kp_nh4no3,Keq_nh4cl
 430   real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
 431   real(r8), intent(in), dimension(Ncation) :: zc,MW_c
 432   real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
 433   real(r8), intent(in), dimension(Nanion)  :: za,MW_a
 434   real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
 435   real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
 436   real(r8), intent(inout), dimension(nbin_a_max) :: delta_hcl_max,water_a,num_a
 437   real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a,mass_wet_a,MDRH
 438   real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
 439   real(r8), intent(inout), dimension(nbin_a_max) :: aH2O_a,vol_wet_a,gam_ratio,growth_factor
 440   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
 441   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, total_species
 442   real(r8), intent(inout) :: tot_cl_in
 443   real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
 444   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
 445   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
 446   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
 447   real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
 448   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
 449   real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
 450   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,flux_l
 451   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: Heff,volatile_s
 452   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
 453   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
 454   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
 455   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
 456   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
 457   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
 458   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
 459   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 460   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 461   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
 462   real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
 463   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
 464
 465   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
 466
 467
 468
 469   ! local variables
 470   character(len=500) :: tmp_str
 471   integer ibin, iv, jp, ieqblm_ASTEM, islow_intermassxfer               ! RAZ 2/2/2015: bugfix
 472   integer, dimension(nbin_a_max) :: idry_case3a
 473
 474   real(r8) :: dtmax, t_new, t_old, t_out, XT,kelvin_nh4no3
 475   real(r8) :: sum1, sum2, sum3, sum4, sum4a, sum4b, h_flux_s
 476   real(r8) :: phi_nh4no3_s, phi_nh4cl_s,kelvin_nh4cl,Keq_nh4no3
 477   real(r8) :: sumkg_nh3,sumkg_hno3,sumkg_hcl                            ! RAZ 2/2/2015: bugfix
 478   real(r8), dimension(nbin_a_max) :: kgfrac_nh3,kgfrac_hno3,kgfrac_hcl  ! RAZ 2/2/2015: bugfix
 479   real(r8), dimension(ngas_volatile) :: sum_phi_volatile_s, sum_phi_volatile_l, sum_phi_volatile
 480   real(r8), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,df_gas_l
 481   real(r8), dimension(ngas_volatile,nbin_a_max) :: h_s_i_m
 482   real(r8), dimension(3,nbin_a_max) :: electrolyte_sum
 483
 484
 485   ! initialize time
 486   t_old = 0.0
 487   t_out = dtchem
 488
 489   ! reset ASTEM time steps and MESA iterations counters to zero
 490   mosaic_vars_aa%isteps_ASTEM = 0
 491   do ibin = 1, nbin_a
 492      mosaic_vars_aa%iter_MESA(ibin) = 0
 493   enddo
 494
 495 ! RAZ 2/2/2015: begin bugfix
 496   sumkg_nh3   = 0.0
 497   sumkg_hno3  = 0.0
 498   sumkg_hcl   = 0.0
 499   do ibin = 1, nbin_a
 500      sumkg_nh3   = sumkg_nh3   + kg(inh3_g,ibin)
 501      sumkg_hno3  = sumkg_hno3  + kg(ihno3_g,ibin)
 502      sumkg_hcl   = sumkg_hcl   + kg(ihcl_g,ibin)
 503   enddo
 504   do ibin = 1, nbin_a
 505      kgfrac_nh3(ibin)  = kg(inh3_g,ibin)/sumkg_nh3
 506      kgfrac_hno3(ibin) = kg(ihno3_g,ibin)/sumkg_hno3
 507      kgfrac_hcl(ibin)  = kg(ihcl_g,ibin)/sumkg_hcl
 508   enddo
 509 ! RAZ 2/2/2015: end bugfix
 510
 511
 512   !--------------------------------
 513   ! overall integration loop begins over dtchem seconds
 514
 515 10 mosaic_vars_aa%isteps_ASTEM = mosaic_vars_aa%isteps_ASTEM + 1
 516
 517   ! compute new fluxes
 518   phi_nh4no3_s = 0.0
 519   phi_nh4cl_s  = 0.0
 520   ieqblm_ASTEM = mYES                   ! reset to default
 521
 522   do 501 ibin = 1, nbin_a
 523
 524      idry_case3a(ibin) = mNO                    ! reset to default
 525      ! default fluxes and other stuff
 526      do iv = 1, ngas_ioa
 527         sfc_a(iv)                  = gas(iv)
 528         df_gas_s(iv,ibin)          = 0.0
 529         df_gas_l(iv,ibin)          = 0.0
 530         flux_s(iv,ibin)            = 0.0
 531         flux_l(iv,ibin)            = 0.0
 532         Heff(iv,ibin)              = 0.0
 533         volatile_s(iv,ibin)        = 0.0
 534         phi_volatile_s(iv,ibin)    = 0.0
 535         phi_volatile_l(iv,ibin)    = 0.0
 536         integrate(iv,jsolid,ibin)  = mNO        ! reset to default
 537         integrate(iv,jliquid,ibin) = mNO        ! reset to default
 538      enddo
 539
 540 ! RAZ 2/2/2015: begin bugfix
 541 ! Added this block here to prevent aer going negative in "absorb_tiny_******" subroutines
 542 ! update estimated possible condensation for each bin - used to calculate "tiny" amounts
 543      if(jaerosolstate(ibin) .ne. no_aerosol)then
 544         delta_nh3_max(ibin) = 0.1*gas(inh3_g)*kgfrac_nh3(ibin)
 545         delta_hno3_max(ibin)= 0.1*gas(ihno3_g)*kgfrac_hno3(ibin)
 546         delta_hcl_max(ibin) = 0.1*gas(ihcl_g)*kgfrac_hcl(ibin)
 547      endif
 548 ! RAZ 2/2/2015: end bugfix
 549
 550
 551      if(jaerosolstate(ibin) .eq. all_solid)then
 552         jphase(ibin) = jsolid
 553         call ASTEM_flux_dry(ibin, phi_nh4no3_s, phi_nh4cl_s, ieqblm_ASTEM,           &
 554              idry_case3a, sfc_a, df_gas_s, flux_s, phi_volatile_s, integrate, aer, kg,   &
 555              gas, electrolyte, epercent, Keq_sg)
 556
 557      elseif(jaerosolstate(ibin) .eq. all_liquid)then
 558         jphase(ibin) = jliquid
 559         call ASTEM_flux_wet(ibin, ieqblm_ASTEM, sfc_a, df_gas_s, df_gas_l,            &
 560              jaerosolstate, flux_s, Heff, phi_volatile_s, phi_volatile_l, integrate,   &
 561              jphase, aer, kg, gas, jhyst_leg, electrolyte, kel, activity, mc,             &
 562              delta_nh3_max, delta_hno3_max, delta_hcl_max, Keq_nh4cl, Keq_nh4no3,     &
 563              num_a, electrolyte_sum, mass_dry_a, mass_soluble_a, water_a, aH2O,        &
 564              kelvin_nh4no3, kelvin_nh4cl, ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a,    &
 565              na_Ma, nc_Mc, xeq_c, mw_electrolyte, Kp_nh4cl, Kp_nh4no3, Keq_gl, Keq_ll,   &
 566              MW_c, MW_a, total_species, tot_cl_in, molality0,                          &
 567              kappa_nonelectro, mosaic_vars_aa                                      )
 568
 569      elseif(jaerosolstate(ibin) .eq. mixed)then
 570         call ASTEM_flux_mix(ibin, phi_nh4no3_s, phi_nh4cl_s, ieqblm_ASTEM,           &
 571              idry_case3a, sfc_a, df_gas_s, df_gas_l, jaerosolstate, flux_s, Heff,       &
 572              phi_volatile_s, phi_volatile_l, integrate, jphase, aer, kg, gas, jhyst_leg, &
 573              electrolyte, epercent, kel, activity, mc, delta_nh3_max, delta_hno3_max,   &
 574              delta_hcl_max, Keq_nh4cl, Keq_nh4no3, num_a, electrolyte_sum, mass_dry_a, &
 575              mass_soluble_a, water_a, aH2O, kelvin_nh4no3, kelvin_nh4cl, ma, gam,       &
 576              log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte,     &
 577              Kp_nh4cl, Kp_nh4no3, Keq_ll, Keq_gl, Keq_sg, MW_c, MW_a, total_species,     &
 578              tot_cl_in, molality0, kappa_nonelectro, mosaic_vars_aa )   ! jphase(ibin) will be determined in this subr.
 579
 580      endif
 581
 582 501 continue
 583
 584   if(ieqblm_ASTEM .eq. mYES)goto 30     ! all bins have reached eqblm, so quit.
 585
 586
 587 ! RAZ 2/2/2015: new algorithm begin
 588 ! check if extremely slow inter-particle mass transfer is occurring after 20 ASTEM steps
 589   islow_intermassxfer = mNO ! default
 590   if(mosaic_vars_aa%isteps_ASTEM .gt. 20)then
 591     islow_intermassxfer = mYes ! default
 592
 593     do iv = 2, 4 ! HNO3, HCl, NH3
 594       sum_phi_volatile_s(iv) = sum(abs(phi_volatile_s(iv,1:nbin_a)))
 595       sum_phi_volatile_l(iv) = sum(abs(phi_volatile_l(iv,1:nbin_a)))
 596       sum_phi_volatile(iv) = sum_phi_volatile_s(iv) + sum_phi_volatile_l(iv)
 597
 598       if(gas(iv) .gt. 0.01 .and. sum_phi_volatile(iv) .gt. 0.01)islow_intermassxfer = mNO
 599
 600     enddo
 601   endif
 602
 603   if(islow_intermassxfer .eq. mYES)goto 30 ! extremely slow interparticle massxfer, so quit.
 604 ! RAZ 2/2/2015: new algorithm end
 605
 606
 607   !-------------------------
 608   ! calculate maximum possible internal time-step
 609 11 call ASTEM_calculate_dtmax( dtchem,  dtmax, jaerosolstate, idry_case3a, df_gas_s,   &
 610         flux_s, volatile_s, phi_volatile_l, integrate, aer, kg, gas, electrolyte,        &
 611         h_s_i_m, mosaic_vars_aa )
 612   t_new = t_old + dtmax ! update time
 613   if(t_new .gt. t_out)then      ! check if the new time step is too large
 614      dtmax = t_out - t_old
 615      t_new = t_out*1.01
 616   endif
 617
 618
 619   !------------------------------------------
 620   ! do internal time-step (dtmax) integration
 621
 622   do 20 iv = 2, 4
 623
 624      sum1 = 0.0
 625      sum2 = 0.0
 626      sum3 = 0.0
 627      sum4 = 0.0
 628      sum4a= 0.0
 629      sum4b= 0.0
 630
 631      do 21 ibin = 1, nbin_a
 632         if(jaerosolstate(ibin) .eq. no_aerosol)goto 21
 633
 634         jp = jliquid
 635         sum1 = sum1 + aer(iv,jp,ibin)/   &
 636              (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin)*integrate(iv,jp,ibin))
 637
 638         sum2 = sum2 + kg(iv,ibin)*integrate(iv,jp,ibin)/   &
 639              (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin)*integrate(iv,jp,ibin))
 640
 641         jp = jsolid
 642         sum3 = sum3 + aer(iv,jp,ibin)
 643
 644         if(flux_s(iv,ibin) .gt. 0.)then
 645            h_flux_s = dtmax*flux_s(iv,ibin)
 646            sum4a = sum4a + h_flux_s
 647            aer(iv,jp,ibin) = aer(iv,jp,ibin) + h_flux_s
 648         elseif(flux_s(iv,ibin) .lt. 0.)then
 649            h_flux_s = min(h_s_i_m(iv,ibin),dtmax)*flux_s(iv,ibin)
 650            sum4b = sum4b + h_flux_s
 651            aer(iv,jp,ibin) = aer(iv,jp,ibin) + h_flux_s
 652            aer(iv,jp,ibin) = max(aer(iv,jp,ibin), 0.0d0)
 653         endif
 654
 655 21   continue
 656
 657      sum4 = sum4a + sum4b
 658
 659
 660      ! first update gas concentration
 661      gas(iv) = (total_species(iv) - (sum1 + sum3 + sum4) )/   &
 662           (1. + dtmax*sum2)
 663      gas(iv) = max(gas(iv), 0.0d0)
 664
 665      !        if(gas(iv) .lt. 0.)write(6,*) gas(iv)
 666
 667      ! now update aer concentration in the liquid phase
 668      do 22 ibin = 1, nbin_a
 669
 670         if(integrate(iv,jliquid,ibin) .eq. mYES)then
 671            aer(iv,jliquid,ibin) =   &
 672                 (aer(iv,jliquid,ibin) + dtmax*kg(iv,ibin)*gas(iv))/   &
 673                 (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin))
 674         endif
 675
 676 22   continue
 677
 678
 679 20 continue
 680   !------------------------------------------
 681   ! sub-step integration done
 682
 683
 684   !------------------------------------------
 685   ! now update aer(jtotal) and update internal phase equilibrium
 686   ! also do integration of species by mass balance if necessary
 687   !
 688   do 40 ibin = 1, nbin_a
 689      if(jaerosolstate(ibin) .eq. no_aerosol)goto 40
 690
 691      if(jphase(ibin) .eq. jsolid)then
 692         call form_electrolytes(jsolid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)  ! degas excess nh3 (if present)
 693      elseif(jphase(ibin) .eq. jliquid)then
 694         call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in) ! degas excess nh3 (if present)
 695      elseif(jphase(ibin) .eq. jtotal)then
 696         call form_electrolytes(jsolid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)  ! degas excess nh3 (if present)
 697         call form_electrolytes(jliquid,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in) ! degas excess nh3 (if present)
 698      endif
 699
 700      !========================
 701      ! now update jtotal
 702      do iv = 2, ngas_ioa
 703         aer(iv,jtotal,ibin)=aer(iv,jsolid,ibin)+aer(iv,jliquid,ibin)
 704      enddo
 705      !========================
 706
 707
 708      call form_electrolytes(jtotal,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in) ! for MDRH diagnosis
 709
 710
 711
 712      ! update internal phase equilibrium
 713      if(jhyst_leg(ibin) .eq. jhyst_lo)then
 714         call ASTEM_update_phase_eqblm(ibin, jaerosolstate,     &
 715              jphase, aer, jhyst_leg, electrolyte, epercent, activity, mc, num_a,         &
 716              mass_dry_a, mass_wet_a, mass_soluble_a, vol_dry_a, vol_wet_a, water_a,    &
 717              aH2O_a, aH2O, ma, gam, log_gamZ, zc, za,    &
 718              gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac,         &
 719              dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll, growth_factor, MDRH, MDRH_T,      &
 720              molality0, rtol_mesa, jsalt_present, jsalt_index, jsulf_poor, jsulf_rich, &
 721              phi_salt_old,                                                   &
 722              kappa_nonelectro, mosaic_vars_aa )
 723         if (mosaic_vars_aa%f_mos_fail > 0) then
 724            return
 725         endif
 726      else
 727         call do_full_deliquescence(ibin,aer,electrolyte)                ! simply do liquid <-- total
 728      endif
 729
 730
 731 40 continue
 732   !------------------------------------------
 733
 734   ! update time
 735   t_old = t_new
 736
 737   if(mosaic_vars_aa%isteps_ASTEM .ge. mosaic_vars_aa%nmax_ASTEM)then
 738      mosaic_vars_aa%jASTEM_fail = mosaic_vars_aa%jASTEM_fail + 1
 739      write(tmp_str,*)'ASTEM internal steps exceeded', mosaic_vars_aa%nmax_ASTEM
 740      call mosaic_warn_mess(trim(adjustl(tmp_str)))
 741
 742      write(tmp_str,*)'ibin =', ibin
 743      call mosaic_warn_mess(trim(adjustl(tmp_str)))
 744
 745      if(iprint_input .eq. mYES)then
 746         !          call print_input
 747         iprint_input = mNO
 748      endif
 749      goto 30
 750   elseif(t_new .lt. t_out)then
 751      goto 10
 752   endif
 753
 754
 755   ! check if end of dtchem reached
 756   if(t_new .lt. 0.9999*t_out) goto 10
 757
 758 30 mosaic_vars_aa%cumul_steps_ASTEM = mosaic_vars_aa%cumul_steps_ASTEM + mosaic_vars_aa%isteps_ASTEM
 759    mosaic_vars_aa%isteps_ASTEM_max = max( mosaic_vars_aa%isteps_ASTEM_max, mosaic_vars_aa%isteps_ASTEM )
 760   !================================================
 761   ! end of overall integration loop over dtchem seconds
 762
 763
 764   !
 765   ! call subs to calculate fluxes over mixed-phase particles to update H+ ions,
 766   ! which were wiped off during update_phase_eqblm
 767   do ibin = 1, nbin_a
 768
 769      if(jaerosolstate(ibin) .eq. mixed)then
 770         if( electrolyte(jnh4no3,jsolid,ibin).gt. 0.0 .or.   &
 771              electrolyte(jnh4cl, jsolid,ibin).gt. 0.0 )then
 772            call ASTEM_flux_mix(ibin, phi_nh4no3_s, phi_nh4cl_s, ieqblm_ASTEM,        &
 773                 idry_case3a, sfc_a, df_gas_s, df_gas_l, jaerosolstate, flux_s, Heff,    &
 774                 phi_volatile_s, phi_volatile_l, integrate, jphase, aer, kg, gas,        &
 775                 jhyst_leg, electrolyte, epercent, kel, activity, mc, delta_nh3_max,     &
 776                 delta_hno3_max, delta_hcl_max, Keq_nh4cl, Keq_nh4no3, num_a,          &
 777                 electrolyte_sum, mass_dry_a, mass_soluble_a, water_a, aH2O,           &
 778                 kelvin_nh4no3, kelvin_nh4cl, ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, &
 779                 na_Ma, nc_Mc, xeq_c, mw_electrolyte, Kp_nh4cl, Kp_nh4no3, Keq_ll,       &
 780                 Keq_gl, Keq_sg, MW_c, MW_a, total_species, tot_cl_in, molality0,        &
 781                 kappa_nonelectro, mosaic_vars_aa )              ! jphase(ibin) will be determined in this subr.
 782         else
 783            jphase(ibin) = jliquid
 784            call ASTEM_flux_wet(ibin, ieqblm_ASTEM, sfc_a, df_gas_s, df_gas_l,         &
 785                 jaerosolstate, flux_s, Heff, phi_volatile_s, phi_volatile_l,          &
 786                 integrate, jphase, aer, kg, gas, jhyst_leg, electrolyte, kel, activity,   &
 787                 mc, delta_nh3_max, delta_hno3_max, delta_hcl_max, Keq_nh4cl,          &
 788                 Keq_nh4no3, num_a, electrolyte_sum, mass_dry_a, mass_soluble_a,       &
 789                 water_a, aH2O, kelvin_nh4no3, kelvin_nh4cl, ma, gam, log_gamZ, zc, za,    &
 790                 gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, Kp_nh4cl,        &
 791                 Kp_nh4no3, Keq_gl, Keq_ll, MW_c, MW_a, total_species, tot_cl_in,        &
 792                 molality0, kappa_nonelectro, mosaic_vars_aa )
 793         endif
 794      endif
 795
 796   enddo
 797
 798
 799   return
 800 end subroutine ASTEM_semi_volatiles
 801
 802
 803
 804 !***********************************************************************
 805 ! part of ASTEM: computes max time step for gas-aerosol integration
 806 !
 807 ! author: Rahul A. Zaveri
 808 ! update: jan 2005
 809 !-----------------------------------------------------------------------
 810 subroutine ASTEM_calculate_dtmax( dtchem, dtmax, jaerosolstate, idry_case3a,       &
 811      df_gas_s, flux_s, volatile_s, phi_volatile_l, integrate, aer, kg, gas, electrolyte, &
 812      h_s_i_m, mosaic_vars_aa )
 813
 814   use module_data_mosaic_aero, only: r8, nbin_a_max,                       &
 815        ngas_aerchtot, ngas_volatile, nelectrolyte,                         &
 816        naer, ngas_ioa, mYES, jliquid, jsolid, no_aerosol,                  &
 817        nbin_a, alpha_astem,                                                &
 818        jnh4no3, ino3_a, jnh4cl, inh4_a, icl_a,                             &
 819        mosaic_vars_aa_type
 820
 821
 822
 823   ! subr arguments
 824   integer, intent(in), dimension(nbin_a_max) :: jaerosolstate,idry_case3a
 825   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
 826
 827   real(r8), intent(in)  :: dtchem
 828   real(r8), intent(out) :: dtmax
 829   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
 830   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s
 831   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,volatile_s
 832   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
 833   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
 834   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: h_s_i_m
 835   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 836   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 837
 838   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
 839
 840   ! local variables
 841   character(len=500) :: tmp_str
 842   integer ::  ibin, iv
 843   real(r8) :: alpha, h_gas, h_sub_max,h_gas_i(ngas_ioa), h_gas_l, h_gas_s
 844   real(r8) :: sum_kg_phi, sum_kg_phi_pos, sum_kg_phi_neg, sumflux_s     ! RAZ 2/2/2015: revised algorithm
 845   real(r8), dimension(ngas_volatile) :: sum_bin_s,sum_vdf_s,sum_vol_s
 846   real(r8), dimension(ngas_volatile) :: avg_df_gas_s
 847
 848   h_sub_max = dtchem/5.0        ! sec RAZ 2/14/2014
 849
 850   ! GAS-SIDE
 851
 852   ! solid-phase
 853   ! calculate h_gas_i and h_gas_l
 854
 855   h_gas_s = 2.e16
 856
 857   do 5 iv = 2, ngas_ioa
 858      h_gas_i(iv) = 1.e16
 859      sumflux_s = 0.0
 860      do ibin = 1, nbin_a
 861         if(flux_s(iv,ibin) .gt. 0.0)then
 862            sumflux_s = sumflux_s + flux_s(iv,ibin)
 863         endif
 864      enddo
 865
 866      if(sumflux_s .gt. 0.0)then
 867         h_gas_i(iv) = 0.1*gas(iv)/sumflux_s
 868         h_gas_s     = min(h_gas_s, h_gas_i(iv))
 869      endif
 870
 871 5 continue
 872
 873
 874   ! liquid-phase
 875   ! calculate h_gas_s and h_gas_l
 876
 877   h_gas_l = 2.e16
 878   do 6 iv = 2, ngas_ioa
 879      h_gas_i(iv) = 1.e16
 880      sum_kg_phi = 0.0
 881      sum_kg_phi_pos = 0.0
 882      sum_kg_phi_neg = 0.0
 883      do ibin = 1, nbin_a
 884         if(integrate(iv,jliquid,ibin) .eq. mYES)then
 885
 886 !           sum_kg_phi = sum_kg_phi +   &
 887 !                abs(phi_volatile_l(iv,ibin))*kg(iv,ibin)
 888
 889 ! RAZ 2/2/2015: revised algorithm: begin
 890           if(phi_volatile_l(iv,ibin) .gt. 0.0)then
 891            sum_kg_phi_pos = sum_kg_phi_pos + abs(phi_volatile_l(iv,ibin))*kg(iv,ibin)
 892           else
 893            sum_kg_phi_neg = sum_kg_phi_neg + abs(phi_volatile_l(iv,ibin))*kg(iv,ibin)
 894           endif
 895 ! RAZ 2/2/2015: revised algorithm: end
 896
 897         endif
 898      enddo
 899
 900      sum_kg_phi = max(sum_kg_phi_pos, sum_kg_phi_neg) ! RAZ 2/2/2015: revised algorithm
 901
 902      if(sum_kg_phi .gt. 0.0)then
 903         h_gas_i(iv) = alpha_astem/sum_kg_phi
 904         h_gas_l     = min(h_gas_l, h_gas_i(iv))
 905      endif
 906
 907 6 continue
 908
 909   h_gas = min(h_gas_s, h_gas_l)
 910   h_gas = min(h_gas, h_sub_max)
 911
 912
 913
 914
 915   ! AEROSOL-SIDE: solid-phase
 916
 917   ! first load volatile_solid array
 918   do ibin = 1, nbin_a
 919
 920      volatile_s(ino3_a,ibin) = electrolyte(jnh4no3,jsolid,ibin)
 921      volatile_s(inh4_a,ibin) = electrolyte(jnh4cl,jsolid,ibin) +   &
 922           electrolyte(jnh4no3,jsolid,ibin)
 923
 924      if(idry_case3a(ibin) .eq. mYES)then
 925         volatile_s(icl_a,ibin)  = aer(icl_a,jsolid,ibin)
 926      else
 927         volatile_s(icl_a,ibin)  = electrolyte(jnh4cl,jsolid,ibin)
 928      endif
 929
 930   enddo
 931
 932
 933   ! next calculate weighted avg_df_gas_s
 934   do iv = 2, ngas_ioa
 935
 936      sum_bin_s(iv) = 0.0
 937      sum_vdf_s(iv) = 0.0
 938      sum_vol_s(iv) = 0.0
 939
 940      do ibin = 1, nbin_a
 941         if(flux_s(iv,ibin) .lt. 0.)then ! aer -> gas
 942            sum_bin_s(iv) = sum_bin_s(iv) + 1.0
 943            sum_vdf_s(iv) = sum_vdf_s(iv) +   &
 944                 volatile_s(iv,ibin)*df_gas_s(iv,ibin)
 945            sum_vol_s(iv) = sum_vol_s(iv) + volatile_s(iv,ibin)
 946         endif
 947      enddo
 948
 949      if(sum_vol_s(iv) .gt. 0.0)then
 950         avg_df_gas_s(iv) = sum_vdf_s(iv)/sum_vol_s(iv)
 951      else
 952         avg_df_gas_s(iv) = 1.0 ! never used, but set to 1.0 just to be safe
 953      endif
 954
 955   enddo
 956
 957
 958   ! calculate h_s_i_m
 959
 960
 961   do 20 ibin = 1, nbin_a
 962
 963      if(jaerosolstate(ibin) .eq. no_aerosol) goto 20
 964
 965      do 10 iv = 2, ngas_ioa
 966
 967         if(flux_s(iv,ibin) .lt. 0.)then                         ! aer -> gas
 968
 969            alpha = abs(avg_df_gas_s(iv))/   &
 970                 (volatile_s(iv,ibin)*sum_bin_s(iv))
 971            alpha = min(alpha, 1.0d0)
 972
 973            if(idry_case3a(ibin) .eq. mYES)alpha = 1.0
 974
 975            h_s_i_m(iv,ibin) =   &
 976                 -alpha*volatile_s(iv,ibin)/flux_s(iv,ibin)
 977
 978         endif
 979
 980 10   continue
 981
 982
 983 20 continue
 984
 985
 986   dtmax = min(dtchem, h_gas)
 987
 988 !  dtmax = h_sub_max
 989
 990
 991   if(dtmax .eq. 0.0)then
 992      write(tmp_str,*)' dtmax = ', dtmax
 993      call mosaic_warn_mess(trim(adjustl(tmp_str)))
 994   endif
 995
 996   return
 997 end subroutine ASTEM_calculate_dtmax
 998
 999
1000
1001 !***********************************************************************
1002 ! part of ASTEM: updates solid-liquid partitioning after each gas-aerosol
1003 ! mass transfer step
1004 !
1005 ! author: Rahul A. Zaveri
1006 ! update: sep 2015
1007 !
1008 ! 9/3/2015 RAZ: Bugfix - fixed phase state calculations for aerosols that dont contain any salts,
1009 !               but can still contain water due to presence of BC, OC, SOA, and OIN, which are now
1010 !               allowed to absorb some water.
1011 !-----------------------------------------------------------------------
1012 subroutine ASTEM_update_phase_eqblm(ibin, jaerosolstate,       &
1013      jphase, aer, jhyst_leg, electrolyte, epercent, activity, mc, num_a, mass_dry_a,    &
1014      mass_wet_a, mass_soluble_a, vol_dry_a, vol_wet_a, water_a, aH2O_a, aH2O,           &
1015      ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc,                         &
1016      xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_sl, MW_c, MW_a, Keq_ll,       &
1017      growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index,     &
1018      jsulf_poor, jsulf_rich, phi_salt_old,                                              &
1019      kappa_nonelectro, mosaic_vars_aa )
1020
1021   use module_data_mosaic_aero,  only: r8, nbin_a_max, nelectrolyte, Ncation, naer,      &
1022        jtotal, nsalt, all_solid, jsolid, all_liquid, jliquid, jhyst_lo, jhyst_up,       &
1023        Nanion, nrxn_aer_ll, nrxn_aer_sl, MDRH_T_NUM,  &
1024        jsulf_poor_NUM, jsulf_rich_NUM,                                                  &
1025        ptol_mol_astem, mhyst_force_lo,  mhyst_force_up,                                 &
1026        jcacl2, jcano3, mhyst_method,                                                    &
1027        mosaic_vars_aa_type
1028
1029   use module_mosaic_ext,  only: do_full_deliquescence, adjust_solid_aerosol,            &
1030        MESA_PTC, calculate_XT, aerosol_water, adjust_liquid_aerosol,                    &
1031        compute_activities
1032
1033
1034   ! subr arguments
1035   integer, intent(in):: ibin
1036   integer, intent(in), dimension(nsalt) :: jsalt_index
1037   integer, intent(inout), dimension(nsalt) :: jsalt_present
1038   integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
1039   integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
1040   integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich
1041
1042   real(r8), intent(in) :: aH2O,rtol_mesa
1043   real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
1044   real(r8), intent(in), dimension(Ncation) :: zc,MW_c
1045   real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
1046   real(r8), intent(in), dimension(Nanion)  :: za,MW_a
1047   real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
1048   real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,mass_wet_a
1049   real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,vol_dry_a
1050   real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,water_a,gam_ratio
1051   real(r8), intent(inout), dimension(nbin_a_max) :: aH2O_a,growth_factor,MDRH
1052   real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
1053   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1054   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1055   real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
1056   real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
1057   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
1058   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1059   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1060   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1061   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1062   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1063   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
1064   real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
1065   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
1066
1067   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
1068
1069   ! local variables
1070   integer jsalt_dum, js, j_index, je
1071   real(r8) :: CRH, XT, sum_dum
1072
1073
1074   !! EFFI calculate percent composition
1075   sum_dum = 0.0
1076   do je = 1, nelectrolyte
1077      sum_dum = sum_dum + electrolyte(je,jtotal,ibin)
1078   enddo
1079
1080   if(sum_dum .eq. 0.)sum_dum = 1.0
1081
1082   do je = 1, nelectrolyte
1083      epercent(je,jtotal,ibin) = 100.*electrolyte(je,jtotal,ibin)/sum_dum
1084   enddo
1085   !! EFFI
1086
1087
1088   ! calculate overall sulfate ratio
1089   call calculate_XT(ibin,jtotal,XT,aer)         ! calc updated XT
1090
1091
1092 !! begin new algorithm - 6/3/2015 RAZ
1093     jsalt_dum = 0       ! 9/3/2015 RAZ
1094     do js = 1, nsalt
1095        jsalt_present(js) = 0                            ! default value - salt absent
1096
1097        if(epercent(js,jtotal,ibin) .gt. ptol_mol_astem)then
1098           jsalt_present(js) = 1                         ! salt present
1099           jsalt_dum = jsalt_dum + jsalt_index(js)       ! 9/3/2015 RAZ
1100        endif
1101     enddo
1102
1103
1104     if( (epercent(jcano3,jtotal,ibin) .gt. ptol_mol_astem) .or. &
1105         (epercent(jcacl2,jtotal,ibin) .gt. ptol_mol_astem) )then
1106       CRH = 0.0  ! no crystrallization or efflorescence point
1107     else
1108       CRH = 0.35 ! default value
1109     endif
1110
1111
1112     ! now diagnose MDRH
1113     if(jsalt_dum .eq. 0)then                    ! no salts or acids are present ! 9/3/2015 RAZ: updated algorithm for jsalt_dum = 0
1114
1115        CRH = 0.0
1116        MDRH(ibin) = 0.0
1117        jaerosolstate(ibin) = all_solid
1118        jphase(ibin)    = jsolid
1119        jhyst_leg(ibin) = jhyst_lo
1120        call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
1121        water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &    ! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
1122         electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
1123        return
1124
1125     elseif(XT .lt. 1. .and. XT .gt. 0.0)then    ! excess sulfate, always liquid, MDRH=0.0
1126        MDRH(ibin) = 0.0
1127     elseif(XT .ge. 2.0 .or. XT .lt. 0.0)then    ! sulfate poor
1128        j_index = jsulf_poor(jsalt_dum)          ! 9/3/2015 RAZ
1129        MDRH(ibin) = MDRH_T(j_index)
1130     else                                        ! sulfate rich
1131        j_index = jsulf_rich(jsalt_dum)          ! 9/3/2015 RAZ
1132        MDRH(ibin) = MDRH_T(j_index)
1133     endif
1134
1135     CRH = min(CRH, MDRH(ibin)/100.0)            ! 6/3/2015 RAZ
1136
1137 !! end new algorithm - 6/3/2015 RAZ
1138
1139
1140     ! modified step 1: 9/3/2015 RAZ
1141     ! step 1: check if aH2O is below CRH (crystallization or efflorescence point)
1142     if( aH2O_a(ibin).lt.CRH )then
1143        jaerosolstate(ibin) = all_solid
1144        jphase(ibin)    = jsolid
1145        jhyst_leg(ibin) = jhyst_lo
1146        call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
1147        water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &    ! 9/3/2015 RAZ: water due to nonelectrolytes (OC, BC, SOA, OIN)
1148         electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
1149        return
1150     endif
1151
1152     ! step 2: check mhyst_method
1153     if(mhyst_method == mhyst_force_up .or. jhyst_leg(ibin) == jhyst_up) then ! 9/3/2015 RAZ: either forced up OR (new) already fully deliquesced (may be metastable), so continue on upper leg
1154        call do_full_deliquescence(ibin,aer,electrolyte) ! this call is probably not necessary, but do it just to be safe
1155        jaerosolstate(ibin) = all_liquid
1156        jhyst_leg(ibin) = jhyst_up
1157        jphase(ibin) = jliquid
1158        water_a(ibin) = aerosol_water(jtotal,ibin,jaerosolstate,jphase,jhyst_leg,   &
1159           electrolyte,aer,kappa_nonelectro,num_a,mass_dry_a,mass_soluble_a,aH2O,molality0)
1160
1161        if(water_a(ibin) .le. 0.0)then     ! one last attempt to catch bad input
1162           jaerosolstate(ibin) = all_solid ! no soluble material present
1163           jphase(ibin)    = jsolid
1164           jhyst_leg(ibin) = jhyst_lo
1165           call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
1166        else
1167           call adjust_liquid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent)
1168           call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,      &
1169                   electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a, &
1170                   aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)
1171        endif
1172
1173        return
1174     endif
1175
1176
1177     ! step 3: diagnose phase state based on MDRH
1178     if(aH2O*100. .lt. MDRH(ibin)) then
1179        jaerosolstate(ibin) = all_solid
1180        jphase(ibin) = jsolid
1181        call adjust_solid_aerosol(ibin,jphase,aer,jhyst_leg,electrolyte,epercent,water_a)
1182        return
1183     endif
1184
1185
1186     ! step 4: none of the above means it must be sub-saturated or mixed-phase
1187 10  if(jphase(ibin) .eq. jsolid)then
1188       call do_full_deliquescence(ibin,aer,electrolyte)
1189       call MESA_PTC( ibin, jaerosolstate, jphase, aer,          &
1190           jhyst_leg, electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a, &
1191           mass_soluble_a, vol_dry_a, vol_wet_a, water_a, aH2O,          &
1192           ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma,     &
1193           nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_sl, MW_c, MW_a,    &
1194           Keq_ll, growth_factor, molality0, rtol_mesa, jsalt_present,                  &
1195           phi_salt_old, kappa_nonelectro, mosaic_vars_aa )
1196     else
1197       call MESA_PTC( ibin, jaerosolstate, jphase, aer,          &
1198           jhyst_leg, electrolyte, epercent, activity, mc, num_a, mass_dry_a, mass_wet_a, &
1199           mass_soluble_a, vol_dry_a, vol_wet_a, water_a, aH2O,          &
1200           ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma,     &
1201           nc_Mc, xeq_c, mw_electrolyte, mw_aer_mac, dens_aer_mac, Keq_sl, MW_c, MW_a,    &
1202           Keq_ll, growth_factor, molality0, rtol_mesa, jsalt_present,                  &
1203           phi_salt_old, kappa_nonelectro, mosaic_vars_aa )
1204     endif
1205     return
1206
1207   end subroutine ASTEM_update_phase_eqblm
1208
1209
1210
1211 !==================================================================
1212 !
1213 ! LIQUID PARTICLES
1214 !
1215 !***********************************************************************
1216 ! part of ASTEM: computes fluxes over wet aerosols
1217 !
1218 ! author: Rahul A. Zaveri
1219 ! update: Jan 2007
1220 !-----------------------------------------------------------------------
1221 subroutine ASTEM_flux_wet(ibin, ieqblm_ASTEM, sfc_a, df_gas_s, df_gas_l,              &
1222      jaerosolstate, flux_s, Heff, phi_volatile_s, phi_volatile_l, integrate, jphase,    &
1223      aer, kg, gas, jhyst_leg, electrolyte, kel, activity, mc, delta_nh3_max,              &
1224      delta_hno3_max, delta_hcl_max, Keq_nh4cl, Keq_nh4no3, num_a, electrolyte_sum,     &
1225      mass_dry_a, mass_soluble_a, water_a, aH2O, kelvin_nh4no3, kelvin_nh4cl, ma, gam,    &
1226      log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c, mw_electrolyte, Kp_nh4cl,    &
1227      Kp_nh4no3, Keq_gl, Keq_ll, MW_c, MW_a, total_species, tot_cl_in, molality0,         &
1228      kappa_nonelectro, mosaic_vars_aa                                            )
1229
1230   use module_data_mosaic_aero, only: r8, nbin_a_max,                             &
1231        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
1232        Ncation, naer, jliquid, jsolid, mNO, mYES, Nanion, nrxn_aer_gl, nrxn_aer_ll,   &
1233        jcaco3, inh4_a, inh3_g, ihno3_g, ino3_a, ihcl_g, icl_a, jnh4no3, jnh4cl,       &
1234        mosaic_vars_aa_type
1235
1236   use module_mosaic_ext,  only: compute_activities, ions_to_electrolytes,           &
1237        absorb_tiny_nh4no3, absorb_tiny_nh4cl, absorb_tiny_hno3, absorb_tiny_hcl
1238
1239
1240   ! subr arguments
1241   integer, intent(in) :: ibin
1242   integer, intent(inout) :: ieqblm_ASTEM
1243   integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
1244   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1245
1246   real(r8), intent(in) :: aH2O
1247   real(r8), intent(inout) :: Keq_nh4cl,Keq_nh4no3,kelvin_nh4no3,kelvin_nh4cl
1248   real(r8), intent(inout) :: Kp_nh4cl,Kp_nh4no3
1249   real(r8), intent(in), dimension(Ncation) :: zc,MW_c
1250   real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
1251   real(r8), intent(in), dimension(Nanion)  :: za,MW_a
1252   real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
1253   real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
1254   real(r8), intent(inout), dimension(nbin_a_max) :: delta_hcl_max
1255   real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_dry_a,gam_ratio
1256   real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a,water_a
1257   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1258   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, total_species
1259   real(r8), intent(inout) :: tot_cl_in
1260   real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
1261   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1262   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1263   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1264   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s
1265   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l
1266   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,Heff
1267   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
1268   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
1269   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1270   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
1271   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1272   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1273   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1274   real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
1275   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1276   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1277   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
1278
1279   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
1280
1281   ! local variables
1282   character(len=500) :: tmp_str
1283   integer iv, iadjust, iadjust_intermed
1284   real(r8) :: XT, g_nh3_hno3, g_nh3_hcl, a_nh4_no3, a_nh4_cl
1285
1286   call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,   &
1287        nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)    ! for water content calculation
1288   call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg, &
1289        electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,    &
1290        water_a,aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)
1291
1292   if(water_a(ibin) .eq. 0.0)then
1293      write(tmp_str,*)'Water is zero in liquid phase'
1294      call mosaic_warn_mess(trim(adjustl(tmp_str)))
1295      write(tmp_str,*)'Stopping in ASTEM_flux_wet'
1296      call mosaic_warn_mess(trim(adjustl(tmp_str)))
1297      mosaic_vars_aa%zero_water_flag = .true.
1298   endif
1299
1300   !-------------------------------------------------------------------
1301   ! CASE 1: caco3 > 0 absorb acids (and indirectly degas co2)
1302
1303   if(electrolyte(jcaco3,jsolid,ibin) .gt. 0.0)then
1304      call ASTEM_flux_wet_case1(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,          &
1305           phi_volatile_s,integrate,jphase,kg,gas,mc,Keq_ll)
1306      return
1307   endif
1308
1309   !-------------------------------------------------------------------
1310   ! CASE 2: Sulfate-Rich Domain
1311
1312 ! if(XT.lt.1.9999 .and. XT.ge.0.)then  ! RAZ 11/10/2014
1313   if(XT.lt.2.0 .and. XT.ge.0.)then  ! RAZ 11/10/2014
1314      call ASTEM_flux_wet_case2(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,            &
1315           phi_volatile_l,integrate,gas,kel,mc,water_a,ma,gam,gam_ratio,Keq_ll,   &
1316           Keq_gl)
1317      return
1318   endif
1319
1320   !-------------------------------------------------------------------
1321
1322   if( (gas(inh3_g)+aer(inh4_a,jliquid,ibin)) .lt. 1.e-25)goto 10  ! no ammonia in the system
1323
1324   !-------------------------------------------------------------------
1325   ! CASE 3: nh4no3 and/or nh4cl maybe active
1326   ! do some small adjustments (if needed) before deciding case 3
1327
1328   iadjust = mNO         ! default
1329   iadjust_intermed = mNO        ! default
1330
1331   ! nh4no3
1332   g_nh3_hno3 = gas(inh3_g)*gas(ihno3_g)
1333   a_nh4_no3  = aer(inh4_a,jliquid,ibin)*aer(ino3_a,jliquid,ibin)
1334
1335   if(g_nh3_hno3 .gt. 0. .and. a_nh4_no3 .eq. 0.)then
1336      call absorb_tiny_nh4no3(ibin,aer,gas,electrolyte,delta_nh3_max,             &
1337           delta_hno3_max,electrolyte_sum)
1338      iadjust = mYES
1339      iadjust_intermed = mYES
1340   endif
1341
1342   if(iadjust_intermed .eq. mYES)then
1343      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
1344           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
1345      iadjust_intermed = mNO     ! reset
1346   endif
1347
1348   ! nh4cl
1349   g_nh3_hcl = gas(inh3_g)*gas(ihcl_g)
1350   a_nh4_cl  = aer(inh4_a,jliquid,ibin)*aer(icl_a,jliquid,ibin)
1351
1352   if(g_nh3_hcl .gt. 0. .and. a_nh4_cl .eq. 0.)then
1353      call absorb_tiny_nh4cl(ibin,aer,gas,electrolyte,delta_nh3_max,delta_hcl_max,&
1354           electrolyte_sum)
1355      iadjust = mYES
1356      iadjust_intermed = mYES
1357   endif
1358
1359   if(iadjust_intermed .eq. mYES)then
1360      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
1361           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
1362   endif
1363
1364   if(iadjust .eq. mYES)then
1365      call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte,&
1366           activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,       &
1367           log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                 ! update after adjustments
1368   endif
1369
1370
1371   ! all adjustments done...
1372
1373   !--------
1374   kelvin_nh4no3 = kel(inh3_g,ibin)*kel(ihno3_g,ibin)
1375   Keq_nh4no3 = kelvin_nh4no3*activity(jnh4no3,ibin)*Kp_nh4no3   ! = [NH3]s * [HNO3]s
1376
1377   kelvin_nh4cl = kel(inh3_g,ibin)*kel(ihcl_g,ibin)
1378   Keq_nh4cl = kelvin_nh4cl*activity(jnh4cl,ibin)*Kp_nh4cl       ! = [NH3]s * [HCl]s
1379
1380   call ASTEM_flux_wet_case3(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,phi_volatile_l,&
1381        integrate,kg,gas,kel,mc,Keq_nh4cl,Keq_nh4no3,water_a,ma,gam,gam_ratio,    &
1382        Keq_ll,Keq_gl,aer,total_species,tot_cl_in,activity,electrolyte)
1383
1384   return
1385
1386
1387   !-------------------------------------------------------------------
1388   ! CASE 4: ammonia = 0. hno3 and hcl exchange may happen here
1389   ! do small adjustments (if needed) before deciding case 4
1390
1391 10 iadjust = mNO                ! default
1392   iadjust_intermed = mNO        ! default
1393
1394   ! hno3
1395   if(gas(ihno3_g).gt.0. .and. aer(ino3_a,jliquid,ibin).eq.0. .and.   &
1396        aer(icl_a,jliquid,ibin) .gt. 0.0)then
1397      call absorb_tiny_hno3(ibin,aer,gas,delta_hno3_max) ! and degas tiny hcl
1398      iadjust = mYES
1399      iadjust_intermed = mYES
1400   endif
1401
1402   if(iadjust_intermed .eq. mYES)then
1403      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
1404           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
1405      iadjust_intermed = mNO     ! reset
1406   endif
1407
1408   ! hcl
1409   if(gas(ihcl_g).gt.0. .and. aer(icl_a,jliquid,ibin) .eq. 0. .and.   &
1410        aer(ino3_a,jliquid,ibin) .gt. 0.0)then
1411      call absorb_tiny_hcl(ibin,aer,gas,delta_hcl_max)   ! and degas tiny hno3
1412      iadjust = mYES
1413      iadjust_intermed = mYES
1414   endif
1415
1416   if(iadjust_intermed .eq. mYES)then
1417      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
1418           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
1419   endif
1420
1421   if(iadjust .eq. mYES)then
1422      call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte,&
1423           activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,       &
1424           log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                 ! update after adjustments
1425   endif
1426
1427   ! all adjustments done...
1428
1429   call ASTEM_flux_wet_case4(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,phi_volatile_l,&
1430        integrate,kg,gas,kel,mc,water_a,ma,gam,Keq_ll,Keq_gl)
1431
1432
1433   return
1434 end subroutine ASTEM_flux_wet
1435
1436
1437
1438 !***********************************************************************
1439 ! part of ASTEM: subroutines for flux_wet cases
1440 !
1441 ! author: Rahul A. Zaveri
1442 ! update: Jan 2007
1443 !-----------------------------------------------------------------------
1444
1445 ! CASE 1: CaCO3 > 0 absorb all acids (and indirectly degas co2)
1446
1447 subroutine ASTEM_flux_wet_case1(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,         &
1448      phi_volatile_s,integrate,jphase,kg,gas,mc,Keq_ll)
1449
1450   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1451        Ncation,mYES, jsolid,mNO,nrxn_aer_ll,                                       &
1452        jc_h,ihno3_g,ihcl_g
1453
1454
1455   ! subr arguments
1456   integer, intent(in):: ibin
1457   integer, intent(inout) :: ieqblm_ASTEM
1458   integer, intent(inout), dimension(nbin_a_max) :: jphase
1459   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1460
1461   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1462   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
1463   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1464   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
1465   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
1466   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
1467   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1468
1469   ! local variables
1470   integer iv
1471
1472   mc(jc_h,ibin) = sqrt(Keq_ll(3))
1473
1474   ! same as dry case1
1475   if(gas(ihno3_g) .gt. 1.e-6)then
1476      sfc_a(ihno3_g) = 0.0
1477      df_gas_s(ihno3_g,ibin) = gas(ihno3_g)
1478      phi_volatile_s(ihno3_g,ibin) = 1.0
1479      flux_s(ihno3_g,ibin) = kg(ihno3_g,ibin)*df_gas_s(ihno3_g,ibin)
1480      integrate(ihno3_g,jsolid,ibin) = mYES
1481      jphase(ibin) = jsolid
1482      ieqblm_ASTEM = mNO
1483   endif
1484
1485   if(gas(ihcl_g) .gt. 1.e-6)then
1486      sfc_a(ihcl_g)  = 0.0
1487      df_gas_s(ihcl_g,ibin) = gas(ihcl_g)
1488      phi_volatile_s(ihcl_g,ibin) = 1.0
1489      flux_s(ihcl_g,ibin) = kg(ihcl_g,ibin)*df_gas_s(ihcl_g,ibin)
1490      integrate(ihcl_g,jsolid,ibin)  = mYES
1491      jphase(ibin) = jsolid
1492      ieqblm_ASTEM = mNO
1493   endif
1494
1495   return
1496 end subroutine ASTEM_flux_wet_case1
1497
1498
1499
1500 !--------------------------------------------------------------------
1501 ! CASE 2: Sulfate-Rich Domain
1502
1503 subroutine ASTEM_flux_wet_case2(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,           &
1504      phi_volatile_l,integrate,gas,kel,mc,water_a,ma,gam,gam_ratio,Keq_ll,Keq_gl)
1505
1506   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1507        Ncation,mYES, jliquid,mNO,Nanion,nelectrolyte,nrxn_aer_gl,nrxn_aer_ll,      &
1508        jc_h,jc_nh4,inh3_g,jhno3,ja_no3,ihno3_g,jhcl,ja_cl,ihcl_g
1509
1510
1511
1512   ! subr arguments
1513   integer, intent(in) :: ibin
1514   integer, intent(inout) :: ieqblm_ASTEM
1515   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1516
1517   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
1518   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1519   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
1520   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1521   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1522   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l,Heff
1523   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1524   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kel
1525   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: gam
1526   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1527   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1528   ! local variables
1529   real(r8) :: dum_hno3, dum_hcl, dum_nh3
1530
1531
1532   sfc_a(inh3_g)  = kel(inh3_g,ibin)*   &
1533        gam_ratio(ibin)*mc(jc_nh4,ibin)*Keq_ll(3)/   &
1534        (mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1535
1536   sfc_a(ihno3_g) = kel(ihno3_g,ibin)*   &
1537        mc(jc_h,ibin)*ma(ja_no3,ibin)*gam(jhno3,ibin)**2/   &
1538        Keq_gl(3)
1539
1540   sfc_a(ihcl_g)  = kel(ihcl_g,ibin)*   &
1541        mc(jc_h,ibin)*ma(ja_cl,ibin)*gam(jhcl,ibin)**2/   &
1542        Keq_gl(4)
1543
1544   dum_hno3 = max(sfc_a(ihno3_g), gas(ihno3_g))
1545   dum_hcl  = max(sfc_a(ihcl_g), gas(ihcl_g))
1546   dum_nh3  = max(sfc_a(inh3_g), gas(inh3_g))
1547
1548
1549   ! compute relative driving forces
1550   if(dum_hno3 .gt. 0.0)then
1551      df_gas_l(ihno3_g,ibin) = gas(ihno3_g) - sfc_a(ihno3_g)
1552      phi_volatile_l(ihno3_g,ibin)= df_gas_l(ihno3_g,ibin)/dum_hno3
1553   else
1554      phi_volatile_l(ihno3_g,ibin)= 0.0
1555   endif
1556
1557   if(dum_hcl .gt. 0.0)then
1558      df_gas_l(ihcl_g,ibin)  = gas(ihcl_g)  - sfc_a(ihcl_g)
1559      phi_volatile_l(ihcl_g,ibin) = df_gas_l(ihcl_g,ibin)/dum_hcl
1560   else
1561      phi_volatile_l(ihcl_g,ibin) = 0.0
1562   endif
1563
1564   if(dum_nh3 .gt. 0.0)then
1565      df_gas_l(inh3_g,ibin)  = gas(inh3_g)  - sfc_a(inh3_g)
1566      phi_volatile_l(inh3_g,ibin) = df_gas_l(inh3_g,ibin)/dum_nh3
1567   else
1568      phi_volatile_l(inh3_g,ibin) = 0.0
1569   endif
1570
1571
1572   !      if(phi_volatile_l(ihno3_g,ibin) .le. rtol_eqb_astem .and.
1573   !     &   phi_volatile_l(ihcl_g,ibin)  .le. rtol_eqb_astem .and.
1574   !     &   phi_volatile_l(inh3_g,ibin)  .le. rtol_eqb_astem)then
1575   !
1576   !        return
1577   !
1578   !      endif
1579
1580
1581   ! compute Heff
1582   if(dum_hno3 .gt. 0.0)then
1583      Heff(ihno3_g,ibin)=   &
1584           kel(ihno3_g,ibin)*gam(jhno3,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1585           (water_a(ibin)*Keq_gl(3))
1586      integrate(ihno3_g,jliquid,ibin)= mYES
1587      ieqblm_ASTEM = mNO
1588   endif
1589
1590   if(dum_hcl .gt. 0.0)then
1591      Heff(ihcl_g,ibin)=   &
1592           kel(ihcl_g,ibin)*gam(jhcl,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1593           (water_a(ibin)*Keq_gl(4))
1594      integrate(ihcl_g,jliquid,ibin) = mYES
1595      ieqblm_ASTEM = mNO
1596   endif
1597
1598   if(dum_nh3 .gt. 0.0)then
1599      Heff(inh3_g,ibin) =   &
1600           kel(inh3_g,ibin)*gam_ratio(ibin)*1.e-9*Keq_ll(3)/   &
1601           (water_a(ibin)*mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1602      integrate(inh3_g,jliquid,ibin) = mYES
1603      ieqblm_ASTEM = mNO
1604   endif
1605
1606
1607   return
1608 end subroutine ASTEM_flux_wet_case2
1609
1610
1611
1612 !---------------------------------------------------------------------
1613 ! CASE 3: nh4no3 and/or nh4cl may be active
1614
1615 subroutine ASTEM_flux_wet_case3(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,           &
1616      phi_volatile_l,integrate,kg,gas,kel,mc,Keq_nh4cl,Keq_nh4no3,water_a,ma,gam, &
1617      gam_ratio,Keq_ll,Keq_gl,aer,total_species,tot_cl_in,activity,electrolyte)
1618
1619   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1620        Ncation,mYES, jliquid,mNO,Nanion,nelectrolyte,nrxn_aer_gl,nrxn_aer_ll,naer, &
1621        inh3_g,ihcl_g,ihno3_g,ja_no3,jhno3,jc_h,ja_cl,jhcl,jc_nh4
1622
1623   use module_mosaic_ext, only: quadratic,equilibrate_acids
1624
1625
1626   ! subr arguments
1627   integer, intent(in) :: ibin
1628   integer, intent(inout) :: ieqblm_ASTEM
1629   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1630
1631   real(r8), intent(inout) :: Keq_nh4cl,Keq_nh4no3
1632   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
1633   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1634   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, total_species
1635   real(r8), intent(inout) :: tot_cl_in
1636   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1637   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1638   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l,Heff
1639   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1640   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
1641   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: gam,activity
1642   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1643   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1644   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1645   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1646   ! local variables
1647   real(r8) :: a, b, c, dum_hno3, dum_hcl, dum_nh3
1648   ! function
1649   !real(r8) :: quadratic
1650
1651   a =   kg(inh3_g,ibin)
1652   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
1653        + kg(ihno3_g,ibin)*gas(ihno3_g)   &
1654        + kg(ihcl_g,ibin)*gas(ihcl_g)
1655   c = -(kg(ihno3_g,ibin)*Keq_nh4no3 + kg(ihcl_g,ibin)*Keq_nh4cl)
1656
1657   sfc_a(inh3_g)  = quadratic(a,b,c)
1658   sfc_a(ihno3_g) = Keq_nh4no3/max(sfc_a(inh3_g),1.d-20)
1659   sfc_a(ihcl_g)  = Keq_nh4cl/max(sfc_a(inh3_g),1.d-20)
1660
1661
1662   ! diagnose mH+
1663   if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
1664      mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
1665           (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
1666   elseif(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
1667      mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
1668           (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
1669   else
1670      call equilibrate_acids(ibin,aer,gas,electrolyte,activity,mc,water_a,       &
1671        total_species,tot_cl_in,ma,gam,Keq_ll,Keq_gl)    ! hno3 and/or hcl may be > 0 in the gas phase
1672      mc(jc_h,ibin)  = max(mc(jc_h,ibin), sqrt(Keq_ll(3)))
1673
1674      sfc_a(inh3_g)  = kel(inh3_g,ibin)*   &
1675           gam_ratio(ibin)*mc(jc_nh4,ibin)*Keq_ll(3)/   &
1676           (mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1677
1678      sfc_a(ihno3_g) = kel(ihno3_g,ibin)*   &
1679           mc(jc_h,ibin)*ma(ja_no3,ibin)*gam(jhno3,ibin)**2/   &
1680           Keq_gl(3)
1681      sfc_a(ihcl_g)  = kel(ihcl_g,ibin)*   &
1682           mc(jc_h,ibin)*ma(ja_cl,ibin)*gam(jhcl,ibin)**2/   &
1683           Keq_gl(4)
1684   endif
1685
1686   dum_hno3 = max(sfc_a(ihno3_g), gas(ihno3_g))
1687   dum_hcl  = max(sfc_a(ihcl_g), gas(ihcl_g))
1688   dum_nh3  = max(sfc_a(inh3_g), gas(inh3_g))
1689
1690   ! compute relative driving forces
1691   if(dum_hno3 .gt. 0.0)then
1692      df_gas_l(ihno3_g,ibin) = gas(ihno3_g) - sfc_a(ihno3_g)
1693      phi_volatile_l(ihno3_g,ibin)= df_gas_l(ihno3_g,ibin)/dum_hno3
1694   else
1695      phi_volatile_l(ihno3_g,ibin)= 0.0
1696   endif
1697
1698   if(dum_hcl .gt. 0.0)then
1699      df_gas_l(ihcl_g,ibin)  = gas(ihcl_g)  - sfc_a(ihcl_g)
1700      phi_volatile_l(ihcl_g,ibin) = df_gas_l(ihcl_g,ibin)/dum_hcl
1701   else
1702      phi_volatile_l(ihcl_g,ibin) = 0.0
1703   endif
1704
1705   if(dum_nh3 .gt. 0.0)then
1706      df_gas_l(inh3_g,ibin)  = gas(inh3_g)  - sfc_a(inh3_g)
1707      phi_volatile_l(inh3_g,ibin) = df_gas_l(inh3_g,ibin)/dum_nh3
1708   else
1709      phi_volatile_l(inh3_g,ibin) = 0.0
1710   endif
1711
1712
1713
1714   !      if(phi_volatile_l(ihno3_g,ibin) .le. rtol_eqb_astem .and.
1715   !     &   phi_volatile_l(ihcl_g,ibin)  .le. rtol_eqb_astem .and.
1716   !     &   phi_volatile_l(inh3_g,ibin)  .le. rtol_eqb_astem)then
1717   !
1718   !        return
1719   !
1720   !      endif
1721
1722
1723   ! compute Heff
1724   if(dum_hno3 .gt. 0.0)then
1725      Heff(ihno3_g,ibin)=   &
1726           kel(ihno3_g,ibin)*gam(jhno3,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1727           (water_a(ibin)*Keq_gl(3))
1728      integrate(ihno3_g,jliquid,ibin)= mYES
1729      ieqblm_ASTEM = mNO
1730   endif
1731
1732   if(dum_hcl .gt. 0.0)then
1733      Heff(ihcl_g,ibin)=   &
1734           kel(ihcl_g,ibin)*gam(jhcl,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1735           (water_a(ibin)*Keq_gl(4))
1736      integrate(ihcl_g,jliquid,ibin) = mYES
1737      ieqblm_ASTEM = mNO
1738   endif
1739
1740   if(dum_nh3 .gt. 0.0)then
1741      Heff(inh3_g,ibin) =   &
1742           kel(inh3_g,ibin)*gam_ratio(ibin)*1.e-9*Keq_ll(3)/   &
1743           (water_a(ibin)*mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1744      integrate(inh3_g,jliquid,ibin) = mYES
1745      ieqblm_ASTEM = mNO
1746   endif
1747
1748   return
1749 end subroutine ASTEM_flux_wet_case3
1750
1751
1752
1753 !--------------------------------------------------------------------
1754 ! CASE 3a: only NH4NO3 (aq) active
1755
1756 subroutine ASTEM_flux_wet_case3a(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,          &
1757      phi_volatile_l,integrate,kg,gas,kel,mc,Keq_nh4no3,water_a,ma,gam,gam_ratio, &
1758      Keq_ll,Keq_gl)    ! NH4NO3 (aq)
1759
1760   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1761        Ncation,mYES, jliquid,mNO,Nanion,nelectrolyte,nrxn_aer_gl,nrxn_aer_ll,      &
1762        inh3_g,ihno3_g,ja_no3,jhno3,jc_h
1763
1764   use module_mosaic_ext, only: quadratic
1765
1766
1767   ! subr arguments
1768   integer, intent(in) :: ibin
1769   integer, intent(inout) :: ieqblm_ASTEM
1770   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1771
1772   real(r8), intent(inout) :: Keq_nh4no3
1773   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
1774   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1775   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
1776   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1777   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1778   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l,Heff
1779   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1780   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
1781   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: gam
1782   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1783   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1784   ! local variables
1785   real(r8) :: a, b, c, dum_hno3, dum_nh3
1786   ! function
1787   !real(r8) :: quadratic
1788
1789
1790   a =   kg(inh3_g,ibin)
1791   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
1792        + kg(ihno3_g,ibin)*gas(ihno3_g)
1793   c = -(kg(ihno3_g,ibin)*Keq_nh4no3)
1794
1795   sfc_a(inh3_g)  = quadratic(a,b,c)
1796   sfc_a(ihno3_g) = Keq_nh4no3/sfc_a(inh3_g)
1797
1798
1799   ! diagnose mH+
1800   if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
1801      mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
1802           (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
1803   else
1804      mc(jc_h,ibin) = sqrt(Keq_ll(3))
1805   endif
1806
1807
1808   ! compute Heff
1809   dum_hno3 = max(sfc_a(ihno3_g), gas(ihno3_g))
1810   dum_nh3  = max(sfc_a(inh3_g), gas(inh3_g))
1811
1812   ! compute relative driving forces
1813   if(dum_hno3 .gt. 0.0)then
1814      df_gas_l(ihno3_g,ibin) = gas(ihno3_g) - sfc_a(ihno3_g)
1815      phi_volatile_l(ihno3_g,ibin)= df_gas_l(ihno3_g,ibin)/dum_hno3
1816   else
1817      phi_volatile_l(ihno3_g,ibin)= 0.0
1818   endif
1819
1820   if(dum_nh3 .gt. 0.0)then
1821      df_gas_l(inh3_g,ibin)  = gas(inh3_g)  - sfc_a(inh3_g)
1822      phi_volatile_l(inh3_g,ibin) = df_gas_l(inh3_g,ibin)/dum_nh3
1823   else
1824      phi_volatile_l(inh3_g,ibin) = 0.0
1825   endif
1826
1827
1828   !      if(phi_volatile_l(ihno3_g,ibin) .le. rtol_eqb_astem .and.
1829   !     &   phi_volatile_l(inh3_g,ibin)  .le. rtol_eqb_astem)then
1830   !
1831   !        return
1832   !
1833   !      endif
1834
1835
1836   ! compute Heff
1837   Heff(ihno3_g,ibin)=   &
1838        kel(ihno3_g,ibin)*gam(jhno3,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1839        (water_a(ibin)*Keq_gl(3))
1840   integrate(ihno3_g,jliquid,ibin)= mYES
1841
1842
1843   Heff(inh3_g,ibin) =   &
1844        kel(inh3_g,ibin)*gam_ratio(ibin)*1.e-9*Keq_ll(3)/   &
1845        (water_a(ibin)*mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1846   integrate(inh3_g,jliquid,ibin) = mYES
1847
1848
1849   ieqblm_ASTEM = mNO
1850
1851
1852   return
1853 end subroutine ASTEM_flux_wet_case3a
1854
1855
1856
1857 !--------------------------------------------------------------------
1858 ! CASE 3b: only NH4Cl (aq) active
1859
1860 subroutine ASTEM_flux_wet_case3b(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,          &
1861      phi_volatile_l,integrate,kg,gas,kel,mc,Keq_nh4cl,water_a,ma,gam,gam_ratio,  &
1862      Keq_ll,Keq_gl)     ! NH4Cl (aq)
1863
1864   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1865        Ncation,mYES, jliquid,mNO,Nanion,nelectrolyte,nrxn_aer_gl,nrxn_aer_ll,      &
1866        inh3_g,ihcl_g,ja_cl,jhcl,jc_h
1867
1868   use module_mosaic_ext, only: quadratic
1869
1870
1871   ! subr arguments
1872   integer,intent(in) :: ibin
1873   integer, intent(inout) :: ieqblm_ASTEM
1874   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1875
1876   real(r8), intent(inout) :: Keq_nh4cl
1877   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
1878   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1879   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
1880   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1881   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1882   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l,Heff
1883   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1884   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
1885   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: gam
1886   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
1887   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1888   ! local variables
1889   real(r8) :: a, b, c, dum_hcl, dum_nh3
1890   ! function
1891   !real(r8) :: quadratic
1892
1893
1894   a =   kg(inh3_g,ibin)
1895   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
1896        + kg(ihcl_g,ibin)*gas(ihcl_g)
1897   c = -(kg(ihcl_g,ibin)*Keq_nh4cl)
1898
1899   sfc_a(inh3_g)  = quadratic(a,b,c)
1900   sfc_a(ihcl_g)  = Keq_nh4cl /sfc_a(inh3_g)
1901
1902
1903   ! diagnose mH+
1904   if(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
1905      mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
1906           (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
1907   else
1908      mc(jc_h,ibin) = sqrt(Keq_ll(3))
1909   endif
1910
1911
1912   ! compute Heff
1913   dum_hcl  = max(sfc_a(ihcl_g), gas(ihcl_g))
1914   dum_nh3  = max(sfc_a(inh3_g), gas(inh3_g))
1915
1916
1917   ! compute relative driving forces
1918   if(dum_hcl .gt. 0.0)then
1919      df_gas_l(ihcl_g,ibin)  = gas(ihcl_g)  - sfc_a(ihcl_g)
1920      phi_volatile_l(ihcl_g,ibin) = df_gas_l(ihcl_g,ibin)/dum_hcl
1921   else
1922      phi_volatile_l(ihcl_g,ibin) = 0.0
1923   endif
1924
1925   if(dum_nh3 .gt. 0.0)then
1926      df_gas_l(inh3_g,ibin)  = gas(inh3_g)  - sfc_a(inh3_g)
1927      phi_volatile_l(inh3_g,ibin) = df_gas_l(inh3_g,ibin)/dum_nh3
1928   else
1929      phi_volatile_l(inh3_g,ibin) = 0.0
1930   endif
1931
1932
1933
1934   !      if(phi_volatile_l(ihcl_g,ibin)  .le. rtol_eqb_astem .and.
1935   !     &   phi_volatile_l(inh3_g,ibin)  .le. rtol_eqb_astem)then
1936   !
1937   !        return
1938   !
1939   !      endif
1940
1941
1942
1943   ! compute Heff
1944   Heff(ihcl_g,ibin)=   &
1945        kel(ihcl_g,ibin)*gam(jhcl,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
1946        (water_a(ibin)*Keq_gl(4))
1947   integrate(ihcl_g,jliquid,ibin) = mYES
1948
1949
1950   Heff(inh3_g,ibin) =   &
1951        kel(inh3_g,ibin)*gam_ratio(ibin)*1.e-9*Keq_ll(3)/   &
1952        (water_a(ibin)*mc(jc_h,ibin)*Keq_ll(2)*Keq_gl(2))
1953   integrate(inh3_g,jliquid,ibin) = mYES
1954
1955
1956   ieqblm_ASTEM = mNO
1957
1958
1959
1960   return
1961 end subroutine ASTEM_flux_wet_case3b
1962
1963
1964
1965 !-----------------------------------------------------------------------
1966 ! CASE 4: NH3 = 0 (in gas and aerosol). hno3 and hcl exchange may happen here
1967
1968 subroutine ASTEM_flux_wet_case4(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,           &
1969      phi_volatile_l,integrate,kg,gas,kel,mc,water_a,ma,gam,Keq_ll,Keq_gl)
1970
1971   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
1972        Ncation,mYES, jliquid,mNO,Nanion,nelectrolyte,nrxn_aer_gl,nrxn_aer_ll,      &
1973        jhno3,ja_no3,ihno3_g,jhcl,ja_cl,ihcl_g,jc_h
1974
1975
1976   ! subr arguments
1977   integer, intent(in) :: ibin
1978   integer, intent(inout) :: ieqblm_ASTEM
1979   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
1980
1981   real(r8), intent(inout), dimension(nbin_a_max) :: water_a
1982   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1983   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
1984   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
1985   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
1986   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l,Heff
1987   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
1988   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
1989   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max):: gam
1990   real(r8), intent(inout), dimension(Ncation,nbin_a_max) ::mc
1991   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
1992   ! local variables
1993   real(r8) :: dum_numer, dum_denom, gas_eqb_ratio, dum_hno3, dum_hcl
1994
1995
1996   dum_numer = kel(ihno3_g,ibin)*Keq_gl(4)*ma(ja_no3,ibin)*   &
1997        gam(jhno3,ibin)**2
1998   dum_denom = kel(ihcl_g,ibin)*Keq_gl(3)*ma(ja_cl ,ibin)*   &
1999        gam(jhcl,ibin)**2
2000
2001
2002   if(dum_denom .eq. 0.0 .or. dum_numer .eq. 0.0)then
2003      mc(jc_h,ibin) = sqrt(Keq_ll(3))
2004      return
2005   endif
2006
2007   gas_eqb_ratio = dum_numer/dum_denom   ! Ce,hno3/Ce,hcl
2008
2009
2010   ! compute equilibrium surface concentrations
2011   sfc_a(ihcl_g) =   &
2012        ( kg(ihno3_g,ibin)*gas(ihno3_g) + kg(ihcl_g,ibin)*gas(ihcl_g) )/   &
2013        ( kg(ihcl_g,ibin) + gas_eqb_ratio*kg(ihno3_g,ibin) )
2014   sfc_a(ihno3_g)= gas_eqb_ratio*sfc_a(ihcl_g)
2015
2016
2017   ! diagnose mH+
2018   if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
2019      mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
2020           (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
2021   elseif(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
2022      mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
2023           (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
2024   else
2025      mc(jc_h,ibin) = sqrt(Keq_ll(3))
2026   endif
2027
2028
2029   ! compute Heff
2030   dum_hno3 = max(sfc_a(ihno3_g), gas(ihno3_g))
2031   dum_hcl  = max(sfc_a(ihcl_g), gas(ihcl_g))
2032
2033   ! compute relative driving forces
2034   if(dum_hno3 .gt. 0.0)then
2035      df_gas_l(ihno3_g,ibin) = gas(ihno3_g) - sfc_a(ihno3_g)
2036      phi_volatile_l(ihno3_g,ibin)= df_gas_l(ihno3_g,ibin)/dum_hno3
2037   else
2038      phi_volatile_l(ihno3_g,ibin)= 0.0
2039   endif
2040
2041   if(dum_hcl .gt. 0.0)then
2042      df_gas_l(ihcl_g,ibin)  = gas(ihcl_g)  - sfc_a(ihcl_g)
2043      phi_volatile_l(ihcl_g,ibin)= df_gas_l(ihcl_g,ibin)/dum_hcl
2044   else
2045      phi_volatile_l(ihcl_g,ibin)= 0.0
2046   endif
2047
2048
2049   !      if(phi_volatile_l(ihno3_g,ibin) .le. rtol_eqb_astem .and.
2050   !     &   phi_volatile_l(ihcl_g,ibin)  .le. rtol_eqb_astem)then
2051   !
2052   !        return
2053   !
2054   !      endif
2055
2056
2057
2058   ! compute Heff
2059   Heff(ihno3_g,ibin)=   &
2060        kel(ihno3_g,ibin)*gam(jhno3,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
2061        (water_a(ibin)*Keq_gl(3))
2062   integrate(ihno3_g,jliquid,ibin)= mYES
2063
2064
2065   Heff(ihcl_g,ibin)=   &
2066        kel(ihcl_g,ibin)*gam(jhcl,ibin)**2*mc(jc_h,ibin)*1.e-9/   &
2067        (water_a(ibin)*Keq_gl(4))
2068   integrate(ihcl_g,jliquid,ibin) = mYES
2069
2070
2071   ieqblm_ASTEM = mNO
2072
2073
2074
2075   return
2076 end subroutine ASTEM_flux_wet_case4
2077
2078
2079
2080 !===========================================================
2081 !
2082 ! DRY PARTICLES
2083 !
2084 !===========================================================
2085 !***********************************************************************
2086 ! part of ASTEM: computes gas-aerosol fluxes over dry aerosols
2087 !
2088 ! author: Rahul A. Zaveri
2089 ! update: dec 2006
2090 !-----------------------------------------------------------------------
2091 subroutine ASTEM_flux_dry(ibin, phi_nh4no3_s,phi_nh4cl_s,ieqblm_ASTEM,           &
2092      idry_case3a,sfc_a,df_gas_s,flux_s,phi_volatile_s,integrate,aer,kg,gas,      &
2093      electrolyte,epercent,Keq_sg)
2094
2095   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2096        naer,jsolid, nrxn_aer_sg,nelectrolyte,                                      &
2097        jcaco3,jcacl2,jnacl,ihno3_g,jnh4cl,ihcl_g,inh3_g,jnh4no3
2098
2099   use module_mosaic_ext, only: calculate_XT
2100
2101
2102   ! subr arguments
2103   integer, intent(in) :: ibin
2104   integer, intent(inout) :: ieqblm_ASTEM
2105   integer, intent(inout), dimension(nbin_a_max) :: idry_case3a
2106   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2107
2108   real(r8), intent(out) :: phi_nh4no3_s,phi_nh4cl_s
2109   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2110   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2111   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2112   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s
2113   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s
2114   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2115   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2116   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2117   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2118   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
2119   ! local variables
2120   integer iv
2121   real(r8) :: XT, prod_nh4no3, prod_nh4cl, volatile_cl
2122
2123
2124
2125   call calculate_XT(ibin,jsolid,XT,aer)
2126
2127   !-----------------------------------------------------------------
2128   ! CASE 1:  caco3 > 0 absorb all acids (and indirectly degas co2)
2129
2130   if(electrolyte(jcaco3,jsolid,ibin) .gt. 0.0)then
2131      call ASTEM_flux_dry_case1(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,          &
2132           phi_volatile_s,integrate,kg,gas)
2133
2134      return
2135   endif
2136
2137   !-----------------------------------------------------------------
2138   ! CASE 2: Sulfate-Rich Domain
2139
2140 ! if(XT.lt.1.9999 .and. XT.ge.0.)then   ! excess sulfate (acidic) ! RAZ 11/10/2014
2141   if(XT.lt.2.0 .and. XT.ge.0.)then   ! excess sulfate (acidic) ! RAZ 11/10/2014
2142      call ASTEM_flux_dry_case2(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,          &
2143           phi_volatile_s,integrate,kg,gas)
2144
2145      return
2146   endif
2147
2148   !-------------------------------------------------------------------
2149   ! CASE 3: hno3 and hcl exchange may happen here and nh4cl may form/evaporate
2150
2151   volatile_cl  = electrolyte(jnacl,jsolid,ibin) +   &
2152        electrolyte(jcacl2,jsolid,ibin)
2153
2154
2155   if(volatile_cl .gt. 0.0 .and. gas(ihno3_g).gt. 0.0 )then
2156
2157      call ASTEM_flux_dry_case3a(ibin,ieqblm_ASTEM,idry_case3a,sfc_a,df_gas_s,    &
2158           flux_s,phi_volatile_s,integrate,aer,kg,gas)
2159
2160      prod_nh4cl = max( (gas(inh3_g)*gas(ihcl_g)-Keq_sg(2)), 0.0d0) +   &
2161           electrolyte(jnh4cl, jsolid,ibin)
2162
2163      if(prod_nh4cl .gt. 0.0)then
2164         call ASTEM_flux_dry_case3b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s, &
2165              flux_s,phi_volatile_s,integrate,aer,kg,gas,electrolyte,epercent,    &
2166              Keq_sg)
2167      endif
2168
2169      return
2170   endif
2171
2172   !-----------------------------------------------------------------
2173   ! CASE 4: nh4no3 or nh4cl or both may be active
2174
2175   prod_nh4no3 = max( (gas(inh3_g)*gas(ihno3_g)-Keq_sg(1)), 0.0d0) +   &
2176        electrolyte(jnh4no3,jsolid,ibin)
2177   prod_nh4cl  = max( (gas(inh3_g)*gas(ihcl_g) -Keq_sg(2)), 0.0d0) +   &
2178        electrolyte(jnh4cl, jsolid,ibin)
2179
2180   if(prod_nh4no3 .gt. 0.0 .or. prod_nh4cl .gt. 0.0)then
2181      call ASTEM_flux_dry_case4(ibin,phi_nh4no3_s,phi_nh4cl_s,ieqblm_ASTEM,sfc_a, &
2182           df_gas_s,flux_s,phi_volatile_s,integrate,kg,gas,electrolyte,epercent,  &
2183           Keq_sg,aer)
2184      return
2185   endif
2186
2187   !-----------------------------------------------------------------
2188
2189   return
2190 end subroutine ASTEM_flux_dry
2191
2192
2193
2194 !----------------------------------------------------------------------
2195
2196 !***********************************************************************
2197 ! part of ASTEM: subroutines for flux_dry cases
2198 !
2199 ! author: Rahul A. Zaveri
2200 ! update: dec 2006
2201 !-----------------------------------------------------------------------
2202
2203 ! CASE 1:  caco3 > 0 absorb all acids (and indirectly degas co2)
2204
2205 subroutine ASTEM_flux_dry_case1(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,         &
2206      phi_volatile_s,integrate,kg,gas)
2207
2208   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2209        mYES,jsolid,mNO, ihno3_g,ihcl_g
2210
2211
2212   ! subr arguments
2213   integer, intent(in) :: ibin
2214   integer, intent(inout) :: ieqblm_ASTEM
2215   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2216
2217   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2218   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2219   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2220   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2221   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2222
2223
2224   if(gas(ihno3_g) .gt. 1.e-6)then
2225      sfc_a(ihno3_g) = 0.0
2226      df_gas_s(ihno3_g,ibin) = gas(ihno3_g)
2227      phi_volatile_s(ihno3_g,ibin) = 1.0
2228      flux_s(ihno3_g,ibin) = kg(ihno3_g,ibin)*df_gas_s(ihno3_g,ibin)
2229      integrate(ihno3_g,jsolid,ibin) = mYES
2230      ieqblm_ASTEM = mNO
2231   endif
2232
2233   if(gas(ihcl_g) .gt. 1.e-6)then
2234      sfc_a(ihcl_g)  = 0.0
2235      df_gas_s(ihcl_g,ibin) = gas(ihcl_g)
2236      phi_volatile_s(ihcl_g,ibin) = 1.0
2237      flux_s(ihcl_g,ibin)  = kg(ihcl_g,ibin)*df_gas_s(ihcl_g,ibin)
2238      integrate(ihcl_g,jsolid,ibin)  = mYES
2239      ieqblm_ASTEM = mNO
2240   endif
2241
2242
2243   return
2244 end subroutine ASTEM_flux_dry_case1
2245
2246
2247
2248 !---------------------------------------------------------------------
2249 ! CASE 2: Sulfate-Rich Domain
2250
2251 subroutine ASTEM_flux_dry_case2(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,         &
2252      phi_volatile_s,integrate,kg,gas) ! TOUCH
2253
2254   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2255        mYES,jsolid,mNO, inh3_g
2256
2257
2258
2259   ! subr arguments
2260   integer, intent(in) :: ibin
2261   integer, intent(inout) :: ieqblm_ASTEM
2262   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2263
2264   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2265   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2266   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2267   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2268   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2269
2270
2271   if(gas(inh3_g).gt.1.e-6)then
2272      sfc_a(inh3_g) = 0.0
2273      df_gas_s(inh3_g,ibin) = gas(inh3_g)
2274      phi_volatile_s(inh3_g,ibin)  = 1.0
2275      flux_s(inh3_g,ibin) = kg(inh3_g,ibin)*gas(inh3_g)
2276      integrate(inh3_g,jsolid,ibin) = mYES
2277      ieqblm_ASTEM = mNO
2278   endif
2279
2280
2281   return
2282 end subroutine ASTEM_flux_dry_case2
2283
2284
2285
2286 !---------------------------------------------------------------------
2287 ! CASE 3a: degas hcl from nacl or cacl2 by flux_s balance with hno3
2288
2289 subroutine ASTEM_flux_dry_case3a(ibin,ieqblm_ASTEM,idry_case3a,sfc_a,df_gas_s,   &
2290      flux_s,phi_volatile_s,integrate,aer,kg,gas)
2291
2292   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2293        naer,jsolid,mYES,mNO,                                                       &
2294        ihno3_g,icl_a,ihcl_g
2295
2296
2297   ! subr arguments
2298   integer, intent(in) :: ibin
2299   integer, intent(inout) :: ieqblm_ASTEM
2300   integer, intent(inout), dimension(nbin_a_max) :: idry_case3a  ! changed "out" to "inout" RAZ 11/11/2014
2301   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2302
2303   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2304   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2305   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2306   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2307   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2308   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2309
2310
2311   if(gas(ihno3_g) .gt. 1.e-6)then
2312      sfc_a(ihno3_g) = 0.0
2313      sfc_a(ihcl_g)  = gas(ihcl_g) + aer(icl_a,jsolid,ibin)
2314
2315      df_gas_s(ihno3_g,ibin) = gas(ihno3_g)
2316      df_gas_s(ihcl_g,ibin)  = -aer(icl_a,jsolid,ibin)
2317
2318      flux_s(ihno3_g,ibin) = kg(ihno3_g,ibin)*gas(ihno3_g)
2319      flux_s(ihcl_g,ibin)  = -flux_s(ihno3_g,ibin)
2320
2321      phi_volatile_s(ihno3_g,ibin) = 1.0
2322      phi_volatile_s(ihcl_g,ibin)=df_gas_s(ihcl_g,ibin)/sfc_a(ihcl_g)
2323
2324      integrate(ihno3_g,jsolid,ibin) = mYES
2325      integrate(ihcl_g,jsolid,ibin)  = mYES
2326
2327      idry_case3a(ibin) = mYES
2328      ieqblm_ASTEM = mNO
2329   endif
2330
2331   return
2332 end subroutine ASTEM_flux_dry_case3a
2333
2334
2335
2336
2337 !---------------------------------------------------------------------
2338 ! CASE 3b: nh4cl may form/evaporate here
2339
2340 subroutine ASTEM_flux_dry_case3b(ibin, phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,  &
2341      flux_s,phi_volatile_s,integrate,aer,kg,gas,electrolyte,epercent,Keq_sg)      ! TOUCH
2342
2343   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2344        naer,nsalt, jsolid,mYES,mNO,nrxn_aer_sg,                                    &
2345        rtol_eqb_ASTEM,ptol_mol_ASTEM,                                              &
2346        nelectrolyte,jnh4cl,ihcl_g,inh3_g,icl_a
2347
2348   use module_mosaic_ext, only: quadratic,degas_solid_nh4cl
2349
2350
2351   ! subr arguments
2352   integer, intent(in) :: ibin
2353   integer, intent(inout) ::ieqblm_ASTEM
2354   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2355
2356   real(r8), intent(out) :: phi_nh4cl_s
2357   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2358   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2359   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2360   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2361   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2362   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2363   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2364   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2365   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
2366   ! local variables
2367   integer iactive_nh4cl, js
2368   real(r8) :: a, b, c
2369   real(r8) :: sum_dum
2370   ! function
2371   !real(r8) :: quadratic
2372
2373
2374
2375   !! EFFI calculate percent composition
2376   sum_dum = 0.0
2377   do js = 1, nsalt
2378      sum_dum = sum_dum + electrolyte(js,jsolid,ibin)
2379   enddo
2380
2381   if(sum_dum .eq. 0.)sum_dum = 1.0
2382
2383   epercent(jnh4cl,jsolid,ibin) = 100.*electrolyte(jnh4cl,jsolid,ibin)/sum_dum
2384   !! EFFI
2385
2386
2387
2388   !-------------------
2389   ! set default values for flags
2390   iactive_nh4cl  = 1
2391
2392
2393   ! compute relative driving force
2394   phi_nh4cl_s = (gas(inh3_g)*gas(ihcl_g) - Keq_sg(2))/   &
2395        max(gas(inh3_g)*gas(ihcl_g),Keq_sg(2))
2396
2397
2398   !-------------------
2399   ! now determine if nh4cl is active or significant
2400   ! nh4cl
2401   if( abs(phi_nh4cl_s) .lt. rtol_eqb_ASTEM )then
2402      iactive_nh4cl = 0
2403   elseif(gas(inh3_g)*gas(ihcl_g) .lt. Keq_sg(2) .and.   &
2404        epercent(jnh4cl, jsolid,ibin) .le. ptol_mol_ASTEM)then
2405      iactive_nh4cl = 0
2406      if(epercent(jnh4cl, jsolid,ibin) .gt. 0.0)then
2407         call degas_solid_nh4cl(ibin,aer,gas,electrolyte,Keq_sg)
2408      endif
2409   endif
2410
2411
2412   ! check the outcome
2413   if(iactive_nh4cl .eq. 0)return
2414
2415
2416   !-----------------
2417   ! nh4cl is active
2418
2419
2420   a =   kg(inh3_g,ibin)
2421   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
2422        + kg(ihcl_g,ibin)*gas(ihcl_g)
2423   c = -(kg(ihcl_g,ibin)*Keq_sg(2))
2424
2425   sfc_a(inh3_g) = quadratic(a,b,c)
2426   sfc_a(ihcl_g) = Keq_sg(2)/sfc_a(inh3_g)
2427
2428   df_gas_s(ihcl_g,ibin) = gas(ihcl_g) - sfc_a(ihcl_g)
2429   df_gas_s(inh3_g,ibin) = gas(inh3_g) - sfc_a(inh3_g)
2430
2431   flux_s(inh3_g,ibin) = kg(inh3_g,ibin)*df_gas_s(inh3_g,ibin)
2432   flux_s(ihcl_g,ibin) = flux_s(ihcl_g,ibin) + flux_s(inh3_g,ibin)
2433
2434   phi_volatile_s(inh3_g,ibin) = phi_nh4cl_s
2435
2436   if(flux_s(ihcl_g,ibin) .gt. 0.0)then
2437      df_gas_s(ihcl_g,ibin) = flux_s(ihcl_g,ibin)/kg(ihcl_g,ibin)        ! recompute df_gas
2438      phi_volatile_s(ihcl_g,ibin) = phi_nh4cl_s
2439   else
2440      sfc_a(ihcl_g)  = gas(ihcl_g) + aer(icl_a,jsolid,ibin)
2441      df_gas_s(ihcl_g,ibin) = -aer(icl_a,jsolid,ibin)
2442      phi_volatile_s(ihcl_g,ibin)=df_gas_s(ihcl_g,ibin)/sfc_a(ihcl_g)  ! not to be used
2443   endif
2444
2445   integrate(inh3_g,jsolid,ibin) = mYES
2446   integrate(ihcl_g,jsolid,ibin) = mYES  ! integrate HCl with explicit euler
2447
2448   ieqblm_ASTEM = mNO
2449
2450   return
2451 end subroutine ASTEM_flux_dry_case3b
2452
2453
2454
2455
2456 !---------------------------------------------------------------------
2457 ! Case 4: NH4NO3 and/or NH4Cl may be active
2458
2459 subroutine ASTEM_flux_dry_case4(ibin, phi_nh4no3_s,phi_nh4cl_s,ieqblm_ASTEM,     &
2460      sfc_a,df_gas_s,flux_s,phi_volatile_s,integrate,kg,gas,electrolyte,epercent, &
2461      Keq_sg,aer)       ! TOUCH
2462
2463   use module_data_mosaic_aero, only: r8, nbin_a_max,                             &
2464        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
2465        nsalt,jsolid,nrxn_aer_sg,naer,                                            &
2466        rtol_eqb_ASTEM,ptol_mol_ASTEM,                                            &
2467        jnh4no3,jnh4cl,ihno3_g,inh3_g,ihcl_g
2468
2469   use module_mosaic_ext, only: quadratic,degas_solid_nh4no3,degas_solid_nh4cl
2470
2471
2472   ! subr arguments
2473   integer, intent(in) :: ibin
2474   integer, intent(inout) :: ieqblm_ASTEM
2475   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2476
2477   real(r8), intent(out) :: phi_nh4no3_s,phi_nh4cl_s
2478   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2479   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2480   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2481   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2482   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2483   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2484   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2485   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2486   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
2487   ! local variables
2488   integer iactive_nh4no3, iactive_nh4cl, iactive, js
2489   real(r8) :: a, b, c
2490   real(r8) :: sum_dum
2491   ! function
2492   !real(r8) :: quadratic
2493
2494
2495
2496   !! EFFI calculate percent composition
2497   sum_dum = 0.0
2498   do js = 1, nsalt
2499      sum_dum = sum_dum + electrolyte(js,jsolid,ibin)
2500   enddo
2501
2502   if(sum_dum .eq. 0.)sum_dum = 1.0
2503
2504   epercent(jnh4no3,jsolid,ibin) = 100.*electrolyte(jnh4no3,jsolid,ibin)/sum_dum
2505   epercent(jnh4cl, jsolid,ibin) = 100.*electrolyte(jnh4cl, jsolid,ibin)/sum_dum
2506   !! EFFI
2507
2508
2509   !-------------------
2510   ! set default values for flags
2511   iactive_nh4no3 = 1
2512   iactive_nh4cl  = 2
2513
2514
2515   ! compute diagnostic products and ratios
2516   phi_nh4no3_s = (gas(inh3_g)*gas(ihno3_g) - Keq_sg(1))/   &
2517        max(gas(inh3_g)*gas(ihno3_g),Keq_sg(1))
2518   phi_nh4cl_s  = (gas(inh3_g)*gas(ihcl_g) - Keq_sg(2))/   &
2519        max(gas(inh3_g)*gas(ihcl_g),Keq_sg(2))
2520
2521
2522   !-------------------
2523   ! now determine if nh4no3 and/or nh4cl are active or significant
2524
2525   ! nh4no3
2526   if( abs(phi_nh4no3_s) .lt. rtol_eqb_ASTEM )then
2527      iactive_nh4no3 = 0
2528   elseif(gas(inh3_g)*gas(ihno3_g) .lt. Keq_sg(1) .and.   &
2529        epercent(jnh4no3,jsolid,ibin) .le. ptol_mol_ASTEM)then
2530      iactive_nh4no3 = 0
2531      if(epercent(jnh4no3,jsolid,ibin) .gt. 0.0)then
2532         call degas_solid_nh4no3(ibin,aer,gas,electrolyte,Keq_sg)
2533      endif
2534   endif
2535
2536   ! nh4cl
2537   if( abs(phi_nh4cl_s) .lt. rtol_eqb_ASTEM )then
2538      iactive_nh4cl = 0
2539   elseif(gas(inh3_g)*gas(ihcl_g) .lt. Keq_sg(2) .and.   &
2540        epercent(jnh4cl, jsolid,ibin) .le. ptol_mol_ASTEM)then
2541      iactive_nh4cl = 0
2542      if(epercent(jnh4cl, jsolid,ibin) .gt. 0.0)then
2543         call degas_solid_nh4cl(ibin,aer,gas,electrolyte,Keq_sg)
2544      endif
2545   endif
2546
2547
2548   iactive = iactive_nh4no3 + iactive_nh4cl
2549
2550   ! check the outcome
2551   if(iactive .eq. 0)return
2552
2553
2554   goto (1,2,3),iactive
2555
2556   !---------------------------------
2557   ! only nh4no3 solid is active
2558 1 call ASTEM_flux_dry_case4a(ibin,phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s,      &
2559        flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)
2560   return
2561
2562
2563   !-----------------
2564   ! only nh4cl solid is active
2565 2 call ASTEM_flux_dry_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,&
2566        phi_volatile_s,integrate,kg,gas,Keq_sg)
2567   return
2568
2569
2570   !-----------------
2571   ! both nh4no3 and nh4cl are active
2572 3 call ASTEM_flux_dry_case4ab(ibin,phi_nh4no3_s,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,  &
2573        df_gas_s,flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)
2574
2575
2576   return
2577 end subroutine ASTEM_flux_dry_case4
2578
2579
2580
2581 !---------------------------------------------------------------------
2582 ! Case 4a
2583
2584 subroutine ASTEM_flux_dry_case4a(ibin, phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s, &
2585      flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg) ! NH4NO3 solid
2586
2587   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2588        jsolid,mYES,mNO, nrxn_aer_sg,                                               &
2589        ihno3_g,inh3_g
2590
2591   use module_mosaic_ext, only: quadratic
2592
2593
2594   ! subr arguments
2595   integer, intent(in) :: ibin
2596   integer, intent(inout) :: ieqblm_ASTEM
2597   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2598
2599   real(r8), intent(in) :: phi_nh4no3_s
2600   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2601   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2602   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2603   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2604   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2605   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2606   ! local variables
2607   real(r8) :: a, b, c
2608   ! function
2609   !real(r8) :: quadratic
2610
2611
2612
2613   a =   kg(inh3_g,ibin)
2614   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
2615        + kg(ihno3_g,ibin)*gas(ihno3_g)
2616   c = -(kg(ihno3_g,ibin)*Keq_sg(1))
2617
2618   sfc_a(inh3_g)  = quadratic(a,b,c)
2619   sfc_a(ihno3_g) = Keq_sg(1)/sfc_a(inh3_g)
2620
2621   integrate(ihno3_g,jsolid,ibin) = mYES
2622   integrate(inh3_g,jsolid,ibin)  = mYES
2623
2624   df_gas_s(ihno3_g,ibin)=gas(ihno3_g)-sfc_a(ihno3_g)
2625   df_gas_s(inh3_g,ibin) =gas(inh3_g) -sfc_a(inh3_g)
2626
2627   phi_volatile_s(ihno3_g,ibin)= phi_nh4no3_s
2628   phi_volatile_s(inh3_g,ibin) = phi_nh4no3_s
2629
2630   flux_s(ihno3_g,ibin) = kg(ihno3_g,ibin)*df_gas_s(ihno3_g,ibin)
2631   flux_s(inh3_g,ibin)  = flux_s(ihno3_g,ibin)
2632
2633   ieqblm_ASTEM = mNO
2634
2635   return
2636 end subroutine ASTEM_flux_dry_case4a
2637
2638
2639
2640 !----------------------------------------------------------------
2641 ! Case 4b
2642
2643 subroutine ASTEM_flux_dry_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,  &
2644      flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg) ! NH4Cl solid
2645
2646   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2647        mYES,jsolid, mNO,nrxn_aer_sg,                                               &
2648        inh3_g,ihcl_g
2649
2650   use module_mosaic_ext, only: quadratic
2651
2652
2653   ! subr arguments
2654   integer, intent(in) :: ibin
2655   integer, intent(inout) :: ieqblm_ASTEM
2656   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2657
2658   real(r8), intent(in) :: phi_nh4cl_s
2659   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2660   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2661   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2662   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2663   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2664   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2665   ! local variables
2666   real(r8) :: a, b, c
2667   ! function
2668   !real(r8) :: quadratic
2669
2670
2671   a =   kg(inh3_g,ibin)
2672   b = - kg(inh3_g,ibin)*gas(inh3_g)   &
2673        + kg(ihcl_g,ibin)*gas(ihcl_g)
2674   c = -(kg(ihcl_g,ibin)*Keq_sg(2))
2675
2676   sfc_a(inh3_g) = quadratic(a,b,c)
2677   sfc_a(ihcl_g) = Keq_sg(2) /sfc_a(inh3_g)
2678
2679   integrate(ihcl_g,jsolid,ibin) = mYES
2680   integrate(inh3_g,jsolid,ibin) = mYES
2681
2682   df_gas_s(ihcl_g,ibin) = gas(ihcl_g)-sfc_a(ihcl_g)
2683   df_gas_s(inh3_g,ibin) = gas(inh3_g)-sfc_a(inh3_g)
2684
2685   phi_volatile_s(ihcl_g,ibin) = phi_nh4cl_s
2686   phi_volatile_s(inh3_g,ibin) = phi_nh4cl_s
2687
2688   flux_s(ihcl_g,ibin) = kg(ihcl_g,ibin)*df_gas_s(ihcl_g,ibin)
2689   flux_s(inh3_g,ibin) = flux_s(ihcl_g,ibin)
2690
2691   ieqblm_ASTEM = mNO
2692
2693   return
2694 end subroutine ASTEM_flux_dry_case4b
2695
2696
2697
2698
2699 !-------------------------------------------------------------------
2700 ! Case 4ab
2701
2702 subroutine ASTEM_flux_dry_case4ab(ibin, phi_nh4no3_s, phi_nh4cl_s,ieqblm_ASTEM,  &
2703      sfc_a,df_gas_s,flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)    ! NH4NO3 + NH4Cl (solid)
2704
2705   use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_aerchtot, ngas_volatile, &
2706        mNO,nrxn_aer_sg,&
2707        ihcl_g,ihno3_g,inh3_g
2708
2709   use module_mosaic_ext, only: quadratic
2710
2711
2712   ! subr arguments
2713   integer, intent(in) :: ibin
2714   integer, intent(inout) :: ieqblm_ASTEM
2715   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2716
2717   real(r8), intent(in) :: phi_nh4no3_s,phi_nh4cl_s
2718   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2719   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
2720   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2721   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,flux_s
2722   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2723   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
2724   ! local variables
2725   real(r8) :: a,b,c,flux_nh3_est, flux_nh3_max, ratio_flux
2726   ! function
2727   !real(r8) :: quadratic
2728
2729   call ASTEM_flux_dry_case4a(ibin,phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s,      &
2730        flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)
2731   call ASTEM_flux_dry_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,       &
2732        flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)
2733
2734
2735   ! estimate nh3 flux and adjust hno3 and/or hcl if necessary
2736
2737   flux_nh3_est = flux_s(ihno3_g,ibin)+flux_s(ihcl_g,ibin)
2738   flux_nh3_max = kg(inh3_g,ibin)*gas(inh3_g)
2739
2740
2741   if(flux_nh3_est .le. flux_nh3_max)then
2742
2743      flux_s(inh3_g,ibin) = flux_nh3_est                 ! all ok - no adjustments needed
2744      sfc_a(inh3_g)       = gas(inh3_g) -                           &  ! recompute sfc_a(ihno3_g)
2745           flux_s(inh3_g,ibin)/kg(inh3_g,ibin)
2746      phi_volatile_s(inh3_g,ibin) = max(abs(phi_nh4no3_s),   &
2747           abs(phi_nh4cl_s))
2748
2749   else                  ! reduce hno3 and hcl flux_ses as necessary so that nh3 flux_s = flux_s_nh3_max
2750
2751      ratio_flux          = flux_nh3_max/flux_nh3_est
2752      flux_s(inh3_g,ibin) = flux_nh3_max
2753      flux_s(ihno3_g,ibin)= flux_s(ihno3_g,ibin)*ratio_flux
2754      flux_s(ihcl_g,ibin) = flux_s(ihcl_g,ibin) *ratio_flux
2755
2756      sfc_a(inh3_g) = 0.0
2757      sfc_a(ihno3_g)= gas(ihno3_g) -                        &  ! recompute sfc_a(ihno3_g)
2758           flux_s(ihno3_g,ibin)/kg(ihno3_g,ibin)
2759      sfc_a(ihcl_g) = gas(ihcl_g)  -                        &  ! recompute sfc_a(ihcl_g)
2760           flux_s(ihcl_g,ibin)/kg(ihcl_g,ibin)
2761
2762      df_gas_s(inh3_g,ibin) =gas(inh3_g) -sfc_a(inh3_g)
2763      df_gas_s(ihno3_g,ibin)=gas(ihno3_g)-sfc_a(ihno3_g)
2764      df_gas_s(ihcl_g,ibin) =gas(ihcl_g) -sfc_a(ihcl_g)
2765
2766      phi_volatile_s(inh3_g,ibin) = max(abs(phi_nh4no3_s),   &
2767           abs(phi_nh4cl_s))
2768
2769
2770   endif
2771
2772   ieqblm_ASTEM = mNO
2773
2774   return
2775 end subroutine ASTEM_flux_dry_case4ab
2776
2777
2778
2779 !=======================================================================
2780 !
2781 ! MIXED-PHASE PARTICLES
2782 !
2783 !***********************************************************************
2784 ! part of ASTEM: computes gas-aerosol fluxes over mixed-phase aerosols
2785 !
2786 ! author: Rahul A. Zaveri
2787 ! update: apr 2006
2788 !-----------------------------------------------------------------------
2789
2790 subroutine ASTEM_flux_mix(ibin, phi_nh4no3_s, phi_nh4cl_s, ieqblm_ASTEM, idry_case3a, &
2791      sfc_a, df_gas_s, df_gas_l, jaerosolstate, flux_s, Heff, phi_volatile_s,            &
2792      phi_volatile_l, integrate, jphase, aer, kg, gas, jhyst_leg, electrolyte, epercent,   &
2793      kel, activity, mc, delta_nh3_max, delta_hno3_max, delta_hcl_max, Keq_nh4cl,        &
2794      Keq_nh4no3, num_a, electrolyte_sum, mass_dry_a, mass_soluble_a, water_a, aH2O,     &
2795      kelvin_nh4no3, kelvin_nh4cl, ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma,      &
2796      nc_Mc, xeq_c, mw_electrolyte, Kp_nh4cl, Kp_nh4no3, Keq_ll, Keq_gl, Keq_sg, MW_c,     &
2797      MW_a, total_species, tot_cl_in, molality0, kappa_nonelectro, mosaic_vars_aa )
2798
2799   use module_data_mosaic_aero, only: r8, nbin_a_max,                             &
2800        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
2801        Ncation, naer, jliquid, nsalt, jsolid, mNO, mYES, jtotal, Nanion, nrxn_aer_gl,      &
2802        nrxn_aer_ll, nrxn_aer_sg,                                                   &
2803        jcaco3, jcacl2, jnacl, ihno3_g, jnh4cl, ihcl_g, inh3_g, jnh4no3, ja_no3, jc_h,      &
2804        ja_cl, jhcl, icl_a, inh4_a, ino3_a, jhno3, mosaic_vars_aa_type
2805
2806   use module_mosaic_ext,  only: compute_activities, ions_to_electrolytes,           &
2807        absorb_tiny_nh4cl, degas_tiny_nh4cl, absorb_tiny_nh4no3, degas_tiny_nh4no3,   &
2808        absorb_tiny_hno3, absorb_tiny_hcl
2809
2810
2811   ! subr arguments
2812   integer, intent(in) :: ibin
2813   integer, intent(inout) :: ieqblm_ASTEM
2814   integer, intent(inout), dimension(nbin_a_max) :: idry_case3a,jaerosolstate
2815   integer, intent(inout), dimension(nbin_a_max) :: jphase,jhyst_leg
2816   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
2817
2818   real(r8), intent(out) :: phi_nh4no3_s, phi_nh4cl_s
2819   real(r8), intent(in) :: aH2O
2820   real(r8), intent(inout) :: Keq_nh4cl,Keq_nh4no3,kelvin_nh4no3,Kp_nh4cl
2821   real(r8), intent(inout) :: kelvin_nh4cl,Kp_nh4no3
2822   real(r8), intent(in), dimension(Ncation) :: zc,MW_c
2823   real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
2824   real(r8), intent(in), dimension(Nanion)  :: za,MW_a
2825   real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
2826   real(r8), intent(inout), dimension(nbin_a_max) :: delta_nh3_max,delta_hno3_max
2827   real(r8), intent(inout), dimension(nbin_a_max) :: delta_hcl_max,mass_soluble_a
2828   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,num_a,mass_dry_a
2829   real(r8), intent(inout), dimension(nbin_a_max) :: gam_ratio
2830   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2831   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, total_species
2832   real(r8), intent(inout) :: tot_cl_in
2833   real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
2834   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
2835   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
2836   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
2837   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
2838   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s
2839   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l
2840   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,Heff
2841   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
2842   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
2843   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg, kel
2844   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)  :: activity,gam
2845   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
2846   real(r8), intent(inout), dimension(Ncation,nbin_a_max) ::mc
2847   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
2848   real(r8), intent(inout), dimension(3,nbin_a_max) :: electrolyte_sum
2849   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2850   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2851   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
2852   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
2853
2854   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
2855
2856   ! local variables
2857   character(len=500) :: tmp_str
2858   integer iv, iadjust, iadjust_intermed, js
2859   real(r8) :: XT,g_nh3_hno3,g_nh3_hcl,a_nh4_no3,a_nh4_cl,a_no3,a_cl,prod_nh4no3
2860   real(r8) :: volatile_cl,sum_dum,prod_nh4cl
2861
2862
2863   call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,   &
2864        nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)    ! for water content calculation
2865   call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte,   &
2866        activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,log_gamZ, &
2867        gam_ratio,Keq_ll,molality0,kappa_nonelectro)
2868
2869   if(water_a(ibin) .eq. 0.0)then
2870      write(tmp_str,*)'Water is zero in liquid phase'
2871      call mosaic_warn_mess(trim(adjustl(tmp_str)))
2872      write(tmp_str,*)'Stopping in ASTEM_flux_wet'
2873      call mosaic_warn_mess(trim(adjustl(tmp_str)))
2874      mosaic_vars_aa%zero_water_flag = .true.
2875   endif
2876
2877
2878   !! EFFI calculate percent composition
2879   sum_dum = 0.0
2880   do js = 1, nsalt
2881      sum_dum = sum_dum + electrolyte(js,jsolid,ibin)
2882   enddo
2883
2884   if(sum_dum .eq. 0.)sum_dum = 1.0
2885
2886   epercent(jcaco3,jsolid,ibin) = 100.*electrolyte(jcaco3,jsolid,ibin)/sum_dum
2887   !! EFFI
2888
2889
2890
2891   !-----------------------------------------------------------------
2892   ! MIXED CASE 1:  caco3 > 0 absorb all acids (and indirectly degas co2)
2893
2894   if(epercent(jcaco3,jsolid,ibin) .gt. 0.0)then
2895      jphase(ibin) = jliquid
2896      call ASTEM_flux_wet_case1(ibin,ieqblm_ASTEM,sfc_a,df_gas_s,flux_s,          &
2897           phi_volatile_s,integrate,jphase,kg,gas,mc,Keq_ll)
2898      return
2899   endif
2900
2901   !-----------------------------------------------------------------
2902   ! MIXED CASE 2: Sulfate-Rich Domain
2903
2904 ! if(XT.lt.1.9999 .and. XT.ge.0.)then   ! excess sulfate (acidic) ! RAZ 11/10/2014
2905   if(XT.lt.2.0 .and. XT.ge.0.)then   ! excess sulfate (acidic) ! RAZ 11/10/2014
2906      jphase(ibin) = jliquid
2907      call ASTEM_flux_wet_case2(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,            &
2908           phi_volatile_l,integrate,gas,kel,mc,water_a,ma,gam,gam_ratio,Keq_ll,   &
2909           Keq_gl)
2910      return
2911   endif
2912
2913   !-------------------------------------------------------------------
2914   ! MIXED CASE 3: hno3 and hcl exchange may happen here and nh4cl may form/evaporate
2915
2916   volatile_cl  = electrolyte(jnacl,jsolid,ibin) +   &
2917        electrolyte(jcacl2,jsolid,ibin)
2918
2919
2920   if(volatile_cl .gt. 0.0 .and. gas(ihno3_g).gt. 0.0 )then
2921
2922      call ASTEM_flux_dry_case3a(ibin,ieqblm_ASTEM,idry_case3a,sfc_a,df_gas_s,    &
2923           flux_s,phi_volatile_s,integrate,aer,kg,gas)
2924
2925      prod_nh4cl = max( (gas(inh3_g)*gas(ihcl_g)-Keq_sg(2)), 0.0d0) +   &
2926           electrolyte(jnh4cl, jsolid,ibin)
2927
2928      if(prod_nh4cl .gt. 0.0)then
2929         call ASTEM_flux_dry_case3b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s, &
2930              flux_s,phi_volatile_s,integrate,aer,kg,gas,electrolyte,epercent,    &
2931              Keq_sg)
2932      endif
2933
2934      jphase(ibin) = jsolid
2935
2936      return
2937   endif
2938
2939   !-------------------------------------------------------------------
2940   ! MIXED CASE 4: nh4no3 or nh4cl or both may be active
2941
2942   if( electrolyte(jnh4no3,jsolid,ibin).gt.0. .and.   &
2943        electrolyte(jnh4cl,jsolid,ibin) .gt.0. )then
2944      jphase(ibin) = jsolid
2945      call ASTEM_flux_dry_case4(ibin,phi_nh4no3_s,phi_nh4cl_s,ieqblm_ASTEM,sfc_a, &
2946           df_gas_s,flux_s,phi_volatile_s,integrate,kg,gas,electrolyte,epercent,  &
2947           Keq_sg,aer)
2948
2949      if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
2950         mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
2951              (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
2952      elseif(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
2953         mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
2954              (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
2955      else
2956         mc(jc_h,ibin) = sqrt(Keq_ll(3))
2957      endif
2958
2959      return
2960
2961   elseif( electrolyte(jnh4no3,jsolid,ibin).gt.0. )then
2962      ! do small adjustments for nh4cl aq
2963      g_nh3_hcl= gas(inh3_g)*gas(ihcl_g)
2964      a_nh4_cl = aer(inh4_a,jliquid,ibin)*aer(icl_a,jliquid,ibin)
2965
2966      iadjust = mNO              ! initialize
2967      if(g_nh3_hcl .gt. 0.0 .and. a_nh4_cl .eq. 0.0)then
2968         call absorb_tiny_nh4cl(ibin,aer,gas,electrolyte,delta_nh3_max,           &
2969              delta_hcl_max,electrolyte_sum)
2970         iadjust = mYES
2971      elseif(g_nh3_hcl .eq. 0.0 .and. a_nh4_cl .gt. 0.0)then
2972         call degas_tiny_nh4cl(ibin,aer,gas,electrolyte)
2973         iadjust = mYES
2974      endif
2975
2976      if(iadjust .eq. mYES)then
2977         call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,   &
2978              na_Ma,nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)      ! update after adjustments
2979         call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,         &
2980              electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,    &
2981              aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                      ! update after adjustments
2982      endif
2983
2984      call ASTEM_flux_mix_case4a(ibin,phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s,   &
2985           df_gas_l,flux_s,Heff,phi_volatile_s,phi_volatile_l,integrate,jphase,kg,&
2986           gas,electrolyte,epercent,kel,activity,mc,Keq_nh4cl,water_a,            &
2987           kelvin_nh4cl,ma,gam,gam_ratio,Kp_nh4cl,Keq_ll,Keq_gl,Keq_sg,aer)      ! nh4no3 solid + nh4cl aq
2988      jphase(ibin) = jtotal
2989      return
2990
2991   elseif( electrolyte(jnh4cl,jsolid,ibin).gt.0.)then
2992      ! do small adjustments for nh4no3 aq
2993      g_nh3_hno3= gas(inh3_g)*gas(ihno3_g)
2994      a_nh4_no3 = aer(inh4_a,jliquid,ibin)*aer(ino3_a,jliquid,ibin)
2995
2996      iadjust = mNO              ! initialize
2997      if(g_nh3_hno3 .gt. 0.0 .and. a_nh4_no3 .eq. 0.0)then
2998         call absorb_tiny_nh4no3(ibin,aer,gas,electrolyte,delta_nh3_max,          &
2999              delta_hno3_max,electrolyte_sum)
3000         iadjust = mYES
3001      elseif(g_nh3_hno3 .eq. 0.0 .and. a_nh4_no3 .gt. 0.0)then
3002         call degas_tiny_nh4no3(ibin,aer,gas,electrolyte)
3003         iadjust = mYES
3004      endif
3005
3006      if(iadjust .eq. mYES)then
3007         call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,   &
3008              na_Ma,nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)      ! update after adjustments
3009         call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,         &
3010              electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,    &
3011              aH2O,ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                      ! update after adjustments
3012      endif
3013
3014      kelvin_nh4no3 = kel(inh3_g,ibin)*kel(ihno3_g,ibin)
3015      Keq_nh4no3 = kelvin_nh4no3*activity(jnh4no3,ibin)*Kp_nh4no3        ! = [NH3]s * [HNO3]s
3016
3017      call ASTEM_flux_mix_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,    &
3018           df_gas_l,flux_s,Heff,phi_volatile_s,phi_volatile_l,integrate,jphase,kg,&
3019           gas,electrolyte,epercent,kel,activity,mc,Keq_nh4no3,water_a,           &
3020           kelvin_nh4no3,ma,gam,gam_ratio,Keq_ll,Keq_gl,Kp_nh4no3,Keq_sg,aer)    ! nh4cl solid + nh4no3 aq
3021      jphase(ibin) = jtotal
3022      return
3023   endif
3024
3025
3026   !-------------------------------------------------------------------
3027
3028   if( (gas(inh3_g)+aer(inh4_a,jliquid,ibin)) .lt. 1.e-25)goto 10  ! no ammonia in the system
3029
3030   !-------------------------------------------------------------------
3031   ! MIXED CASE 5: liquid nh4no3 and/or nh4cl maybe active
3032   ! do some small adjustments (if needed) before deciding case 3
3033
3034   iadjust = mNO         ! default
3035   iadjust_intermed = mNO        ! default
3036
3037   ! nh4no3
3038   g_nh3_hno3 = gas(inh3_g)*gas(ihno3_g)
3039   a_nh4_no3  = aer(inh4_a,jliquid,ibin)*aer(ino3_a,jliquid,ibin)
3040
3041   if(g_nh3_hno3 .gt. 0. .and. a_nh4_no3 .eq. 0.)then
3042      call absorb_tiny_nh4no3(ibin,aer,gas,electrolyte,delta_nh3_max,             &
3043           delta_hno3_max,electrolyte_sum)
3044      iadjust = mYES
3045      iadjust_intermed = mYES
3046   endif
3047
3048   if(iadjust_intermed .eq. mYES)then
3049      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
3050           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
3051      iadjust_intermed = mNO     ! reset
3052   endif
3053
3054   ! nh4cl
3055   g_nh3_hcl = gas(inh3_g)*gas(ihcl_g)
3056   a_nh4_cl  = aer(inh4_a,jliquid,ibin)*aer(icl_a,jliquid,ibin)
3057
3058   if(g_nh3_hcl .gt. 0. .and. a_nh4_cl .eq. 0.)then
3059      call absorb_tiny_nh4cl(ibin,aer,gas,electrolyte,delta_nh3_max,delta_hcl_max,&
3060           electrolyte_sum)
3061      iadjust = mYES
3062      iadjust_intermed = mYES
3063   endif
3064
3065   if(iadjust_intermed .eq. mYES)then
3066      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
3067           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
3068   endif
3069
3070   if(iadjust .eq. mYES)then
3071      call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,electrolyte,&
3072           activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,ma,gam,       &
3073           log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                 ! update after adjustments
3074   endif
3075
3076
3077   ! all adjustments done...
3078
3079   !--------
3080   kelvin_nh4no3 = kel(inh3_g,ibin)*kel(ihno3_g,ibin)
3081   Keq_nh4no3 = kelvin_nh4no3*activity(jnh4no3,ibin)*Kp_nh4no3   ! = [NH3]s * [HNO3]s
3082
3083   kelvin_nh4cl = kel(inh3_g,ibin)*kel(ihcl_g,ibin)
3084   Keq_nh4cl = kelvin_nh4cl*activity(jnh4cl,ibin)*Kp_nh4cl       ! = [NH3]s * [HCl]s
3085
3086   call ASTEM_flux_wet_case3(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,phi_volatile_l,&
3087        integrate,kg,gas,kel,mc,Keq_nh4cl,Keq_nh4no3,water_a,ma,gam,gam_ratio,    &
3088        Keq_ll,Keq_gl,aer,total_species,tot_cl_in,activity,electrolyte)
3089   jphase(ibin) = jliquid
3090
3091   return
3092
3093
3094   !-------------------------------------------------------------------
3095   ! MIXED CASE 6: ammonia = 0. liquid hno3 and hcl exchange may happen here
3096   ! do small adjustments (if needed) before deciding case 4
3097
3098 10 iadjust = mNO                ! default
3099   iadjust_intermed = mNO        ! default
3100
3101   ! hno3
3102   if(gas(ihno3_g).gt.0. .and. aer(ino3_a,jliquid,ibin).eq.0. .and.   &
3103        aer(icl_a,jliquid,ibin) .gt. 0.0)then
3104      call absorb_tiny_hno3(ibin,aer,gas,delta_hno3_max)        ! and degas tiny hcl
3105      iadjust = mYES
3106      iadjust_intermed = mYES
3107   endif
3108
3109   if(iadjust_intermed .eq. mYES)then
3110      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
3111           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
3112      iadjust_intermed = mNO     ! reset
3113   endif
3114
3115   ! hcl
3116   if(gas(ihcl_g).gt.0. .and. aer(icl_a,jliquid,ibin) .eq. 0. .and.   &
3117        aer(ino3_a,jliquid,ibin) .gt. 0.0)then
3118      call absorb_tiny_hcl(ibin,aer,gas,delta_hcl_max)                 ! and degas tiny hno3
3119      iadjust = mYES
3120      iadjust_intermed = mYES
3121   endif
3122
3123   if(iadjust_intermed .eq. mYES)then
3124      call ions_to_electrolytes(jliquid,ibin,XT,aer,electrolyte,zc,za,xeq_a,na_Ma,&
3125           nc_Mc,xeq_c,mw_electrolyte,MW_c,MW_a)         ! update after adjustments
3126   endif
3127
3128   if(iadjust .eq. mYES)then
3129      call compute_activities(ibin,jaerosolstate,jphase,aer,jhyst_leg,            &
3130           electrolyte,activity,mc,num_a,mass_dry_a,mass_soluble_a,water_a,aH2O,  &
3131           ma,gam,log_gamZ,gam_ratio,Keq_ll,molality0,kappa_nonelectro)                 ! update after adjustments
3132   endif
3133
3134   ! all adjustments done...
3135
3136   call ASTEM_flux_wet_case4(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,phi_volatile_l,&
3137        integrate,kg,gas,kel,mc,water_a,ma,gam,Keq_ll,Keq_gl)
3138   jphase(ibin) = jliquid
3139
3140   return
3141 end subroutine ASTEM_flux_mix
3142
3143 !----------------------------------------------------------------------
3144
3145
3146
3147 !------------------------------------------------------------------
3148 ! Mix Case 4a: NH4NO3 solid maybe active. NH4Cl aq maybe active
3149
3150 subroutine ASTEM_flux_mix_case4a(ibin, phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s, &
3151      df_gas_l,flux_s,Heff,phi_volatile_s,phi_volatile_l,integrate,jphase,kg,gas, &
3152      electrolyte,epercent,kel,activity,mc,Keq_nh4cl,water_a,kelvin_nh4cl,ma,gam, &
3153      gam_ratio,Kp_nh4cl,Keq_ll,Keq_gl,Keq_sg,aer)   ! TOUCH
3154
3155   use module_data_mosaic_aero, only: r8, nbin_a_max,                             &
3156        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
3157        Ncation,mYES,jsolid,mNO,jliquid,jtotal,Nanion,nrxn_aer_gl,nrxn_aer_ll,    &
3158        nrxn_aer_sg,naer,                                                         &
3159        rtol_eqb_ASTEM,ptol_mol_ASTEM,                                            &
3160        jnh4no3,ihno3_g,inh3_g,ihcl_g,jnh4cl,ja_no3,jhno3,jc_h,ja_cl,jhcl
3161
3162   use module_mosaic_ext, only: degas_solid_nh4no3
3163
3164
3165   ! subr arguments
3166   integer, intent(in) :: ibin
3167   integer, intent(inout) :: ieqblm_ASTEM
3168   integer, intent(inout), dimension(nbin_a_max) :: jphase
3169   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
3170
3171   real(r8), intent(inout) :: Keq_nh4cl,kelvin_nh4cl,Kp_nh4cl
3172   real(r8), intent(out) :: phi_nh4no3_s
3173   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
3174   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3175   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
3176   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
3177   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
3178   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
3179   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s,df_gas_l
3180   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,Heff
3181   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
3182   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
3183   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg, kel
3184   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)  :: activity,gam
3185   real(r8), intent(inout), dimension(Ncation,nbin_a_max) ::mc
3186   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
3187   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3188   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3189   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
3190   ! local variables
3191   integer iactive_nh4no3, iactive_nh4cl, js
3192   real(r8) :: sum_dum
3193
3194
3195   ! set default values for flags
3196   iactive_nh4no3 = mYES
3197   iactive_nh4cl  = mYES
3198
3199
3200   !! EFFI calculate percent composition
3201   sum_dum = 0.0
3202   do js = 1, nelectrolyte
3203      sum_dum = sum_dum + electrolyte(js,jsolid,ibin)
3204   enddo
3205
3206   if(sum_dum .eq. 0.)sum_dum = 1.0
3207
3208   epercent(jnh4no3,jsolid,ibin) = 100.*electrolyte(jnh4no3,jsolid,ibin)/sum_dum
3209   !! EFFI
3210
3211
3212
3213   ! nh4no3 (solid)
3214   phi_nh4no3_s = (gas(inh3_g)*gas(ihno3_g) - Keq_sg(1))/   &
3215        max(gas(inh3_g)*gas(ihno3_g),Keq_sg(1))
3216
3217   ! nh4cl (liquid)
3218   kelvin_nh4cl = kel(inh3_g,ibin)*kel(ihcl_g,ibin)
3219   Keq_nh4cl = kelvin_nh4cl*activity(jnh4cl,ibin)*Kp_nh4cl       ! = [NH3]s * [HCl]s
3220
3221
3222   !-------------------
3223   ! now determine if nh4no3 and/or nh4cl are active or significant
3224   ! nh4no3 solid
3225   if( abs(phi_nh4no3_s) .le. rtol_eqb_ASTEM )then
3226      iactive_nh4no3 = mNO
3227   elseif(gas(inh3_g)*gas(ihno3_g) .lt. Keq_sg(1) .and.   &
3228        epercent(jnh4no3,jsolid,ibin) .le. ptol_mol_ASTEM)then
3229      iactive_nh4no3 = mNO
3230      if(epercent(jnh4no3,jsolid,ibin) .gt. 0.0)then
3231         call degas_solid_nh4no3(ibin,aer,gas,electrolyte,Keq_sg)
3232      endif
3233   endif
3234
3235   ! nh4cl aq
3236   if( gas(inh3_g)*gas(ihcl_g).eq.0. .or. Keq_nh4cl.eq.0. )then
3237      iactive_nh4cl = mNO
3238   endif
3239
3240
3241   !---------------------------------
3242   if(iactive_nh4no3 .eq. mYES)then
3243
3244      jphase(ibin) = jsolid
3245      call ASTEM_flux_dry_case4a(ibin,phi_nh4no3_s,ieqblm_ASTEM,sfc_a,df_gas_s,   &
3246           flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)      ! NH4NO3 (solid)
3247
3248      if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
3249         mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
3250              (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
3251      elseif(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
3252         mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
3253              (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
3254      else
3255         mc(jc_h,ibin) = sqrt(Keq_ll(3))
3256      endif
3257
3258   endif
3259
3260
3261   if(iactive_nh4cl .eq. mYES)then
3262
3263      jphase(ibin) = jliquid
3264      call ASTEM_flux_wet_case3b(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,           &
3265           phi_volatile_l,integrate,kg,gas,kel,mc,Keq_nh4cl,water_a,ma,gam,       &
3266           gam_ratio,Keq_ll,Keq_gl)        ! NH4Cl (liquid)
3267
3268      if(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
3269         mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
3270              (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
3271      else
3272         mc(jc_h,ibin) = sqrt(Keq_ll(3))
3273      endif
3274
3275   endif
3276
3277
3278   if(iactive_nh4cl .eq. mYES .and. iactive_nh4no3 .eq. mYES)then
3279      jphase(ibin) = jtotal
3280   endif
3281
3282
3283
3284   return
3285 end subroutine ASTEM_flux_mix_case4a
3286
3287
3288
3289 !------------------------------------------------------------------
3290 ! Mix Case 4b: NH4Cl solid maybe active. NH4NO3 aq may or maybe active
3291
3292 subroutine ASTEM_flux_mix_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,   &
3293      df_gas_l,flux_s,Heff,phi_volatile_s,phi_volatile_l,integrate,jphase,kg,gas, &
3294      electrolyte,epercent,kel,activity,mc,Keq_nh4no3,water_a,kelvin_nh4no3,ma,   &
3295      gam,gam_ratio,Keq_ll,Keq_gl,Kp_nh4no3,Keq_sg,aer) ! TOUCH
3296
3297   use module_data_mosaic_aero, only: r8, nbin_a_max,                             &
3298        ngas_aerchtot, ngas_volatile, nelectrolyte,                               &
3299        Ncation,mYES,nsalt,jsolid,mNO,jliquid,jtotal,Nanion,nrxn_aer_gl,naer,     &
3300        nrxn_aer_ll,nrxn_aer_sg,                                                  &
3301        rtol_eqb_ASTEM,ptol_mol_ASTEM,                                            &
3302        jnh4cl,ihcl_g,inh3_g,ihno3_g,jnh4no3,ja_cl,jhcl,jc_h,ja_no3,jhno3
3303
3304   use module_mosaic_ext, only: degas_solid_nh4cl
3305
3306   ! subr arguments
3307   integer, intent(in) :: ibin
3308   integer, intent(inout) :: ieqblm_ASTEM
3309   integer, intent(inout), dimension(nbin_a_max) :: jphase
3310   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
3311
3312   real(r8), intent(inout) :: Keq_nh4no3,kelvin_nh4no3,Kp_nh4no3
3313   real(r8), intent(out) :: phi_nh4cl_s
3314   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,gam_ratio
3315   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3316   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a
3317   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
3318   real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
3319   real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
3320   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_s
3321   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: df_gas_l
3322   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,Heff
3323   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
3324   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
3325   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
3326   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max)  :: activity,gam
3327   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
3328   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
3329   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3330   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3331   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
3332   ! local variables
3333   integer iactive_nh4no3, iactive_nh4cl, js
3334   real(r8) :: sum_dum
3335
3336
3337   ! set default values for flags
3338   iactive_nh4cl  = mYES
3339   iactive_nh4no3 = mYES
3340
3341
3342   !! EFFI calculate percent composition
3343   sum_dum = 0.0
3344   do js = 1, nsalt
3345      sum_dum = sum_dum + electrolyte(js,jsolid,ibin)
3346   enddo
3347
3348   if(sum_dum .eq. 0.)sum_dum = 1.0
3349
3350   epercent(jnh4cl,jsolid,ibin) = 100.*electrolyte(jnh4cl,jsolid,ibin)/sum_dum
3351   !! EFFI
3352
3353
3354   ! nh4cl (solid)
3355   phi_nh4cl_s  = (gas(inh3_g)*gas(ihcl_g) - Keq_sg(2))/   &
3356        max(gas(inh3_g)*gas(ihcl_g),Keq_sg(2))
3357
3358   ! nh4no3 (liquid)
3359   kelvin_nh4no3 = kel(inh3_g,ibin)*kel(ihno3_g,ibin)
3360   Keq_nh4no3 = kelvin_nh4no3*activity(jnh4no3,ibin)*Kp_nh4no3   ! = [NH3]s * [HNO3]s
3361
3362
3363   !-------------------
3364   ! now determine if nh4no3 and/or nh4cl are active or significant
3365   ! nh4cl (solid)
3366   if( abs(phi_nh4cl_s) .le. rtol_eqb_ASTEM )then
3367      iactive_nh4cl = mNO
3368   elseif(gas(inh3_g)*gas(ihcl_g) .lt. Keq_sg(2) .and.   &
3369        epercent(jnh4cl,jsolid,ibin) .le. ptol_mol_ASTEM)then
3370      iactive_nh4cl = mNO
3371      if(epercent(jnh4cl,jsolid,ibin) .gt. 0.0)then
3372         call degas_solid_nh4cl(ibin,aer,gas,electrolyte,Keq_sg)
3373      endif
3374   endif
3375
3376   ! nh4no3 (liquid)
3377   if( gas(inh3_g)*gas(ihno3_g).eq.0. .or. Keq_nh4no3.eq.0. )then
3378      iactive_nh4no3 = mNO
3379   endif
3380
3381
3382   !---------------------------------
3383   if(iactive_nh4cl .eq. mYES)then
3384
3385      jphase(ibin) = jsolid
3386      call ASTEM_flux_dry_case4b(ibin,phi_nh4cl_s,ieqblm_ASTEM,sfc_a,df_gas_s,    &
3387           flux_s,phi_volatile_s,integrate,kg,gas,Keq_sg)    ! NH4Cl (solid)
3388
3389      if(sfc_a(ihcl_g).gt.0.0 .and. ma(ja_cl,ibin).gt.0.0)then
3390         mc(jc_h,ibin) = Keq_gl(4)*sfc_a(ihcl_g)/   &
3391              (kel(ihcl_g,ibin)*gam(jhcl,ibin)**2 * ma(ja_cl,ibin))
3392      elseif(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
3393         mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
3394              (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
3395      else
3396         mc(jc_h,ibin) = sqrt(Keq_ll(3))
3397      endif
3398
3399   endif
3400
3401
3402   if(iactive_nh4no3 .eq. mYES)then
3403
3404      jphase(ibin) = jliquid
3405      call ASTEM_flux_wet_case3a(ibin,ieqblm_ASTEM,sfc_a,df_gas_l,Heff,           &
3406           phi_volatile_l,integrate,kg,gas,kel,mc,Keq_nh4no3,water_a,ma,gam,      &
3407           gam_ratio,Keq_ll,Keq_gl)    ! NH4NO3 (liquid)
3408
3409      if(sfc_a(ihno3_g).gt.0.0 .and. ma(ja_no3,ibin).gt.0.0)then
3410         mc(jc_h,ibin) = Keq_gl(3)*sfc_a(ihno3_g)/   &
3411              (kel(ihno3_g,ibin)*gam(jhno3,ibin)**2 * ma(ja_no3,ibin))
3412      else
3413         mc(jc_h,ibin) = sqrt(Keq_ll(3))
3414      endif
3415
3416   endif
3417
3418
3419   if(iactive_nh4cl .eq. mYES .and. iactive_nh4no3 .eq. mYES)then
3420      jphase(ibin) = jtotal
3421   endif
3422
3423
3424
3425   return
3426 end subroutine ASTEM_flux_mix_case4b
3427
3428
3429
3430 !***********************************************************************
3431 ! part of ASTEM: condenses h2so4, msa, and nh3 analytically over dtchem [s]
3432 !
3433 ! author: Rahul A. Zaveri
3434 ! update: jan 2007
3435 !-----------------------------------------------------------------------
3436
3437 subroutine ASTEM_non_volatiles( dtchem,  jaerosolstate, jphase, &
3438      aer, kg, gas, gas_avg, gas_netprod_otrproc,                                      &
3439      jhyst_leg, electrolyte, epercent, kel, activity, mc, delta_nh3_max,                &
3440      delta_hno3_max, delta_hcl_max, num_a, mass_wet_a, mass_dry_a, mass_soluble_a,     &
3441      vol_dry_a, vol_wet_a, water_a, water_a_hyst, water_a_up, aH2O_a, total_species,    &
3442      tot_cl_in,                                                                   &
3443      aH2O, ma, gam, log_gamZ, zc, za, gam_ratio,  &
3444      xeq_a, na_Ma, nc_Mc, xeq_c, a_zsr, mw_electrolyte, partial_molar_vol, sigma_soln,   &
3445      T_K, RH_pc, mw_aer_mac, dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a, MW_c,       &
3446      growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index,     &
3447      jsulf_poor, jsulf_rich, phi_salt_old,                                   &
3448      kappa_nonelectro, mosaic_vars_aa ) ! TOUCH
3449
3450   use module_data_mosaic_aero,  only: r8, nbin_a_max, nbin_a,   &
3451        ngas_aerchtot, ngas_volatile, nelectrolyte,    &
3452        Ncation, naer, no_aerosol, jtotal, mNO, mYES, Nanion, nrxn_aer_ll, nrxn_aer_sl,    &
3453        nsalt, MDRH_T_NUM, jsulf_poor_NUM, jsulf_rich_NUM,                            &
3454        ih2so4_g, imsa_g, inh3_g, ihno3_g, ihcl_g, iso4_a, imsa_a, jcaco3, jcano3, jnano3,  &
3455        jcacl2, jnacl, inh4_a, mosaic_vars_aa_type
3456
3457   use module_mosaic_ext, only: aerosol_phase_state,conform_electrolytes
3458
3459
3460   !Intent ins
3461   integer, intent(in), dimension(nsalt) :: jsalt_index
3462   integer, intent(in), dimension(jsulf_poor_NUM) :: jsulf_poor
3463   integer, intent(in), dimension(jsulf_rich_NUM) :: jsulf_rich
3464
3465   real(r8), intent(in) :: dtchem
3466   real(r8), intent(in) :: aH2O,T_K,RH_pc,rtol_mesa
3467   real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
3468   real(r8), intent(in), dimension(Ncation) :: zc,MW_c
3469   real(r8), intent(in), dimension(Nanion)  :: za,MW_a
3470   real(r8), intent(in),    dimension(ngas_aerchtot) :: gas_netprod_otrproc
3471   real(r8), intent(in), dimension(ngas_aerchtot) :: partial_molar_vol
3472   real(r8), intent(in), dimension(nelectrolyte) :: mw_electrolyte
3473   real(r8), intent(in), dimension (6,nelectrolyte) :: a_zsr
3474   real(r8), intent(in), dimension(naer) :: kappa_nonelectro
3475
3476   !Intent-inouts
3477   integer, intent(inout), dimension(nsalt) :: jsalt_present
3478   integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase,jhyst_leg
3479
3480   real(r8), intent(inout) :: sigma_water
3481   real(r8), intent(inout), dimension(Ncation) :: nc_Mc,xeq_c
3482   real(r8), intent(inout), dimension(Nanion)  :: xeq_a,na_Ma
3483   real(r8), intent(inout), dimension(nbin_a_max) :: num_a,mass_soluble_a,vol_dry_a
3484   real(r8), intent(inout), dimension(nbin_a_max) :: water_a,delta_nh3_max
3485   real(r8), intent(inout), dimension(nbin_a_max) :: delta_hno3_max,delta_hcl_max
3486   real(r8), intent(inout), dimension(nbin_a_max) :: water_a_hyst,water_a_up,aH2O_a
3487   real(r8), intent(inout), dimension(nbin_a_max) :: mass_wet_a,mass_dry_a
3488   real(r8), intent(inout), dimension(nbin_a_max) :: growth_factor,MDRH
3489   real(r8), intent(inout), dimension(nbin_a_max) :: vol_wet_a,gam_ratio,sigma_soln
3490   real(r8), intent(inout), dimension(ngas_volatile) :: total_species
3491   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3492   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas_avg  ! average gas conc. over dtchem time step (nmol/m3)
3493
3494   type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
3495
3496   ! gas_netprod_otrproc = gas net production rate from other processes
3497   !    such as gas-phase chemistry and emissions (nmol/m3/s)
3498   real(r8), intent(inout) :: tot_cl_in
3499   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
3500   real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
3501   real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
3502   real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
3503   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg,kel
3504   real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
3505   real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
3506   real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
3507   real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
3508   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3509   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
3510   real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: epercent
3511   real(r8), intent(inout), dimension(nsalt) :: phi_salt_old
3512
3513   !Local variables
3514   integer ibin,iupdate_phase_state
3515   real(r8) :: decay_h2so4,decay_msa,delta_h2so4,delta_tmsa,delta_nh3,delta_hno3
3516   real(r8) :: delta_hcl,XT,sumkg_h2so4,sumkg_msa,sumkg_nh3,sumkg_hno3,sumkg_hcl
3517   real(r8) :: tmp_kxt, tmp_kxt2, tmp_pok, tmp_pxt, tmp_q1, tmp_q3, tmp_q4
3518   real(r8), dimension(nbin_a) :: delta_so4,delta_msa,delta_nh4
3519   real(r8), dimension(nbin_a) :: new_so4a, old_so4a !BALLI for debugging only
3520
3521   sumkg_h2so4 = 0.0
3522   sumkg_msa   = 0.0
3523   sumkg_nh3   = 0.0
3524   sumkg_hno3  = 0.0
3525   sumkg_hcl   = 0.0
3526   do ibin = 1, nbin_a
3527      sumkg_h2so4 = sumkg_h2so4 + kg(ih2so4_g,ibin)
3528      sumkg_msa   = sumkg_msa   + kg(imsa_g,ibin)
3529      sumkg_nh3   = sumkg_nh3   + kg(inh3_g,ibin)
3530      sumkg_hno3  = sumkg_hno3  + kg(ihno3_g,ibin)
3531      sumkg_hcl   = sumkg_hcl   + kg(ihcl_g,ibin)
3532   enddo
3533
3534
3535
3536   !--------------------------------------
3537   ! H2SO4
3538   tmp_q1 = gas(ih2so4_g)
3539   tmp_pxt = max( gas_netprod_otrproc(ih2so4_g)*dtchem, 0.0_r8 )
3540   tmp_kxt = sumkg_h2so4*dtchem
3541   old_so4a(1:nbin_a) = aer(iso4_a,jtotal,1:nbin_a)   ! added for debug REMOVE IT BALLI AFTER DEBUG
3542   if ( (tmp_q1+tmp_pxt > 1.e-14_r8) .and. &
3543        (tmp_kxt >= 1.0e-20_r8) ) then
3544
3545 !    ! integrate h2so4 condensation analytically
3546 !    decay_h2so4   = exp(-sumkg_h2so4*dtchem)
3547 !    delta_h2so4   = gas(ih2so4_g)*(1.0 - decay_h2so4)
3548 !    gas(ih2so4_g) = gas(ih2so4_g)*decay_h2so4
3549
3550      ! integrate h2so4 condensation + gas-phase production analytically
3551      ! tmp_q1 = mix-rat at t=tcur
3552      ! tmp_q3 = mix-rat at t=tcur+dtchem
3553      ! tmp_q4 = avg mix-rat between t=tcur and t=tcur+dtchem
3554      if (tmp_kxt > 0.001_r8) then
3555         ! use analytical exponential expression
3556         tmp_pok = tmp_pxt/tmp_kxt
3557         tmp_q3 = (tmp_q1 - tmp_pok)*exp(-tmp_kxt) + tmp_pok
3558         tmp_q4 = (tmp_q1 - tmp_pok)*(1.0_r8 - exp(-tmp_kxt))/tmp_kxt + tmp_pok
3559      else
3560         ! use taylors series expansion
3561         tmp_kxt2 = tmp_kxt*tmp_kxt
3562         tmp_q3 = tmp_q1 *(1.0_r8 - tmp_kxt        + tmp_kxt2*0.5_r8) &
3563                + tmp_pxt*(1.0_r8 - tmp_kxt*0.5_r8 + tmp_kxt2/6.0_r8)
3564         tmp_q4 = tmp_q1 *(1.0_r8 - tmp_kxt*0.5_r8 + tmp_kxt2/6.0_r8) &
3565                + tmp_pxt*(0.5_r8 - tmp_kxt/6.0_r8 + tmp_kxt2/24.0_r8)
3566      end if
3567      gas(ih2so4_g) = tmp_q3
3568      gas_avg(ih2so4_g) = tmp_q4
3569      delta_h2so4 = (tmp_q1 + tmp_pxt) - tmp_q3   ! this is the change due to condensation
3570
3571
3572      ! now distribute delta_h2so4 to each bin and conform the particle (may degas by massbal)
3573      do ibin = 1, nbin_a
3574         if(jaerosolstate(ibin) .ne. no_aerosol)then
3575            delta_so4(ibin) = delta_h2so4*kg(ih2so4_g,ibin)/sumkg_h2so4
3576            aer(iso4_a,jtotal,ibin) = aer(iso4_a,jtotal,ibin) +   &
3577                 delta_so4(ibin)
3578         endif
3579      enddo
3580
3581   else
3582      ! h2so4 conc. (after production) is negligible OR
3583      !    uptake by aerosols is negligible
3584      ! in this case, update gas conc. (production) but do not bother to update aerosol conc.
3585      gas(ih2so4_g) = tmp_q1 + tmp_pxt
3586      gas_avg(ih2so4_g) = tmp_q1 + tmp_pxt*0.5_r8
3587      delta_h2so4 = 0.0
3588      do ibin = 1, nbin_a
3589         delta_so4(ibin) = 0.0
3590      enddo
3591
3592   endif
3593 ! debug output (Remove this BALLI after debugging)
3594   new_so4a(1:nbin_a) = aer(iso4_a,jtotal,1:nbin_a)   ! added for debug
3595   ! h2so4 condensation is now complete
3596   !--------------------------------------
3597
3598
3599
3600   ! MSA
3601   if(gas(imsa_g) .gt. 1.e-14)then
3602
3603      ! integrate msa condensation analytically
3604      decay_msa   = exp(-sumkg_msa*dtchem)
3605      delta_tmsa  = gas(imsa_g)*(1.0 - decay_msa)
3606      gas(imsa_g) = gas(imsa_g)*decay_msa
3607
3608      ! now distribute delta_msa to each bin and conform the particle (may degas by massbal)
3609      do ibin = 1, nbin_a
3610         if(jaerosolstate(ibin) .ne. no_aerosol)then
3611            delta_msa(ibin) = delta_tmsa*kg(imsa_g,ibin)/sumkg_msa
3612            aer(imsa_a,jtotal,ibin) = aer(imsa_a,jtotal,ibin) +   &
3613                 delta_msa(ibin)
3614         endif
3615      enddo
3616
3617   else
3618
3619      delta_tmsa = 0.0
3620      do ibin = 1, nbin_a
3621         delta_msa(ibin) = 0.0
3622      enddo
3623
3624   endif
3625   ! msa condensation is now complete
3626   !-------------------------------------
3627
3628
3629
3630   ! compute max allowable nh3, hno3, and hcl condensation
3631   delta_nh3 = gas(inh3_g) *(1.0 - exp(-sumkg_nh3*dtchem))
3632   delta_hno3= gas(ihno3_g)*(1.0 - exp(-sumkg_hno3*dtchem))
3633   delta_hcl = gas(ihcl_g) *(1.0 - exp(-sumkg_hcl*dtchem))
3634
3635   ! compute max possible nh4 condensation for each bin
3636   do ibin = 1, nbin_a
3637      if(jaerosolstate(ibin) .ne. no_aerosol)then
3638         delta_nh3_max(ibin) = delta_nh3*kg(inh3_g,ibin)/sumkg_nh3
3639         delta_hno3_max(ibin)= delta_hno3*kg(ihno3_g,ibin)/sumkg_hno3
3640         delta_hcl_max(ibin) = delta_hcl*kg(ihcl_g,ibin)/sumkg_hcl
3641      endif
3642   enddo
3643
3644
3645   if(delta_h2so4 .eq. 0.0 .and. delta_tmsa .eq. 0.0)then
3646      iupdate_phase_state = mNO
3647      goto 100
3648   endif
3649
3650
3651   ! now condense appropriate amounts of nh3 to each bin  EFFI
3652   do ibin = 1, nbin_a
3653
3654      if(electrolyte(jnacl,jtotal,ibin)  .eq. 0.0 .and.   &
3655           electrolyte(jcacl2,jtotal,ibin) .eq. 0.0 .and.   &
3656           electrolyte(jnano3,jtotal,ibin) .eq. 0.0 .and.   &
3657           electrolyte(jcano3,jtotal,ibin) .eq. 0.0 .and.   &
3658           electrolyte(jcaco3,jtotal,ibin) .eq. 0.0 .and.   &
3659           jaerosolstate(ibin) .ne. no_aerosol)then
3660
3661         delta_nh4(ibin) = min( (2.*delta_so4(ibin)+delta_msa(ibin)),   &
3662              delta_nh3_max(ibin) )
3663
3664         aer(inh4_a,jtotal,ibin) = aer(inh4_a,jtotal,ibin) +        &  ! update aer-phase
3665              delta_nh4(ibin)
3666
3667         gas(inh3_g) = gas(inh3_g) - delta_nh4(ibin)             ! update gas-phase
3668
3669      else
3670
3671         delta_nh4(ibin)     = 0.0
3672
3673      endif
3674
3675   enddo
3676
3677   iupdate_phase_state = mYES
3678
3679
3680   ! recompute phase equilibrium
3681 100 if(iupdate_phase_state .eq. mYES)then
3682      do ibin = 1, nbin_a
3683         if(jaerosolstate(ibin) .ne. no_aerosol)then
3684            call conform_electrolytes(jtotal, ibin, XT, aer, gas, electrolyte,          &
3685                 total_species, tot_cl_in)
3686            call aerosol_phase_state( ibin, jaerosolstate,                     &
3687                 jphase, aer, jhyst_leg, electrolyte, epercent, kel, activity, mc, num_a, &
3688                 mass_wet_a, mass_dry_a, mass_soluble_a, vol_dry_a, vol_wet_a, water_a,   &
3689                 water_a_hyst, water_a_up, aH2O_a, aH2O,                  &
3690                 ma, gam, log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc,   &
3691                 xeq_c,        mw_electrolyte, partial_molar_vol, sigma_soln, T_K,        & ! RAZ deleted a_zsr
3692                 RH_pc, mw_aer_mac, dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a,      &
3693                 MW_c, growth_factor, MDRH, MDRH_T, molality0, rtol_mesa, jsalt_present,  &
3694                 jsalt_index, jsulf_poor, jsulf_rich, phi_salt_old,            &
3695                 kappa_nonelectro, mosaic_vars_aa )
3696         endif
3697      enddo
3698   endif
3699
3700   return
3701 end subroutine ASTEM_non_volatiles
3702
3703
3704
3705 !=================================================================
3706 ! SOA module
3707
3708 !***********************************************************************
3709 ! part of ASTEM: condenses secondary organic species over TSI time interval
3710 ! mechanism adapted from SORGAM
3711 !
3712 ! author: Rahul A. Zaveri
3713 ! update: apr 2005
3714 !-----------------------------------------------------------------------
3715 subroutine ASTEM_secondary_organics(dtchem, jaerosolstate,sfc_a,Heff,            &
3716      phi_volatile_l,integrate,aer,kg,gas,sat_soa,total_species)
3717
3718   use module_data_mosaic_aero, only: nbin_a_max, nbin_a, naer, no_aerosol,   &
3719        ngas_aerchtot, ngas_volatile, jtotal,mYES,                            &
3720        isoa_first
3721
3722
3723   ! subr arguments
3724   integer, intent(in), dimension(nbin_a_max) :: jaerosolstate
3725   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
3726
3727   real(r8), intent(in) :: dtchem
3728   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3729   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, sat_soa, total_species
3730   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
3731   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: Heff
3732   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
3733   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3734   ! local variables
3735   integer ibin, iv, jp,ieqblm, nsteps_max,ieqblm_soa,isteps_SOA
3736   parameter(nsteps_max = 400)
3737   real(r8) :: dtmax, t_new, t_old, t_out
3738   real(r8) :: sum1, sum2
3739
3740
3741   ! initialize time
3742   t_old = 0.0
3743   t_out = dtchem
3744   isteps_SOA = 0
3745
3746
3747   do iv = isoa_first, ngas_volatile
3748      total_species(iv) = gas(iv)
3749      do ibin = 1, nbin_a
3750         if (jaerosolstate(ibin) .eq. no_aerosol) cycle
3751         total_species(iv) = total_species(iv) + aer(iv,jtotal,ibin)
3752      enddo
3753   enddo
3754
3755
3756
3757   ! overall integration loop begins over dtchem seconds
3758 10 isteps_SOA = isteps_SOA + 1
3759
3760   ! compute new fluxes
3761   ieqblm_soa = mYES                     ! reset to default
3762
3763   do 501 ibin = 1, nbin_a
3764      if (jaerosolstate(ibin) .eq. no_aerosol) goto 501
3765
3766      call ASTEM_flux_soa(ibin,sfc_a,Heff,integrate,aer,gas,sat_soa,ieqblm_soa)
3767
3768 501 continue
3769   if(ieqblm_soa .eq. mYES)goto 30 ! all bins have reached equilibrium
3770
3771   !-----------------------
3772
3773
3774 ! calculate maximum possible internal time-step
3775 11 call ASTEM_dtmax_soa(dtchem, dtmax, phi_volatile_l,integrate,kg)
3776   t_new = t_old + dtmax ! update time
3777   if(t_new .gt. t_out)then      ! check if the new time step is too large
3778      dtmax = t_out - t_old
3779      t_new = t_out*1.01
3780   endif
3781
3782
3783
3784
3785   !------------------------------------------
3786   ! do internal time-step (dtmax) integration
3787
3788   jp = jtotal
3789
3790   do 20 iv = isoa_first, ngas_volatile
3791
3792      sum1 = 0.0
3793      sum2 = 0.0
3794
3795      do 21 ibin = 1, nbin_a
3796         if(jaerosolstate(ibin) .eq. no_aerosol)goto 21
3797
3798         sum1 = sum1 + aer(iv,jp,ibin)/   &
3799              (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin)*integrate(iv,jp,ibin))
3800         sum2 = sum2 + kg(iv,ibin)*integrate(iv,jp,ibin)/   &
3801              (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin)*integrate(iv,jp,ibin))
3802
3803 21   continue
3804
3805      ! first update gas concentration
3806      gas(iv) = (total_species(iv) - sum1)/   &
3807           (1. + dtmax*sum2)
3808
3809      ! now update aer concentration in the jp phase
3810      do 22 ibin = 1, nbin_a
3811         if (jaerosolstate(ibin) .eq. no_aerosol) goto 22
3812
3813         if(integrate(iv,jp,ibin) .eq. mYES)then
3814            aer(iv,jp,ibin) =   &
3815                 (aer(iv,jp,ibin) + dtmax*kg(iv,ibin)*gas(iv))/   &
3816                 (1. + dtmax*kg(iv,ibin)*Heff(iv,ibin))
3817         endif
3818
3819 22   continue
3820
3821 20 continue
3822   !------------------------------------------
3823   ! sub-step integration done
3824
3825
3826   ! update jtotal
3827   !      do iv = isoa_first, ngas_volatile
3828   !        aer(iv,jtotal,ibin)=aer(iv,jsolid,ibin)+aer(iv,jliquid,ibin)
3829   !      enddo
3830
3831
3832   ! update time
3833   t_old = t_new
3834
3835   if(t_new .lt. 0.9999*t_out) goto 10
3836   !================================================
3837   ! end of integration
3838
3839 30 continue
3840
3841
3842   return
3843 end subroutine ASTEM_secondary_organics
3844
3845
3846
3847 !***********************************************************************
3848 ! part of ASTEM: computes fluxes of soa species
3849 !
3850 ! author: Rahul A. Zaveri
3851 ! update: apr 2005
3852 !-----------------------------------------------------------------------
3853 subroutine ASTEM_flux_soa(ibin,sfc_a,Heff,integrate,aer,gas,sat_soa,ieqblm_soa)         ! TOUCH
3854
3855   use module_data_mosaic_aero, only: r8, ngas_aerchtot, ngas_volatile, &
3856        nbin_a_max,naer,mNO,jtotal, rtol_eqb_ASTEM,                     &
3857        ioc_a,isoa_first
3858
3859
3860   ! subr arguments
3861   integer, intent(in) :: ibin
3862   integer, intent(inout) :: ieqblm_soa
3863   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
3864
3865   real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
3866   real(r8), intent(inout), dimension(ngas_volatile) :: sfc_a, sat_soa
3867   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: Heff
3868   real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
3869   ! local variables
3870   integer iv, jp
3871   real(r8) :: dum, sum_dum, sum_soa, small_oc
3872   real(r8), dimension(ngas_volatile,nbin_a_max) :: df_gas_o,flux_o,phi_volatile_o
3873
3874   small_oc  = 1.e-15            ! ng/m^3
3875
3876
3877   ! default fluxes and other stuff
3878   do iv = isoa_first, ngas_volatile
3879      sfc_a(iv)               = gas(iv)
3880      df_gas_o(iv,ibin)       = 0.0
3881      flux_o(iv,ibin)         = 0.0
3882      phi_volatile_o(iv,ibin) = 0.0
3883   enddo
3884
3885
3886   jp = jtotal
3887
3888   ! compute mole fractions of soa species
3889   sum_soa = 0.0
3890   do iv = isoa_first, ngas_volatile
3891      sum_soa = sum_soa + aer(iv,jp,ibin)
3892   enddo
3893   sum_soa = sum_soa + aer(ioc_a,jp,ibin)/200.  ! 200 is assumed MW of primary OC
3894
3895
3896   ! check threshold concentration for SOA formation in the absence of primary OC
3897   if(aer(ioc_a,jp,ibin) .eq. 0.0)then
3898      sum_dum = 0.0
3899      do iv = isoa_first, ngas_volatile
3900         sum_dum = sum_dum + (gas(iv)+aer(iv,jp,ibin))/sat_soa(iv)
3901      enddo
3902
3903      if(sum_dum .le. 1.0)then   ! transfer all aer to gas and quit
3904         do iv = isoa_first, ngas_volatile
3905            gas(iv)         = gas(iv) + aer(iv,jp,ibin)
3906            aer(iv,jp,ibin) = 0.0
3907            integrate(iv,jp,ibin) = 0.0
3908         enddo
3909         return
3910      endif
3911
3912      sum_soa = max(sum_soa, 1.d-10)
3913
3914   endif
3915
3916
3917
3918
3919   ! compute Heff
3920   do iv = isoa_first, ngas_volatile
3921
3922      Heff(iv,ibin) = sat_soa(iv)/sum_soa
3923      sfc_a(iv) = aer(iv,jp,ibin)*Heff(iv,ibin)          ! nmol/m^3
3924      df_gas_o(iv,ibin) = gas(iv) - sfc_a(iv)
3925
3926      dum = max(sfc_a(iv),gas(iv))
3927      if(dum .gt. 0.0)then
3928         phi_volatile_o(iv,ibin) = df_gas_o(iv,ibin)/dum
3929      else
3930         phi_volatile_o(iv,ibin) = 0.0
3931      endif
3932
3933      ! check equilibrium
3934      if(abs(phi_volatile_o(iv,ibin)) .le. rtol_eqb_ASTEM)then
3935         integrate(iv,jp,ibin) = 0.0
3936      else
3937         integrate(iv,jp,ibin) = 1.0
3938         ieqblm_soa = mNO
3939      endif
3940
3941   enddo
3942
3943
3944   return
3945 end subroutine ASTEM_flux_soa
3946
3947
3948
3949 !***********************************************************************
3950 ! part of ASTEM: computes fluxes of soa species
3951 !
3952 ! author: Rahul A. Zaveri
3953 ! update: apr 2005
3954 !-----------------------------------------------------------------------
3955 subroutine ASTEM_dtmax_soa(dtchem, dtmax, phi_volatile_l,integrate,kg)          ! TOUCH
3956
3957   use module_data_mosaic_aero, only: r8, ngas_aerchtot, ngas_volatile, &
3958        nbin_a_max,jtotal,mYES, alpha_astem,nbin_a,                     &
3959        isoa_first
3960
3961
3962   ! subr arguments
3963   integer, intent(inout), dimension(ngas_volatile,3,nbin_a_max) :: integrate
3964
3965   real(r8), intent(in)  :: dtchem
3966   real(r8), intent(out) :: dtmax
3967   real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
3968   real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
3969
3970   ! local variables
3971   character(len=500) :: tmp_str
3972   integer ibin, iv, jp
3973   real(r8) :: h_gas, h_gas_i(ngas_volatile), h_sub_max,   &
3974        sum_kg_phi
3975
3976
3977   h_sub_max = dtchem/6. ! sec
3978
3979   jp = jtotal
3980
3981   ! GAS-SIDE
3982   ! calculate h_gas_i and h_gas
3983
3984   h_gas = 2.e16
3985
3986   do 6 iv = isoa_first, ngas_volatile
3987
3988      h_gas_i(iv) = 1.e16
3989      sum_kg_phi = 0.0
3990
3991      do ibin = 1, nbin_a
3992         if(integrate(iv,jtotal,ibin) .eq. mYES)then
3993            sum_kg_phi = sum_kg_phi +   &
3994                 abs(phi_volatile_l(iv,ibin))*kg(iv,ibin)
3995         endif
3996      enddo
3997
3998      if(sum_kg_phi .gt. 0.0)then
3999         h_gas_i(iv) = alpha_astem/sum_kg_phi
4000         h_gas       = min(h_gas, h_gas_i(iv))
4001      endif
4002
4003 6 continue
4004
4005
4006   dtmax = min(h_gas, h_sub_max)
4007
4008
4009   if(dtmax .le. 1.0e-10)then
4010      write(tmp_str,*)' SOA dtmax = ', dtmax
4011      call mosaic_warn_mess(trim(adjustl(tmp_str)))
4012   endif
4013
4014
4015   return
4016 end subroutine ASTEM_dtmax_soa
4017
4018
4019
4020 end module module_mosaic_astem