chem/module_mosaic_box_aerchem.F

   1 module module_mosaic_box_aerchem
   2
   3 use module_data_mosaic_kind, only: r8
   4
   5 implicit none
   6
   7 contains
   8   ! zz01aerchemistry.f (mosaic.25.0)
   9   !********************************************************************************************
  10   !   code history
  11   !   6/3/2015 RAZ  - bound temperature between 220 K and 330 K
  12   !   6/3/2015 RAZ  - bound drh_mutual between 0% and 100%
  13   !   01-may-07 raz - updated CRH and hysteresis treatment for cano3 and cacl2 salts
  14   !   09-jan-07 raz - major clean up of variables and subroutines
  15   !   25-sep-06 raz - added kelvin effect treatment for condensing species
  16   !   22-sep-06 raz - changed "min" to "max" in ratio_AN and ratio_AC definitions
  17   !   21-jul-06 raz - revised and debugged kelvin effect algorithm
  18   !   17-jun-06 raz - added MSA chemistry in particle phase
  19   !   06-jan-05 raz - implemented revised ASTEM algorithm
  20   !   08-oct-05 raz - debugged
  21   !   21-sep-05 raz - revised adaptive time stepping scheme in MESA.
  22   !   28-apr-05 raz - reversed calls to form_cacl2 and form_nacl
  23   !                   fixed caco3 error in subr. electrolytes_to_ions
  24   !                   renamed dens_aer to dens_aer_mac; mw_aer to mw_aer_mac
  25   !   27-apr-05 raz - updated dry_mass calculation approach in MESA_convergence
  26   !   22-apr-05 raz - fixed CaSO4 mass balance problem and updated algorithm to
  27   !                   calculate phi_volatile for nh3, hno3, and hcl.
  28   !   20-apr-05 raz - updated ASCEEM
  29   !   19-apr-05 raz - updated the algorithm to constrain the nh4 concentration
  30   !                   during simultaneous nh3, hno3, and hcl integration such
  31   !                   that it does not exceed the max possible value for a given bin
  32   !   14-apr-05 raz - fixed ASTEM_flux_wet_case3 and ASTEM_flux_dry_case3c
  33   !   11-apr-05 raz - added SOA based on SORGAM mechanism
  34   !   11-jan-05 raz - major updates to many subroutines
  35   !   18-nov-04 rce - make sure that acos argument is between +/-1.0
  36   !   28-jan-04 rce - added subr aerchem_boxtest_output;
  37   !       eliminated some unnecessary "include v33com-"
  38   !   01-dec-03 rce - added "implicit none" to many routines;
  39   !       eliminated some unnecessary "include v33com-"
  40   !   05-oct-03 raz - added hysteresis treatment
  41   !   02-sep-03 raz - implemented ASTEM
  42   !   10-jul-03 raz - changed ix to ixd in interp. subrs fast*_up and fast*_lo
  43   !   08-jul-03 raz - implemented ASTEM (adaptive step time-split
  44   !                   explicit euler method)
  45   !   26-jun-03 raz - updated almost all the subrs. this version contains
  46   !       options for rigorous and fast solvers (including lsode solver)
  47   !
  48   !   07-oct-02 raz - made zx and zm integers in activity coeff subs.
  49   !   16-sep-02 raz - updated many subrs to treat calcium salts
  50   !   19-aug-02 raz - inlcude v33com9a in subr aerosolmtc
  51   !   14-aug-02 rce - "(msectional.eq.0)" changed to "(msectional.le.0)"
  52   !   07-aug-02 rce - this is rahul's latest version from freshair
  53   !       AFTER adding "real mean_molecular_speed" wherever it is used
  54   !   01-apr-02 raz - made final tests and gave the code to jerome
  55   !
  56   !   04--14-dec-01 rce - several minor changes during initial testing/debug
  57   !       in 3d los angeles simulation
  58   !       (see earlier versions for details about these changes)
  59   !-----------------------------------------------------------------------
  60   !23456789012345678901234567890123456789012345678901234567890123456789012
  61
  62   !***********************************************************************
  63   ! MOSAIC (Model for Simulating Aerosol Interactions and Chemistry)
  64   !
  65   ! author: Rahul A. Zaveri
  66   ! update: dec 2004
  67   !-----------------------------------------------------------------------
  68
  69   subroutine mosaic_box_aerchemistry(        aH2O,               T_K,            &!Intent-ins
  70        P_atm,                  RH_pc,        dtchem,                             &
  71        mcall_load_mosaic_parameters,         mcall_print_aer_in, sigmag_a,       &
  72        kappa_nonelectro,                                                         &
  73        jaerosolstate,          aer,                                              &!Intent-inouts
  74        num_a,                  water_a,      gas,                                &
  75        gas_avg,                gas_netprod_otrproc,              Dp_dry_a,       &
  76        dp_wet_a,               jhyst_leg,                                        &
  77        mosaic_vars_aa,                                                           &
  78        mass_dry_a_bgn,         mass_dry_a,                                       &!Intent-outs
  79        dens_dry_a_bgn,         dens_dry_a,   water_a_hyst,       aH2O_a,         &
  80        uptkrate_h2so4,         gam_ratio,    jaerosolstate_bgn                   )
  81
  82     use module_data_mosaic_aero, only:                                             &
  83          nbin_a_max, ngas_aerchtot, ngas_volatile, naer, nsalt,                    &!Parameters
  84          Nanion, Ncation, nrxn_aer_sl, nrxn_aer_ll, nrxn_aer_gl, nrxn_aer_sg,      &!Parameters
  85          MDRH_T_NUM, nelectrolyte,                                                 &!Parameters
  86          jsalt_index, jsulf_poor, jsulf_rich, rtol_mesa, dens_aer_mac,             &
  87          mw_aer_mac, zc, MW_c, za, MW_a, mw_comp_a, dens_comp_a, b_zsr,aw_min,     &
  88          mw_electrolyte, partial_molar_vol, a_zsr, d_mdrh, b_mtem, ref_index_a,    &
  89          Nmax_mesa, nmax_ASTEM, mosaic_vars_aa_type
  90
  91     implicit none
  92
  93     !Intent-ins
  94     integer, intent(in) :: mcall_load_mosaic_parameters, mcall_print_aer_in
  95
  96     real(r8), intent(in) :: aH2O
  97     real(r8), intent(in) :: T_K, P_atm, RH_pc
  98     real(r8), intent(in) :: dtchem
  99
 100     real(r8), intent(in), dimension(nbin_a_max)        :: sigmag_a
 101     real(r8), intent(in), dimension(naer)              :: kappa_nonelectro
 102
 103     !Intent-inouts
 104     integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate
 105     integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
 106
 107     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 108     real(r8), intent(inout), dimension(nbin_a_max)        :: num_a
 109     real(r8), intent(inout), dimension(nbin_a_max)        :: water_a
 110     real(r8), intent(inout), dimension(ngas_aerchtot)     :: gas
 111     real(r8), intent(inout), dimension(ngas_aerchtot)     :: gas_avg  ! average gas conc. over dtchem time step (nmol/m3)
 112     real(r8), intent(in),    dimension(ngas_aerchtot)     :: gas_netprod_otrproc
 113               ! gas_netprod_otrproc = gas net production rate from other processes
 114               !    such as gas-phase chemistry and emissions (nmol/m3/s)
 115               ! this allows the condensation (gasaerexch) routine to apply production and condensation loss
 116               !    together, which is more accurate numerically
 117               ! NOTE - must be >= zero, as numerical method can fail when it is negative
 118               ! NOTE - currently for mosaic, only the value for h2so4 can be non-zero
 119     real(r8), intent(inout), dimension(nbin_a_max)        :: Dp_dry_a, dp_wet_a
 120
 121     ! note - purpose of this data structure is to simplify passing new variables
 122     !        into and out of the many mosaic routines
 123     type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
 124
 125     !Intent-outs
 126     integer, intent(out), dimension(nbin_a_max) :: jaerosolstate_bgn
 127
 128     real(r8), intent(out), dimension(nbin_a_max) :: mass_dry_a_bgn
 129     real(r8), intent(out), dimension(nbin_a_max) :: mass_dry_a
 130     real(r8), intent(out), dimension(nbin_a_max) :: dens_dry_a_bgn
 131     real(r8), intent(out), dimension(nbin_a_max) :: dens_dry_a
 132     real(r8), intent(out), dimension(nbin_a_max) :: water_a_hyst
 133     real(r8), intent(out), dimension(nbin_a_max) :: aH2O_a
 134     real(r8), intent(out), dimension(nbin_a_max) :: gam_ratio
 135     real(r8), intent(out)                        :: uptkrate_h2so4  ! rate of h2so4 uptake by aerosols (1/s)
 136
 137     !Local Variables
 138     integer :: iprint_input, irepeat_mosaic
 139     integer :: mcall_print_aer
 140     integer, dimension(nbin_a_max) :: jphase
 141
 142     integer :: iaer !BALLI- remove this after debugging
 143
 144     real(r8) :: sigma_water,Kp_nh4cl
 145     real(r8) :: Kp_nh4no3,Kp_nh3
 146     real(r8) :: tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in
 147
 148     real(r8), dimension(nbin_a_max) :: mass_soluble_a
 149     real(r8), dimension(ngas_volatile) :: sat_soa,total_species
 150     real(r8), dimension(nrxn_aer_sl) :: Keq_sl
 151     real(r8), dimension(nrxn_aer_ll) :: Keq_ll
 152     real(r8), dimension(nrxn_aer_gl) :: Keq_gl
 153     real(r8), dimension(nrxn_aer_sg) :: Keq_sg
 154     real(r8), dimension(MDRH_T_NUM) :: MDRH_T
 155     real(r8), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
 156     real(r8), dimension(ngas_volatile,nbin_a_max) :: flux_s,flux_l
 157     real(r8), dimension(ngas_volatile,nbin_a_max) :: volatile_s
 158     real(r8), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
 159     real(r8), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
 160     real(r8), dimension(ngas_aerchtot,nbin_a_max) :: kg
 161     real(r8), dimension(nelectrolyte,nbin_a_max) :: activity,gam
 162     real(r8), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
 163     real(r8), dimension(Ncation,nbin_a_max) :: mc
 164     real(r8), dimension(Nanion,nbin_a_max) :: ma
 165     real(r8), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 166
 167
 168     call update_thermodynamic_constants(  aH2O,    T_K,                                & !intent-ins
 169          sat_soa,    aH2O_a,   log_gamZ,  Keq_sl,  sigma_water,  Kp_nh4cl,             & !intent-outs
 170          Kp_nh4no3,  Kp_nh3,   Keq_ll,    Keq_gl,  Keq_sg,       MDRH_T,               &
 171          molality0                                                                     )
 172
 173 ! rc_easter 2013-07-30 -
 174 ! the purpose of the irepeat loop was to provide more accurate cpu timings
 175 ! now that the cnn<-->gas,aer mapping is done earlier, you would have to
 176 !    save the gas,aer,num_a,... arrays then restore them for each repeat cycle
 177     do irepeat_mosaic = 1, 1
 178        mcall_print_aer = mcall_print_aer_in
 179        if (irepeat_mosaic > 1) mcall_print_aer = 0
 180
 181        call initialize_mosaic_variables(                                                & !intent-ins
 182             jaerosolstate, flux_s, flux_l, volatile_s, phi_volatile_s, phi_volatile_l,  & !intent-outs
 183             jphase, kg, electrolyte, activity, mc, mass_dry_a, mass_soluble_a,          &
 184             dens_dry_a, ma, gam, gam_ratio                                              )
 185
 186        mosaic_vars_aa%isteps_astem = 0
 187        mosaic_vars_aa%isteps_astem_max = 0
 188        mosaic_vars_aa%jastem_call = 0
 189        mosaic_vars_aa%jmesa_call = 0
 190        mosaic_vars_aa%jmesa_fail = 0
 191        mosaic_vars_aa%niter_mesa_max = 0
 192        mosaic_vars_aa%nmax_astem = nmax_astem
 193        mosaic_vars_aa%nmax_mesa = nmax_mesa
 194        mosaic_vars_aa%cumul_steps_astem = 0.0_r8
 195        mosaic_vars_aa%niter_mesa = 0.0_r8
 196        uptkrate_h2so4 = 0.0_r8
 197
 198        call overall_massbal_in( aer, gas, gas_netprod_otrproc, dtchem,                  & !intent-ins
 199             total_species, tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in,    & !intent-outs
 200             tot_ca_in )
 201
 202        call MOSAIC_dynamic_solver(      mcall_print_aer,     dtchem,                    & !intent-ins
 203             aH2O,           T_K,        RH_pc,               P_atm,                     &
 204             irepeat_mosaic, tot_cl_in,  sigmag_a,            kappa_nonelectro,          &
 205             jaerosolstate,  flux_s,     flux_l,              volatile_s,                & !intent-inouts
 206             phi_volatile_s, phi_volatile_l,                  jphase,           aer,     &
 207             kg,             gas,        gas_avg,             gas_netprod_otrproc,       &
 208             jhyst_leg,      electrolyte,                     activity,                  &
 209             mc,             sat_soa,    num_a,               Dp_dry_a,         Dp_wet_a,&
 210             mass_dry_a,     mass_soluble_a,                  dens_dry_a,       water_a, &
 211             gam,            log_gamZ,   gam_ratio,           Keq_ll,           Keq_gl,  &
 212             Keq_sg,         Keq_sl,     Kp_nh4cl,            Kp_nh4no3,        ma,      &
 213             sigma_water,    MDRH_T,     molality0,                                      &
 214             total_species,  aH2O_a,     uptkrate_h2so4,                                 &
 215             mosaic_vars_aa,                                                             &
 216             iprint_input,                                                               & !intent-outs
 217             mass_dry_a_bgn, dens_dry_a_bgn,                                             &
 218             water_a_hyst,   jaerosolstate_bgn                                           )
 219
 220        if (mosaic_vars_aa%f_mos_fail > 0) then
 221           return
 222        endif
 223
 224
 225        call overall_massbal_out( iprint_input, 0, mosaic_vars_aa%isteps_ASTEM, aer, gas, &
 226           tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in )
 227
 228     enddo
 229
 230
 231
 232
 233     return
 234   end subroutine mosaic_box_aerchemistry
 235
 236
 237
 238   !***********************************************************************
 239   ! subroutine to calculate just aerosol water
 240   !
 241   ! author: Rahul A. Zaveri
 242   ! update: May 2016
 243   !-----------------------------------------------------------------------
 244   subroutine MOSAIC_aerosol_water_only(      aH2O,         T_K,                  &!Intent-ins
 245        P_atm,                  RH_pc,        dtchem,                             &
 246        kappa_nonelectro,                                                         &
 247        jaerosolstate,          jhyst_leg,                                        &!Intent-inouts
 248        aer,                    num_a,        water_a,      gas,                  &
 249        Dp_dry_a,               Dp_wet_a,                                         &
 250        mosaic_vars_aa,                                                           &
 251        mass_dry_a,             dens_dry_a                                        )!Intent-outs
 252
 253     use module_data_mosaic_aero,  only:                                                 &
 254          nbin_a_max, ngas_aerchtot, ngas_volatile, nelectrolyte,                        &!Parameters
 255          Nanion, Ncation, naer, nbin_a, no_aerosol, jtotal,                             &
 256          jsalt_index, jsulf_poor, jsulf_rich,                                           &
 257          jhyst_lo, jhyst_up, mhyst_method, mhyst_uporlo_jhyst,                          &
 258          mhyst_uporlo_waterhyst, mhyst_force_lo, mhyst_force_up,                        &
 259          mSECTIONAL, mSIZE_FRAMEWORK, MDRH_T_NUM,                                       &
 260          nrxn_aer_gl, nrxn_aer_ll, nrxn_aer_sg, nrxn_aer_sl, nsalt,                     &
 261          zc, za, a_zsr, b_zsr, aw_min,                                                  &
 262          mw_electrolyte, partial_molar_vol,                                             &
 263          dens_aer_mac, mw_aer_mac, dens_comp_a, mw_comp_a, ref_index_a, MW_a, MW_c,     &
 264          density_max_allow, density_min_allow,                                          &
 265          rtol_mesa, nmax_astem, nmax_mesa, mosaic_vars_aa_type
 266
 267     use module_data_mosaic_asecthp, only: isize_of_ibin, itype_of_ibin,dcen_sect       ! TBD
 268
 269     use module_mosaic_ext,        only: aerosol_water_up, aerosol_phase_state, &
 270                                         calc_dry_n_wet_aerosol_props, conform_electrolytes
 271
 272     implicit none
 273
 274     !Intent-ins
 275     real(r8), intent(in) :: dtchem
 276     real(r8), intent(in) :: aH2O, T_K, RH_pc, P_atm
 277
 278     real(r8), intent(in), dimension(naer)       :: kappa_nonelectro
 279
 280     !Intent-inouts
 281     integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
 282     integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate
 283
 284     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 285     real(r8), intent(inout), dimension(nbin_a_max) :: Dp_dry_a, Dp_wet_a
 286     real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a
 287     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
 288     real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a
 289     real(r8), intent(inout), dimension(nbin_a_max) :: num_a
 290     real(r8), intent(inout), dimension(nbin_a_max) :: water_a
 291
 292     type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
 293
 294     !Intent-outs
 295
 296
 297
 298     !Local variables
 299     character(len=256) :: errmsg
 300     integer :: ibin, isize, itype, iv
 301     integer :: irepeat_mosaic
 302     integer, dimension(nbin_a_max) :: jphase
 303     integer, dimension(nbin_a_max) :: jaerosolstate_bgn
 304     integer, dimension(ngas_volatile,3,nbin_a_max) :: integrate
 305     integer, dimension(nsalt) :: jsalt_present
 306
 307     real(r8) :: Keq_nh4cl
 308     real(r8) :: Kp_nh3, Kp_nh4cl, Kp_nh4no3
 309     real(r8) :: sigma_water
 310     real(r8) :: tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in
 311     real(r8) :: XT
 312
 313     real(r8), dimension(MDRH_T_NUM) :: MDRH_T
 314     real(r8), dimension(Nanion ) :: na_Ma, xeq_a
 315     real(r8), dimension(Ncation) :: nc_Mc, xeq_c
 316     real(r8), dimension(Nanion, nbin_a_max) :: ma
 317     real(r8), dimension(Ncation,nbin_a_max) :: mc
 318     real(r8), dimension(nsalt) :: phi_salt_old
 319
 320     real(r8), dimension(nbin_a_max) :: aH2O_a
 321     real(r8), dimension(nbin_a_max) :: area_dry_a, area_wet_a
 322     real(r8), dimension(nbin_a_max) :: delta_hcl_max, delta_nh3_max, delta_hno3_max
 323     real(r8), dimension(nbin_a_max) :: dens_dry_a_bgn, dens_wet_a
 324     real(r8), dimension(nbin_a_max) :: dp_core_a
 325     real(r8), dimension(nbin_a_max) :: gam_ratio
 326     real(r8), dimension(nbin_a_max) :: growth_factor
 327     real(r8), dimension(nbin_a_max) :: mass_dry_a_bgn, mass_soluble_a, mass_wet_a
 328     real(r8), dimension(nbin_a_max) :: MDRH
 329     real(r8), dimension(nbin_a_max) :: sigma_soln
 330     real(r8), dimension(nbin_a_max) :: vol_dry_a, vol_wet_a
 331     real(r8), dimension(nbin_a_max) :: water_a_hyst, water_a_up
 332
 333     real(r8), dimension(ngas_aerchtot) :: gas_netprod_otrproc
 334     real(r8), dimension(ngas_volatile) :: sfc_a
 335     real(r8), dimension(ngas_volatile) :: sat_soa, total_species
 336     real(r8), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l, phi_volatile_s
 337     real(r8), dimension(ngas_volatile,nbin_a_max) :: volatile_s
 338     real(r8), dimension(ngas_volatile,nbin_a_max) :: flux_s, flux_l
 339     real(r8), dimension(ngas_aerchtot,nbin_a_max) :: kel
 340     real(r8), dimension(ngas_aerchtot,nbin_a_max) :: kg
 341
 342     real(r8), dimension(nelectrolyte,nbin_a_max) :: activity, gam
 343     real(r8), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
 344     real(r8), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 345     real(r8), dimension(nelectrolyte,3,nbin_a_max) :: epercent
 346     real(r8), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
 347
 348     real(r8), dimension(nrxn_aer_sl) :: Keq_sl
 349     real(r8), dimension(nrxn_aer_ll) :: Keq_ll
 350     real(r8), dimension(nrxn_aer_gl) :: Keq_gl
 351     real(r8), dimension(nrxn_aer_sg) :: Keq_sg
 352
 353     complex, dimension(nbin_a_max) :: ri_shell_a,ri_avg_a,ri_core_a
 354
 355
 356     gas_netprod_otrproc(1:ngas_volatile) = 0.0_r8
 357
 358     call update_thermodynamic_constants(  aH2O,    T_K,                                & !intent-ins
 359          sat_soa,    aH2O_a,   log_gamZ,  Keq_sl,  sigma_water,  Kp_nh4cl,             & !intent-outs
 360          Kp_nh4no3,  Kp_nh3,   Keq_ll,    Keq_gl,  Keq_sg,       MDRH_T,               &
 361          molality0                                                                     )
 362
 363     irepeat_mosaic = 1
 364
 365     call initialize_mosaic_variables(                                                & !intent-ins
 366          jaerosolstate, flux_s, flux_l, volatile_s, phi_volatile_s, phi_volatile_l,  & !intent-outs
 367          jphase, kg, electrolyte, activity, mc, mass_dry_a, mass_soluble_a,          &
 368          dens_dry_a, ma, gam, gam_ratio                                              )
 369
 370     mosaic_vars_aa%isteps_astem = 0
 371     mosaic_vars_aa%isteps_astem_max = 0
 372     mosaic_vars_aa%jastem_call = 0
 373     mosaic_vars_aa%jmesa_call = 0
 374     mosaic_vars_aa%jmesa_fail = 0
 375     mosaic_vars_aa%niter_mesa_max = 0
 376     mosaic_vars_aa%nmax_astem = nmax_astem
 377     mosaic_vars_aa%nmax_mesa = nmax_mesa
 378     mosaic_vars_aa%cumul_steps_astem = 0.0_r8
 379     mosaic_vars_aa%niter_mesa = 0.0_r8
 380
 381     call overall_massbal_in( aer, gas, gas_netprod_otrproc, dtchem,                  & !intent-ins
 382          total_species, tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in,    & !intent-outs
 383          tot_ca_in )
 384
 385
 386     vol_dry_a = 0.0_r8!*BALLI- ASK dick, if we dont initialize it here the code blows up. In conform_aerosol_number, vol_dry_a do not get any value as num_a(ibin)>0.0
 387
 388     do ibin = 1, nbin_a
 389        call check_aerosol_mass(ibin, jaerosolstate,jphase,aer,num_a, mass_dry_a)
 390        jaerosolstate_bgn(ibin) = jaerosolstate(ibin)
 391
 392        if(jaerosolstate(ibin) .ne. no_aerosol) then!goto 500
 393
 394           call conform_electrolytes(jtotal,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)        ! conforms aer(jtotal) to a valid aerosol
 395           call check_aerosol_mass(ibin,jaerosolstate,jphase,aer,num_a, mass_dry_a) ! check mass again after conform_electrolytes
 396
 397           jaerosolstate_bgn(ibin) = jaerosolstate(ibin)
 398           if(jaerosolstate(ibin) .ne. no_aerosol)then !goto 500    ! ignore this bin
 399
 400              call conform_aerosol_number(ibin,jaerosolstate,aer,num_a,vol_dry_a,Dp_dry_a) ! adjusts number conc so that it conforms with bin mass and diameter
 401
 402              ! when mhyst_method = mhyst_uporlo_waterhyst,
 403              ! initialize water_a_hyst at first time step using the user-input jhyst_leg
 404              !BSINGH - 05/28/2013(RCE updates - if cond structure has been modified)
 405              if (mosaic_vars_aa%it_mosaic == 1) then
 406                 if (mhyst_method == mhyst_uporlo_waterhyst) then
 407                    if(jhyst_leg(ibin) == jhyst_lo)then
 408                       water_a_hyst(ibin) = 0.0
 409                    else
 410                       water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)        ! at 60% RH
 411                       water_a_hyst(ibin) = water_a_up(ibin)
 412                    endif
 413                 else if (mhyst_method == mhyst_force_lo) then
 414                    jhyst_leg(ibin) = jhyst_lo
 415                    water_a_hyst(ibin) = 0.0
 416                 else if (mhyst_method == mhyst_force_up) then
 417                    jhyst_leg(ibin)    = jhyst_up
 418                    water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH
 419                    water_a_hyst(ibin) = water_a_up(ibin)
 420                 end if
 421              end if
 422              !BSINGH - 05/28/2013(RCE updates)
 423           endif
 424        endif
 425        if (irepeat_mosaic == 1) then
 426           mass_dry_a_bgn(ibin) = mass_dry_a(ibin)
 427           if ( (jaerosolstate(ibin) .eq. no_aerosol) .or.   &
 428                (min(mass_dry_a(ibin),vol_dry_a(ibin)) .le. 1.0e-35) ) then
 429              call calc_aerosol_dry_density( ibin,aer,dens_dry_a)
 430              dens_dry_a_bgn(ibin) = dens_dry_a(ibin)
 431           else
 432              dens_dry_a_bgn(ibin) = mass_dry_a(ibin)/vol_dry_a(ibin)
 433           end if
 434           dens_dry_a_bgn(ibin) = max( density_min_allow, &
 435                min( density_max_allow, dens_dry_a_bgn(ibin) ) )
 436        end if
 437
 438        if (jaerosolstate(ibin) .eq. no_aerosol) then
 439           if (msize_framework == msectional) then
 440              isize = isize_of_ibin(ibin)
 441              itype = itype_of_ibin(ibin)
 442              Dp_dry_a(ibin) = dcen_sect(isize,itype)
 443              Dp_wet_a(ibin) = Dp_dry_a(ibin)
 444           end if
 445        end if
 446
 447     enddo
 448
 449
 450
 451   ! compute aerosol phase state
 452   do ibin = 1, nbin_a
 453
 454      if(jaerosolstate(ibin) .ne. no_aerosol)then
 455         call aerosol_phase_state( ibin, jaerosolstate, jphase,  &
 456              aer, jhyst_leg, electrolyte, epercent, kel, activity, mc, num_a, mass_wet_a, &
 457              mass_dry_a, mass_soluble_a, vol_dry_a, vol_wet_a, water_a, water_a_hyst,  &
 458              water_a_up, aH2O_a, aH2O, ma, gam, & !BALLI
 459              log_gamZ, zc, za, gam_ratio, xeq_a, na_Ma, nc_Mc, xeq_c,              & ! RAZ deleted a_zsr
 460              mw_electrolyte, partial_molar_vol, sigma_soln, T_K, RH_pc, mw_aer_mac,    &
 461              dens_aer_mac, sigma_water, Keq_ll, Keq_sl, MW_a, MW_c, growth_factor, MDRH, &
 462              MDRH_T, molality0, rtol_mesa, jsalt_present, jsalt_index, jsulf_poor,     &
 463              jsulf_rich, phi_salt_old,                                      &
 464              kappa_nonelectro, mosaic_vars_aa )
 465
 466         call calc_dry_n_wet_aerosol_props(                                &
 467            ibin, jaerosolstate, aer, electrolyte, water_a, num_a,         &  ! input
 468            dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &  ! input
 469            Dp_dry_a, Dp_wet_a, dp_core_a,                                 &  ! output
 470            area_dry_a, area_wet_a, mass_dry_a, mass_wet_a,                &  ! output
 471            vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a,                  &  ! output
 472            ri_shell_a, ri_core_a, ri_avg_a                                )  ! output
 473      endif
 474   enddo
 475
 476
 477     do ibin = 1, nbin_a
 478        if(jaerosolstate(ibin).ne.no_aerosol) then
 479
 480 !         if (mhyst_method == mhyst_uporlo_jhyst) then ! rce 2017.10.21
 481           if (mhyst_method == mhyst_uporlo_jhyst .or. & ! rce 2017.10.21
 482               mhyst_method == mhyst_uporlo_waterhyst) then
 483              ! rce 2017.10.21 - for both mhyst_uporlo_... use jhyst_leg (from astem) to set water_a_hyst
 484              if(jhyst_leg(ibin) == jhyst_lo)then
 485                 water_a_hyst(ibin) = 0.0
 486              else
 487                 water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH
 488                 water_a_hyst(ibin) = water_a_up(ibin)
 489              endif
 490 !         elseif (mhyst_method == mhyst_uporlo_waterhyst) then ! rce 2017.10.21
 491 !            water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)      ! at 60% RH
 492 !            if (water_a_hyst(ibin) <= 0.5*water_a_up(ibin)) then
 493 !               jhyst_leg(ibin) = jhyst_lo
 494 !               water_a_hyst(ibin) = 0.0
 495 !            else
 496 !               jhyst_leg(ibin) = jhyst_up
 497 !               water_a_hyst(ibin) = water_a_up(ibin)
 498 !            endif
 499 !            !BSINGH - 05/28/2013(RCE updates)
 500           else if (mhyst_method == mhyst_force_lo) then
 501              jhyst_leg(ibin) = jhyst_lo
 502              water_a_hyst(ibin) = 0.0
 503           else if (mhyst_method == mhyst_force_up) then
 504              jhyst_leg(ibin) = jhyst_up
 505              water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH ! rce 2017.10.21
 506              water_a_hyst(ibin) = water_a_up(ibin)
 507              !BSINGH - 05/28/2013(RCE updates ENDS)
 508           else
 509              write(errmsg,*) '*** MOSAIC_aerosol_water - bad mhyst_method =', mhyst_method!BSINGH - 05/28/2013(RCE updates)
 510              call wrf_error_fatal(trim(adjustl(errmsg)))
 511           endif
 512
 513        endif
 514        if ( (jaerosolstate(ibin) .eq. no_aerosol) .or.   &
 515             (min(mass_dry_a(ibin),vol_dry_a(ibin)) .le. 1.0e-35) ) then
 516           call calc_aerosol_dry_density( ibin,aer,dens_dry_a)
 517        end if
 518        dens_dry_a(ibin) = max( density_min_allow, &
 519             min( density_max_allow, dens_dry_a(ibin) ) )
 520
 521     enddo
 522
 523     return
 524   end subroutine MOSAIC_aerosol_water_only
 525
 526
 527
 528   !***********************************************************************
 529   ! interface to dynamic gas-particle exchange solver
 530   !
 531   ! author: Rahul A. Zaveri
 532   ! update: jan 2005
 533   !-----------------------------------------------------------------------
 534
 535   subroutine MOSAIC_dynamic_solver( mcall_print_aer,    dtchem,                    & !intent-ins
 536        aH2O,           T_K,        RH_pc,               P_atm,                     &
 537        irepeat_mosaic, tot_cl_in,  sigmag_a,                                       &
 538        kappa_nonelectro,                                                           &
 539        jaerosolstate,  flux_s,     flux_l,              volatile_s,                & !intent-inouts
 540        phi_volatile_s, phi_volatile_l,                  jphase,           aer,     &
 541        kg,             gas,        gas_avg,             gas_netprod_otrproc,       &
 542        jhyst_leg,      electrolyte,                     activity,                  &
 543        mc,             sat_soa,    num_a,               Dp_dry_a,         Dp_wet_a,&
 544        mass_dry_a,     mass_soluble_a,                  dens_dry_a,       water_a, &
 545        gam,            log_gamZ,   gam_ratio,           Keq_ll,           Keq_gl,  &
 546        Keq_sg,         Keq_sl,     Kp_nh4cl,            Kp_nh4no3,        ma,      &
 547        sigma_water,    MDRH_T,     molality0,                                      &
 548        total_species,  aH2O_a,     uptkrate_h2so4,                                 &
 549        mosaic_vars_aa,                                                             &
 550        iprint_input,                                                               & !intent-outs
 551        mass_dry_a_bgn, dens_dry_a_bgn,                                             &
 552        water_a_hyst,   jaerosolstate_bgn                                           )
 553
 554     use module_data_mosaic_aero,  only: nbin_a_max, ngas_aerchtot, ngas_volatile,        &
 555          nelectrolyte,                                                                   &!Parameters
 556          Ncation, naer, no_aerosol, jtotal, mhyst_uporlo_waterhyst, jhyst_lo,            &!Parameters
 557          density_max_allow, density_min_allow, mSECTIONAL, mON, mASTEM, mLSODE,          &!Parameters
 558          mhyst_uporlo_jhyst, jhyst_up, Nanion, nrxn_aer_gl, nrxn_aer_ll,                &
 559          nrxn_aer_sg, nrxn_aer_sl, nsalt, MDRH_T_NUM,  mhyst_force_lo,  mhyst_force_up,  &
 560          nbin_a, mSIZE_FRAMEWORK, mGAS_AER_XFER, mDYNAMIC_SOLVER, mhyst_method,         &
 561          zc, za, a_zsr, mw_electrolyte, partial_molar_vol, dens_aer_mac,      &
 562          mw_aer_mac,  dens_comp_a, mw_comp_a, ref_index_a, MW_a, MW_c, rtol_mesa,         &
 563          jsalt_index, jsulf_poor, jsulf_rich,                                          &
 564          iso4_a,                                                                       & !balli for debug only remove it
 565          mosaic_vars_aa_type
 566
 567     use module_data_mosaic_asecthp, only: isize_of_ibin,itype_of_ibin,dcen_sect       ! TBD
 568
 569     use module_mosaic_astem,      only: ASTEM
 570
 571     use module_mosaic_ext,        only: aerosol_water_up,calc_dry_n_wet_aerosol_props,&
 572                                         conform_electrolytes, dumpxx
 573 !   use module_print_aer,         only: print_aer
 574     use module_mosaic_lsode,      only: mosaic_lsode
 575
 576     implicit none
 577
 578     !Intent-ins
 579     integer, intent(in) :: mcall_print_aer
 580     integer, intent(in) :: irepeat_mosaic
 581
 582     real(r8), intent(in) :: dtchem
 583     real(r8), intent(in) :: aH2O, T_K, RH_pc, P_atm
 584
 585     real(r8), intent(in), dimension(nbin_a_max) :: sigmag_a
 586     real(r8), intent(in), dimension(naer)       :: kappa_nonelectro
 587
 588     !Intent-inouts
 589     real(r8), intent(inout) :: Kp_nh4cl
 590     real(r8), intent(inout) :: Kp_nh4no3,sigma_water
 591     real(r8), intent(inout) :: tot_cl_in
 592
 593     real(r8), intent(inout), dimension(MDRH_T_NUM) :: MDRH_T
 594
 595     integer, intent(inout), dimension(nbin_a_max) :: jhyst_leg
 596     integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase
 597
 598     real(r8), intent(inout), dimension(nbin_a_max) :: num_a, Dp_dry_a
 599     real(r8), intent(inout), dimension(nbin_a_max) :: Dp_wet_a, gam_ratio
 600     real(r8), intent(inout), dimension(nbin_a_max) :: aH2O_a
 601
 602     real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a
 603     real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a
 604     real(r8), intent(inout), dimension(nbin_a_max) :: mass_soluble_a
 605     real(r8), intent(inout), dimension(nbin_a_max) :: water_a
 606     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
 607     real(r8), intent(inout), dimension(ngas_volatile) :: sat_soa, total_species
 608     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas_avg  ! average gas conc. over dtchem time step (nmol/m3)
 609     real(r8), intent(in),    dimension(ngas_aerchtot) :: gas_netprod_otrproc
 610               ! gas_netprod_otrproc = gas net production rate from other processes
 611               !    such as gas-phase chemistry and emissions (nmol/m3/s)
 612               ! this allows the condensation (gasaerexch) routine to apply production and condensation loss
 613               !    together, which is more accurate numerically
 614               ! NOTE - must be >= zero, as numerical method can fail when it is negative
 615               ! NOTE - currently for mosaic, only the value for h2so4 can be non-zero
 616
 617     real(r8), intent(inout), dimension(nrxn_aer_ll) :: Keq_ll
 618     real(r8), intent(inout), dimension(nrxn_aer_gl) :: Keq_gl
 619     real(r8), intent(inout), dimension(nrxn_aer_sg) :: Keq_sg
 620     real(r8), intent(inout), dimension(nrxn_aer_sl) :: Keq_sl
 621
 622     real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
 623
 624     real(r8), intent(inout), dimension(Ncation,nbin_a_max) :: mc
 625     real(r8), intent(inout), dimension(Nanion,nbin_a_max) :: ma
 626
 627     real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: flux_s,flux_l
 628     real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: volatile_s
 629     real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
 630     real(r8), intent(inout), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
 631     real(r8), intent(inout), dimension(ngas_aerchtot,nbin_a_max) :: kg
 632     real(r8), intent(inout), dimension(nelectrolyte,nbin_a_max) :: activity,gam
 633     real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
 634
 635     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 636     real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 637     real(r8), intent(inout) :: uptkrate_h2so4  ! rate of h2so4 uptake by aerosols (1/s)
 638
 639     type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
 640
 641     !Intent-outs
 642     integer, intent(out) :: iprint_input
 643     integer, intent(out), dimension(nbin_a_max) :: jaerosolstate_bgn
 644
 645     real(r8), intent(out), dimension(nbin_a_max) :: water_a_hyst
 646     real(r8), intent(out), dimension(nbin_a_max) :: mass_dry_a_bgn,dens_dry_a_bgn
 647
 648     !Local variables
 649     character(len=256) :: errmsg
 650     integer ibin, isize, itype, iv
 651
 652     real(r8) :: XT
 653
 654
 655     real(r8), dimension(nbin_a_max) :: area_dry_a,water_a_up
 656     real(r8), dimension(nbin_a_max) :: area_wet_a,mass_wet_a,vol_wet_a,dens_wet_a
 657     real(r8), dimension(nbin_a_max) :: vol_dry_a
 658     real(r8), dimension(nbin_a_max) :: dp_core_a
 659
 660     real(r8), dimension(nelectrolyte,3,nbin_a_max) :: epercent
 661
 662     complex, dimension(nbin_a_max) :: ri_shell_a,ri_avg_a,ri_core_a
 663
 664
 665     vol_dry_a = 0.0_r8!*BALLI- ASK dick, if we dont initialize it here the code blows up. In conform_aerosol_number, vol_dry_a do not get any value as num_a(ibin)>0.0
 666
 667     !BSINGH - Initialize counters
 668     mosaic_vars_aa%jASTEM_fail = 0
 669     mosaic_vars_aa%jASTEM_call       = 0
 670     mosaic_vars_aa%isteps_ASTEM      = 0
 671     mosaic_vars_aa%isteps_ASTEM_max  = 0
 672     mosaic_vars_aa%niter_MESA        = 0.0_r8
 673     mosaic_vars_aa%cumul_steps_ASTEM = 0.0_r8
 674
 675
 676     call dumpxx( 'aa', dtchem, t_k, p_atm, ah2o, &
 677          jaerosolstate, jphase, jhyst_leg, &
 678          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 679          mosaic_vars_aa )
 680
 681
 682     do ibin = 1, nbin_a
 683        call check_aerosol_mass(ibin, jaerosolstate,jphase,aer,num_a, mass_dry_a)
 684        jaerosolstate_bgn(ibin) = jaerosolstate(ibin)
 685
 686        if(jaerosolstate(ibin) .ne. no_aerosol) then!goto 500
 687
 688           !call conform_aerosol_number(ibin,jaerosolstate,aer,num_a,vol_dry_a, Dp_dry_a)     ! adjusts number conc so that it conforms with bin mass and diameter
 689
 690           call conform_electrolytes(jtotal,ibin,XT,aer,gas,electrolyte,total_species,tot_cl_in)        ! conforms aer(jtotal) to a valid aerosol
 691           call check_aerosol_mass(ibin,jaerosolstate,jphase,aer,num_a, mass_dry_a) ! check mass again after conform_electrolytes
 692
 693           jaerosolstate_bgn(ibin) = jaerosolstate(ibin)
 694           if(jaerosolstate(ibin) .ne. no_aerosol)then !goto 500    ! ignore this bin
 695
 696              ! *** moved "call conform_aerosol_number" here instead of above by RAZ
 697              call conform_aerosol_number(ibin,jaerosolstate,aer,num_a,vol_dry_a,Dp_dry_a) ! adjusts number conc so that it conforms with bin mass and diameter
 698
 699              ! when mhyst_method = mhyst_uporlo_waterhyst,
 700              ! initialize water_a_hyst at first time step using the user-input jhyst_leg
 701              !BSINGH - 05/28/2013(RCE updates - if cond structure has been modified)
 702              if (mosaic_vars_aa%it_mosaic == 1) then
 703                 if (mhyst_method == mhyst_uporlo_waterhyst) then
 704                    if(jhyst_leg(ibin) == jhyst_lo)then
 705                       water_a_hyst(ibin) = 0.0
 706                    else
 707                       water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)        ! at 60% RH
 708                       water_a_hyst(ibin) = water_a_up(ibin)
 709                    endif
 710                 else if (mhyst_method == mhyst_force_lo) then
 711                    jhyst_leg(ibin) = jhyst_lo
 712                    water_a_hyst(ibin) = 0.0
 713                 else if (mhyst_method == mhyst_force_up) then
 714                    jhyst_leg(ibin)    = jhyst_up
 715                    water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH
 716                    water_a_hyst(ibin) = water_a_up(ibin)
 717                 end if
 718              end if
 719              !BSINGH - 05/28/2013(RCE updates)
 720           endif
 721        endif
 722        if (irepeat_mosaic == 1) then
 723           mass_dry_a_bgn(ibin) = mass_dry_a(ibin)
 724           if ( (jaerosolstate(ibin) .eq. no_aerosol) .or.   &
 725                (min(mass_dry_a(ibin),vol_dry_a(ibin)) .le. 1.0e-35) ) then
 726              call calc_aerosol_dry_density( ibin,aer,dens_dry_a)
 727              dens_dry_a_bgn(ibin) = dens_dry_a(ibin)
 728           else
 729              dens_dry_a_bgn(ibin) = mass_dry_a(ibin)/vol_dry_a(ibin)
 730           end if
 731           dens_dry_a_bgn(ibin) = max( density_min_allow, &
 732                min( density_max_allow, dens_dry_a_bgn(ibin) ) )
 733        end if
 734
 735        if (jaerosolstate(ibin) .eq. no_aerosol) then
 736           if (msize_framework == msectional) then
 737              isize = isize_of_ibin(ibin)
 738              itype = itype_of_ibin(ibin)
 739              Dp_dry_a(ibin) = dcen_sect(isize,itype)
 740              Dp_wet_a(ibin) = Dp_dry_a(ibin)
 741           end if
 742        end if
 743
 744     enddo
 745
 746     !cc        call save_pregrow_props !3D
 747     !cc        call specialoutaa( iclm_aer, jclm_aer, kclm_aer, 77, ! 3D
 748     !cc     &          'after_conform' )
 749     !
 750     !-------------------------------------
 751     ! do dynamic gas-aerosol mass transfer for dtchem [s]
 752
 753     if(mGAS_AER_XFER .eq. mON)then
 754        !        call wall_loss(dtchem)
 755
 756        if(mDYNAMIC_SOLVER .eq. mASTEM)then
 757           call ASTEM( mcall_print_aer,          dtchem,           &!intent-ins
 758                sigmag_a,  aH2O,     T_K,         RH_pc,        P_atm,                        &
 759                kappa_nonelectro,                                                             &
 760                jaerosolstate, flux_s,            flux_l,       volatile_s, iprint_input,     &!intent -inout
 761                phi_volatile_s,phi_volatile_l,    jphase,       aer,       kg,       gas,     &
 762                gas_avg,       gas_netprod_otrproc,                                           &
 763                jhyst_leg,     electrolyte,       epercent,     activity,  mc,       sat_soa, &
 764                num_a,         Dp_dry_a,          Dp_wet_a,     dp_core_a, mass_dry_a,        &
 765                mass_soluble_a,vol_dry_a,         dens_dry_a,   water_a,   water_a_hyst,      &
 766                water_a_up,    aH2O_a,            total_species,tot_cl_in, ma,       gam,     &
 767                log_gamZ,      gam_ratio,         Keq_ll,       Keq_gl,    Keq_sg,   Kp_nh4cl,&
 768                Kp_nh4no3,     sigma_water,       Keq_sl,       MDRH_T,    molality0,         &
 769                uptkrate_h2so4,                   mosaic_vars_aa,                             &
 770                area_dry_a,    area_wet_a,        mass_wet_a,vol_wet_a,                       &!intent-out
 771                dens_wet_a,    ri_shell_a,        ri_avg_a,     ri_core_a                     )
 772
 773           if (mosaic_vars_aa%f_mos_fail > 0) then
 774              return
 775           endif
 776
 777           !call ASTEM( mcall_print_aer,                                             &
 778           !     iprint_input,jASTEM_call,dtchem,jaerosolstate,isteps_ASTEM,         &
 779           !     iter_MESA,jMESA_call,flux_s,flux_l,volatile_s,phi_volatile_s,       &
 780           !     phi_volatile_l,jphase,aer,kg,gas,jhyst_leg,electrolyte,epercent,    &
 781           !     activity,mc,sat_soa,delta_nh3_max,delta_hno3_max,delta_hcl_max,     &
 782           !     jASTEM_fail,jMESA_fail,isteps_ASTEM_max,nmax_ASTEM,cumul_steps_ASTEM,num_a,    &
 783           !     Dp_dry_a,Dp_wet_a,dp_core_a,area_dry_a,area_wet_a,mass_wet_a,       &
 784           !     mass_dry_a,mass_soluble_a,vol_dry_a,vol_wet_a,dens_dry_a,dens_wet_a,&
 785           !     sigmag_a,water_a,water_a_hyst,water_a_up,aH2O_a,total_species,      &
 786           !     tot_cl_in,aH2O,                                                     &
 787           !     niter_MESA_max,niter_MESA,ma,gam,log_gamZ,zc,za,gam_ratio,xeq_a,    &
 788           !     na_Ma,nc_Mc,xeq_c,a_zsr,mw_electrolyte,partial_molar_vol,Keq_ll,    &
 789           !     Keq_gl,Keq_sg,Kp_nh4cl,Kp_nh4no3,Keq_nh4cl,sigma_soln,T_K,RH_pc,    &
 790           !     mw_aer_mac,dens_aer_mac,sigma_water,Keq_sl,MW_a,MW_c,ri_shell_a,    &
 791           !     dens_comp_a,mw_comp_a,ref_index_a,ri_avg_a,ri_core_a,P_atm,         &
 792           !     growth_factor,MDRH,MDRH_T,molality0,rtol_mesa,jsalt_present,        &
 793           !     jsalt_index,jsulf_poor,jsulf_rich,Nmax_mesa, phi_salt_old,          &
 794           !     zero_water_flag                                                     )
 795        elseif(mDYNAMIC_SOLVER .eq. mLSODE)then
 796
 797           call MOSAIC_LSODE(dtchem)
 798
 799        endif
 800
 801     endif
 802
 803     !-------------------------------------
 804
 805     ! grows or shrinks size depending on mass increase or decrease
 806
 807     do ibin = 1, nbin_a
 808        if(jaerosolstate(ibin) .ne. no_aerosol)then
 809           call conform_aerosol_size( ibin,jaerosolstate,aer,num_a,       Dp_dry_a,  &
 810                vol_dry_a,mw_aer_mac,dens_aer_mac, mosaic_vars_aa )    ! BOX
 811           if (mosaic_vars_aa%f_mos_fail > 0) then
 812              return
 813           endif
 814        endif
 815     enddo
 816
 817
 818     do ibin = 1, nbin_a
 819        if(jaerosolstate(ibin).ne.no_aerosol) then
 820
 821 !         if (mhyst_method == mhyst_uporlo_jhyst) then ! rce 2017.10.21
 822           if (mhyst_method == mhyst_uporlo_jhyst .or. & ! rce 2017.10.21
 823               mhyst_method == mhyst_uporlo_waterhyst) then
 824              ! rce 2017.10.21 - for both mhyst_uporlo_... use jhyst_leg (from astem) to set water_a_hyst
 825              if(jhyst_leg(ibin) == jhyst_lo)then
 826                 water_a_hyst(ibin) = 0.0
 827              else
 828                 water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH
 829                 water_a_hyst(ibin) = water_a_up(ibin)
 830              endif
 831 !         elseif (mhyst_method == mhyst_uporlo_waterhyst) then ! rce 2017.10.21
 832 !            water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)      ! at 60% RH
 833 !            if (water_a_hyst(ibin) <= 0.5*water_a_up(ibin)) then
 834 !               jhyst_leg(ibin) = jhyst_lo
 835 !               water_a_hyst(ibin) = 0.0
 836 !            else
 837 !               jhyst_leg(ibin) = jhyst_up
 838 !               water_a_hyst(ibin) = water_a_up(ibin)
 839 !            endif
 840 !            !BSINGH - 05/28/2013(RCE updates)
 841           else if (mhyst_method == mhyst_force_lo) then
 842              jhyst_leg(ibin) = jhyst_lo
 843              water_a_hyst(ibin) = 0.0
 844           else if (mhyst_method == mhyst_force_up) then
 845              jhyst_leg(ibin) = jhyst_up
 846              water_a_up(ibin)   = aerosol_water_up(ibin,electrolyte,aer,kappa_nonelectro,a_zsr)   ! at 60% RH ! rce 2017.10.21
 847              water_a_hyst(ibin) = water_a_up(ibin)
 848              !BSINGH - 05/28/2013(RCE updates ENDS)
 849           else
 850              write(errmsg,*) '*** MOSAIC_dynamic_solver - bad mhyst_method =', mhyst_method!BSINGH - 05/28/2013(RCE updates)
 851              call wrf_error_fatal(trim(adjustl(errmsg)))
 852           endif
 853
 854           ! compute final mass and density
 855           call calc_dry_n_wet_aerosol_props(                                &
 856                ibin, jaerosolstate, aer, electrolyte, water_a, num_a,         &  ! input
 857                dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &  ! input
 858                Dp_dry_a, Dp_wet_a, dp_core_a,                                 &  ! output
 859                area_dry_a, area_wet_a, mass_dry_a, mass_wet_a,                &  ! output
 860                vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a,                  &  ! output
 861                ri_shell_a, ri_core_a, ri_avg_a                                )  ! output
 862
 863        endif
 864        if ( (jaerosolstate(ibin) .eq. no_aerosol) .or.   &
 865             (min(mass_dry_a(ibin),vol_dry_a(ibin)) .le. 1.0e-35) ) then
 866           call calc_aerosol_dry_density( ibin,aer,dens_dry_a)
 867        end if
 868        dens_dry_a(ibin) = max( density_min_allow, &
 869             min( density_max_allow, dens_dry_a(ibin) ) )
 870
 871     enddo
 872
 873     if (mcall_print_aer == 1 .or. mcall_print_aer == 2) then
 874        !call print_aer(1,jaerosolstate,isteps_ASTEM,iter_MESA,aer,gas,electrolyte,  &
 875        !     mc,num_a,Dp_dry_a,Dp_wet_a,area_dry_a,area_wet_a,mass_wet_a,mass_dry_a,&
 876        !     water_a)
 877     end if
 878
 879
 880     call dumpxx( 'zz', dtchem, t_k, p_atm, ah2o, &
 881          jaerosolstate, jphase, jhyst_leg, &
 882          aer, gas, num_a, water_a, water_a_hyst, dp_dry_a, &
 883          mosaic_vars_aa )
 884
 885
 886     return
 887   end subroutine MOSAIC_dynamic_solver
 888
 889
 890
 891   !***********************************************************************
 892   ! applies first-order wall loss to number and mass
 893   !
 894   ! author: Rahul A. Zaveri
 895   ! update: jun 2003
 896   !-----------------------------------------------------------------------
 897   subroutine wall_loss(dtchem,aer,num_a)
 898     use module_data_mosaic_aero, only: nbin_a_max,naer,jtotal,jsolid,jliquid,   & !Parameters
 899          nbin_a !Input
 900
 901     implicit none
 902     ! subr arguments
 903     real(r8), intent(in) :: dtchem
 904     real(r8), intent(inout), dimension(nbin_a_max) :: num_a
 905     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
 906     ! local variables
 907     integer  :: iaer, ibin
 908     real(r8) :: kwall
 909
 910
 911     kwall =  5.55e-5  ! 1/s
 912
 913     do ibin = 1, nbin_a
 914
 915        do iaer = 1, naer
 916           aer(iaer,jtotal,ibin)  = aer(iaer,jtotal,ibin)*exp(-kwall*dtchem)
 917           aer(iaer,jsolid,ibin)  = aer(iaer,jsolid,ibin)*exp(-kwall*dtchem)
 918           aer(iaer,jliquid,ibin) = aer(iaer,jliquid,ibin)*exp(-kwall*dtchem)
 919        enddo
 920
 921        num_a(ibin) = num_a(ibin)*exp(-kwall*dtchem)
 922
 923     enddo
 924
 925
 926     return
 927   end subroutine wall_loss
 928
 929
 930
 931   !***********************************************************************
 932   ! intializes all the MOSAIC variables to zero or their default values.
 933   !
 934   ! author: Rahul A. Zaveri
 935   ! update: jun 2003
 936   !-----------------------------------------------------------------------
 937   subroutine initialize_mosaic_variables(                                          & !intent-ins
 938        jaerosolstate, flux_s, flux_l, volatile_s, phi_volatile_s, phi_volatile_l,  & !intent-inouts
 939        jphase, kg, electrolyte, activity, mc, mass_dry_a, mass_soluble_a,          &
 940        dens_dry_a, ma, gam, gam_ratio                                              )
 941
 942     use module_data_mosaic_aero, only: nbin_a_max,nbin_a,naer,                     &
 943          ngas_aerchtot, ngas_volatile,                                             &!Parameters
 944          nelectrolyte,Ncation,ngas_ioa,jtotal,jsolid,jliquid,nanion
 945
 946
 947     implicit none
 948
 949     !Subroutine Arguments
 950     integer, intent(out), dimension(nbin_a_max) :: jaerosolstate,jphase
 951
 952     real(r8), intent(out), dimension(nbin_a_max) :: mass_dry_a,gam_ratio
 953     real(r8), intent(out), dimension(nbin_a_max) :: mass_soluble_a,dens_dry_a
 954
 955     real(r8), intent(out), dimension(ngas_aerchtot,nbin_a_max) :: kg
 956     real(r8), intent(out), dimension(ngas_volatile,nbin_a_max) :: flux_s
 957     real(r8), intent(out), dimension(ngas_volatile,nbin_a_max) :: flux_l
 958     real(r8), intent(out), dimension(ngas_volatile,nbin_a_max) :: volatile_s
 959     real(r8), intent(out), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_s
 960     real(r8), intent(out), dimension(ngas_volatile,nbin_a_max) :: phi_volatile_l
 961     real(r8), intent(out), dimension(nelectrolyte,nbin_a_max)  :: activity,gam
 962     real(r8), intent(out), dimension(Ncation,nbin_a_max)       :: mc
 963     real(r8), intent(out), dimension(Nanion,nbin_a_max)        :: ma
 964
 965     real(r8), intent(out), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
 966
 967     ! local variables
 968     integer iaer, ibin, iv, ja, jc, je
 969
 970     phi_volatile_l(:,:) = 0.0_r8 !BALLI** Ask dick about this initialization
 971
 972     ! initialize to zero
 973     do ibin = 1, nbin_a
 974
 975        mass_dry_a(ibin)     = 0.0
 976        mass_soluble_a(ibin) = 0.0
 977        dens_dry_a(ibin)     =-1.0
 978
 979        do je = 1, nelectrolyte
 980           electrolyte(je,jtotal,ibin)  = 0.0
 981           electrolyte(je,jsolid,ibin)  = 0.0
 982           electrolyte(je,jliquid,ibin) = 0.0
 983           activity(je,ibin)            = 0.0
 984           gam(je,ibin)                 = 0.0
 985        enddo
 986
 987        gam_ratio(ibin)   = 0.0
 988
 989        do iv = 1, ngas_ioa
 990           flux_s(iv,ibin)   = 0.0
 991           flux_l(iv,ibin)   = 0.0
 992           kg(iv,ibin)       = 0.0
 993           phi_volatile_s(iv,ibin) = 0.0
 994           phi_volatile_l(iv,ibin) = 0.0
 995           volatile_s(iv,ibin) = 0.0
 996        enddo
 997
 998
 999        jaerosolstate(ibin) = -1     ! initialize to default value
1000        jphase(ibin) = 0
1001
1002        do jc = 1, ncation
1003           mc(jc,ibin) = 0.0
1004        enddo
1005
1006        do ja = 1, nanion
1007           ma(ja,ibin) = 0.0
1008        enddo
1009
1010     enddo   ! ibin
1011
1012
1013
1014     return
1015   end subroutine initialize_mosaic_variables
1016
1017
1018
1019   subroutine overall_massbal_in( aer, gas, gas_netprod_otrproc, dtchem,            & !intent-ins
1020        total_species, tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in,    & !intent-outs
1021        tot_ca_in )
1022
1023
1024     use module_data_mosaic_aero, only: ngas_aerchtot, ngas_volatile,               &!
1025          naer, nbin_a_max, jtotal, nbin_a,                                         &!Input
1026          ih2so4_g,ihno3_g,ihcl_g,inh3_g,iso4_a,ino3_a,icl_a,inh4_a,ina_a,ica_a      !TBD
1027
1028     implicit none
1029
1030     !Subroutine Arguments
1031     real(r8), intent(in), dimension(ngas_aerchtot) :: gas, gas_netprod_otrproc
1032     real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
1033     real(r8), intent(in) :: dtchem
1034
1035     real(r8), intent(out) :: tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in
1036     real(r8), intent(out), dimension(ngas_volatile) ::total_species
1037
1038
1039     !local Variables
1040     integer ibin
1041
1042     tot_so4_in = gas(ih2so4_g)
1043     tot_no3_in = gas(ihno3_g)
1044     tot_cl_in  = gas(ihcl_g)
1045     tot_nh4_in = gas(inh3_g)
1046     tot_na_in  = 0.0
1047     tot_ca_in  = 0.0
1048
1049     tot_so4_in = gas(ih2so4_g) + max( gas_netprod_otrproc(ih2so4_g)*dtchem, 0.0_r8 )
1050
1051
1052     do ibin = 1, nbin_a
1053        tot_so4_in = tot_so4_in + aer(iso4_a,jtotal,ibin)
1054        tot_no3_in = tot_no3_in + aer(ino3_a,jtotal,ibin)
1055        tot_cl_in  = tot_cl_in  + aer(icl_a, jtotal,ibin)
1056        tot_nh4_in = tot_nh4_in + aer(inh4_a,jtotal,ibin)
1057        tot_na_in  = tot_na_in  + aer(ina_a,jtotal,ibin)
1058        tot_ca_in  = tot_ca_in  + aer(ica_a,jtotal,ibin)
1059     enddo
1060
1061
1062     total_species(inh3_g) = tot_nh4_in
1063     total_species(ihno3_g)= tot_no3_in
1064     total_species(ihcl_g) = tot_cl_in
1065
1066
1067     return
1068   end subroutine overall_massbal_in
1069
1070
1071
1072   subroutine overall_massbal_out( iprint_input, mbin, isteps_ASTEM, aer, gas, &
1073     tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in )
1074     !      include 'v33com'
1075     !      include 'v33com3'
1076     !      include 'v33com9a'
1077     !      include 'v33com9b'
1078     use module_data_mosaic_aero, only: ngas_aerchtot,naer,nbin_a_max,jtotal,    &!Parameters
1079          mYES,mNO,                                                                 &!Parameters
1080          nbin_a,                                                                   &!Input
1081          ih2so4_g,ihno3_g,ihcl_g,inh3_g,iso4_a,ino3_a,icl_a,inh4_a,ina_a,ica_a      !TBD
1082
1083     implicit none
1084
1085     ! subr. agrument
1086
1087     integer, intent(in)    ::  mbin, isteps_ASTEM
1088     integer, intent(inout) ::  iprint_input
1089
1090     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1091     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1092     real(r8), intent(inout) :: tot_so4_in, tot_no3_in, tot_cl_in, tot_nh4_in, tot_na_in, tot_ca_in
1093
1094     ! local variables
1095     integer ibin
1096     real(r8) :: tot_so4_out, tot_no3_out, tot_cl_out, tot_nh4_out, tot_na_out, tot_ca_out
1097     real(r8) :: diff_so4, diff_no3, diff_cl, diff_nh4, diff_na, diff_ca
1098     real(r8) :: reldiff_so4, reldiff_no3, reldiff_cl, reldiff_nh4, reldiff_na, reldiff_ca
1099
1100
1101
1102     tot_so4_out = gas(ih2so4_g)
1103     tot_no3_out = gas(ihno3_g)
1104     tot_cl_out  = gas(ihcl_g)
1105     tot_nh4_out = gas(inh3_g)
1106     tot_na_out  = 0.0
1107     tot_ca_out  = 0.0
1108
1109     do ibin = 1, nbin_a
1110        tot_so4_out = tot_so4_out + aer(iso4_a,jtotal,ibin)
1111        tot_no3_out = tot_no3_out + aer(ino3_a,jtotal,ibin)
1112        tot_cl_out  = tot_cl_out  + aer(icl_a,jtotal,ibin)
1113        tot_nh4_out = tot_nh4_out + aer(inh4_a,jtotal,ibin)
1114        tot_na_out  = tot_na_out  + aer(ina_a,jtotal,ibin)
1115        tot_ca_out  = tot_ca_out  + aer(ica_a,jtotal,ibin)
1116     enddo
1117
1118     diff_so4 = tot_so4_out - tot_so4_in
1119     diff_no3 = tot_no3_out - tot_no3_in
1120     diff_cl  = tot_cl_out  - tot_cl_in
1121     diff_nh4 = tot_nh4_out - tot_nh4_in
1122     diff_na  = tot_na_out  - tot_na_in
1123     diff_ca  = tot_ca_out  - tot_ca_in
1124
1125
1126     reldiff_so4 = 0.0
1127     if(tot_so4_in .gt. 1.e-25 .or. tot_so4_out .gt. 1.e-25)then
1128        reldiff_so4 = diff_so4/max(tot_so4_in, tot_so4_out)
1129     endif
1130
1131     reldiff_no3 = 0.0
1132     if(tot_no3_in .gt. 1.e-25 .or. tot_no3_out .gt. 1.e-25)then
1133        reldiff_no3 = diff_no3/max(tot_no3_in, tot_no3_out)
1134     endif
1135
1136     reldiff_cl = 0.0
1137     if(tot_cl_in .gt. 1.e-25 .or. tot_cl_out .gt. 1.e-25)then
1138        reldiff_cl = diff_cl/max(tot_cl_in, tot_cl_out)
1139     endif
1140
1141     reldiff_nh4 = 0.0
1142     if(tot_nh4_in .gt. 1.e-25 .or. tot_nh4_out .gt. 1.e-25)then
1143        reldiff_nh4 = diff_nh4/max(tot_nh4_in, tot_nh4_out)
1144     endif
1145
1146     reldiff_na = 0.0
1147     if(tot_na_in .gt. 1.e-25 .or. tot_na_out .gt. 1.e-25)then
1148        reldiff_na = diff_na/max(tot_na_in, tot_na_out)
1149     endif
1150
1151     reldiff_ca = 0.0
1152     if(tot_ca_in .gt. 1.e-25 .or. tot_ca_out .gt. 1.e-25)then
1153        reldiff_ca = diff_ca/max(tot_ca_in, tot_ca_out)
1154     endif
1155
1156     if( abs(reldiff_so4) .gt. 1.e-4 .or.   &
1157          abs(reldiff_no3) .gt. 1.e-4 .or.   &
1158          abs(reldiff_cl)  .gt. 1.e-4 .or.   &
1159          abs(reldiff_nh4) .gt. 1.e-4 .or.   &
1160          abs(reldiff_na)  .gt. 1.e-4 .or.   &
1161          abs(reldiff_ca)  .gt. 1.e-4)then
1162
1163
1164        if(iprint_input .eq. mYES)then
1165           write(6,*)'*** mbin = ', mbin, '  isteps = ', isteps_ASTEM
1166           write(6,*)'reldiff_so4 = ', reldiff_so4
1167           write(6,*)'reldiff_no3 = ', reldiff_no3
1168           write(6,*)'reldiff_cl  = ', reldiff_cl
1169           write(6,*)'reldiff_nh4 = ', reldiff_nh4
1170           write(6,*)'reldiff_na  = ', reldiff_na
1171           write(6,*)'reldiff_ca  = ', reldiff_ca
1172           !          call print_input
1173           iprint_input = mNO
1174        endif
1175
1176        !         stop
1177
1178     endif
1179
1180
1181     return
1182   end subroutine overall_massbal_out
1183
1184
1185
1186   !***********************************************************************
1187   ! checks if aerosol mass is too low to be of any significance
1188   ! and determine jaerosolstate
1189   !
1190   ! author: Rahul A. Zaveri
1191   ! update: jan 2005
1192   !-----------------------------------------------------------------------
1193   subroutine check_aerosol_mass(ibin, jaerosolstate,jphase,aer,num_a, mass_dry_a )
1194
1195     use module_data_mosaic_aero, only: nbin_a_max,naer,jtotal,no_aerosol,       &!Parameters
1196          mass_cutoff, mw_aer_mac,                                               &!Parameters
1197          iso4_a,ino3_a,icl_a,imsa_a,ico3_a,ica_a,ina_a,inh4_a                       !TBD
1198
1199     implicit none
1200
1201     !Intent-ins
1202     integer, intent(in) :: ibin
1203     real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
1204
1205     !Intent-inouts
1206     integer, intent(inout), dimension(nbin_a_max) :: jaerosolstate,jphase
1207     real(r8), intent(inout), dimension(nbin_a_max) :: num_a, mass_dry_a
1208
1209     !Local variables
1210     integer iaer
1211     real(r8) :: drymass, aer_H
1212
1213     mass_dry_a(ibin) = 0.0
1214
1215     aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
1216          aer(ino3_a,jtotal,ibin) +   &
1217          aer(icl_a,jtotal,ibin)  +   &
1218          aer(imsa_a,jtotal,ibin) +   &
1219          2.*aer(ico3_a,jtotal,ibin))-   &
1220          (2.*aer(ica_a,jtotal,ibin)  +   &
1221          aer(ina_a,jtotal,ibin)  +   &
1222          aer(inh4_a,jtotal,ibin))
1223     aer_H = max(aer_H, 0.0d0)
1224
1225     do iaer = 1, naer
1226        mass_dry_a(ibin) = mass_dry_a(ibin) +   &
1227             aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)  ! ng/m^3(air)
1228     enddo
1229     mass_dry_a(ibin) = mass_dry_a(ibin) + aer_H
1230
1231     drymass = mass_dry_a(ibin)                      ! ng/m^3(air)
1232     mass_dry_a(ibin) = mass_dry_a(ibin)*1.e-15      ! g/cc(air)
1233
1234     if(drymass .lt. mass_cutoff)then                        ! bin mass is too small
1235        jaerosolstate(ibin) = no_aerosol
1236        jphase(ibin) = 0
1237        if(drymass .eq. 0.)num_a(ibin) = 0.0
1238     endif
1239
1240     return
1241   end subroutine check_aerosol_mass
1242
1243
1244
1245   !***********************************************************************
1246   ! checks and conforms number according to the mass and bin size range
1247   !
1248   ! author: Rahul A. Zaveri
1249   ! update: jan 2005
1250   !-----------------------------------------------------------------------
1251   subroutine conform_aerosol_number(ibin,jaerosolstate,aer,num_a,vol_dry_a, Dp_dry_a)
1252
1253     use module_data_mosaic_constants,  only: pi
1254     use module_data_mosaic_aero,  only: nbin_a_max,naer,mSECTIONAL,no_aerosol,  &!Parameters
1255          jtotal, mw_aer_mac,dens_aer_mac,                                       &!Parameters
1256          msize_framework,                                                       &!Input
1257          iso4_a,ino3_a,icl_a,imsa_a,ico3_a,ica_a,ina_a,inh4_a                    !TBD
1258
1259     use module_data_mosaic_asecthp, only:isize_of_ibin,itype_of_ibin,volumlo_sect,&!TBD
1260          volumhi_sect                                                               !TBD
1261
1262     implicit none
1263
1264     !Intent-ins
1265     integer, intent(in) :: ibin
1266     integer, intent(in), dimension(nbin_a_max) :: jaerosolstate
1267
1268     real(r8), intent(in), dimension(nbin_a_max) :: Dp_dry_a
1269     real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
1270
1271     !intent-inout
1272     real(r8), intent(inout), dimension(nbin_a_max) :: num_a,vol_dry_a
1273
1274
1275     !Local variables
1276     integer :: iaer, isize, itype
1277     real(r8) :: num_at_dlo, num_at_dhi, numold
1278     real(r8) :: aer_H
1279     logical, parameter :: nonsect_set_number_always = .false.
1280
1281     ! when msize_framework = munstructured or mmodal,
1282     !    calculate number from volume concentration and mean dry diameter
1283     !       only when num_a(ibin) <= 0.0
1284     !    this should only happen at the very start of the simulation
1285     ! when msize_framework = msectional,
1286     !    check that mean dry diameter falls within the section/bin limits,
1287     !    and adjust number is this is not true
1288     if (msize_framework /= msectional) then
1289        if (num_a(ibin) > 0.0) return
1290     end if
1291
1292     vol_dry_a(ibin)  = 0.0          ! initialize to 0.0
1293
1294     if(jaerosolstate(ibin) .eq. no_aerosol) return
1295
1296
1297     ! calculate dry volume concentration
1298     aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
1299                 aer(ino3_a,jtotal,ibin) +   &
1300                 aer(icl_a,jtotal,ibin)  +   &
1301                 aer(imsa_a,jtotal,ibin) +   &
1302              2.*aer(ico3_a,jtotal,ibin))-   &
1303             (2.*aer(ica_a,jtotal,ibin)  +   &
1304                 aer(ina_a,jtotal,ibin)  +   &
1305                 aer(inh4_a,jtotal,ibin))
1306     aer_H = max(aer_H, 0.0d0)
1307
1308     do iaer = 1, naer
1309        vol_dry_a(ibin) = vol_dry_a(ibin) +   &
1310             aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer)  ! ncc/m^3(air)
1311     enddo
1312     vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H
1313     vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15                           ! cc(aer)/cc(air)
1314
1315
1316     if (msize_framework /= msectional) then
1317        ! unstructured or modal - set (initialize) number only when incoming value is zero
1318        if (num_a(ibin) <= 0.0) then
1319           num_a(ibin) = vol_dry_a(ibin)/((pi/6.0_r8)*Dp_dry_a(ibin)**3)       ! #/cc(air)
1320        end if
1321     else
1322        ! sectional
1323        if (num_a(ibin) <= 0.0) then
1324           ! in this case, num_a has probably not yet been initialized, so do it
1325           num_a(ibin) = vol_dry_a(ibin)/((pi/6.0_r8)*Dp_dry_a(ibin)**3)       ! #/cc(air)
1326        else
1327           ! in this case, check that bin mean size is within bounds
1328           isize = isize_of_ibin( ibin )
1329           itype = itype_of_ibin( ibin )
1330           num_at_dlo = vol_dry_a(ibin)/volumlo_sect(isize,itype)
1331           num_at_dhi = vol_dry_a(ibin)/volumhi_sect(isize,itype)
1332           numold = num_a(ibin)
1333           num_a(ibin) = min( num_a(ibin), num_at_dlo )
1334           num_a(ibin) = max( num_a(ibin), num_at_dhi )
1335        end if
1336     end if
1337
1338
1339     return
1340   end subroutine conform_aerosol_number
1341
1342
1343
1344   !***********************************************************************
1345   ! calculates dry density
1346   !
1347   ! author: Rahul A. Zaveri
1348   ! update: apr 2010
1349   !-----------------------------------------------------------------------
1350   subroutine calc_aerosol_dry_density(ibin,aer,dens_dry_a)
1351     !      include 'v33com9a'
1352
1353     use module_data_mosaic_aero,  only: nbin_a_max,naer,jtotal,                 &!Parameters
1354          inh4_a,ina_a,ica_a,ico3_a,imsa_a,icl_a,ino3_a,iso4_a,                  & !TBD
1355          mw_aer_mac, dens_aer_mac
1356
1357     !use module_data_mosaic_asecthp, only:                                            !BSINGH - not needed
1358
1359     implicit none
1360
1361     !Intent-in
1362     integer, intent(in) :: ibin
1363     real(r8), intent(in), dimension(naer,3,nbin_a_max) :: aer
1364
1365     !Intent-inout
1366     real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a
1367
1368     ! local variables
1369     integer :: iaer
1370     real(r8) :: aer_H
1371     real(r8) :: tmpa, tmp_volu, tmp_mass
1372
1373
1374     ! calculate dry volume concentration
1375     aer_H = ( 2.*max( 0.0_r8, aer(iso4_a,jtotal,ibin) ) +   &
1376                  max( 0.0_r8, aer(ino3_a,jtotal,ibin) ) +   &
1377                  max( 0.0_r8, aer(icl_a,jtotal,ibin) )  +   &
1378                  max( 0.0_r8, aer(imsa_a,jtotal,ibin) ) +   &
1379               2.*max( 0.0_r8, aer(ico3_a,jtotal,ibin) ) )   &
1380           - ( 2.*max( 0.0_r8, aer(ica_a,jtotal,ibin) )  +   &
1381                  max( 0.0_r8, aer(ina_a,jtotal,ibin) )  +   &
1382                  max( 0.0_r8, aer(inh4_a,jtotal,ibin) ) )
1383     aer_H = max( aer_H, 0.0_r8 )
1384
1385     tmp_mass = aer_H
1386     tmp_volu = aer_H   ! assume density=1.0 for H+
1387
1388     do iaer = 1, naer
1389        tmpa = max( 0.0_r8, aer(iaer,jtotal,ibin) ) * mw_aer_mac(iaer)
1390        tmp_mass = tmp_mass + tmpa                     !  ng/m^3(air)
1391        tmp_volu = tmp_volu + tmpa/dens_aer_mac(iaer)  ! ncc/m^3(air)
1392     enddo
1393
1394     ! the 1.0e-20 ng/m3 cutoff here is equivalent to the
1395     !     1.0e-35 g/cm3 cutoff used in mosaic_dynamic_solver
1396     if (min(tmp_mass,tmp_volu) >= 1.0e-20) then
1397        dens_dry_a(ibin) = tmp_mass/tmp_volu   ! g/cc
1398     else
1399        dens_dry_a(ibin) = 1.0
1400     end if
1401
1402     return
1403   end subroutine calc_aerosol_dry_density
1404
1405
1406
1407   !***********************************************************************
1408   ! updates/conforms size (diameter) according to the mass and number
1409   !
1410   ! author: Rahul A. Zaveri
1411   ! update: oct 2005
1412   !-----------------------------------------------------------------------
1413   subroutine conform_aerosol_size( ibin, jaerosolstate, aer, num_a, Dp_dry_a,     &
1414        vol_dry_a, mw_aer_mac, dens_aer_mac, mosaic_vars_aa )   ! TOUCH
1415
1416     !      include 'v33com9a'
1417     use module_data_mosaic_constants,  only : piover6,  third
1418     use module_data_mosaic_aero,  only : nbin_a_max, naer, no_aerosol, jtotal,       &!Parameters
1419          inh4_a, ina_a, ica_a, ico3_a, imsa_a, icl_a, ino3_a, iso4_a,                &
1420          mosaic_vars_aa_type                                                     !TBD
1421
1422     implicit none
1423
1424     ! subr arguments
1425     integer, intent(in):: ibin
1426     integer, intent(in), dimension(nbin_a_max) :: jaerosolstate
1427
1428     real(r8), intent(in), dimension(nbin_a_max) :: num_a
1429     real(r8), intent(in), dimension(naer) :: mw_aer_mac,dens_aer_mac
1430     real(r8), intent(inout), dimension(nbin_a_max) ::        Dp_dry_a,vol_dry_a
1431     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1432
1433     type (mosaic_vars_aa_type), intent(inout) :: mosaic_vars_aa
1434
1435     ! local variables
1436     integer iaer
1437     real(r8) :: num_at_dlo, num_at_dhi
1438     real(r8) :: aer_H
1439
1440
1441     vol_dry_a(ibin)  = 0.0          ! initialize to 0.0
1442
1443     if(jaerosolstate(ibin) .eq. no_aerosol) return
1444
1445     aer_H = (2.*aer(iso4_a,jtotal,ibin) +   &
1446                 aer(ino3_a,jtotal,ibin) +   &
1447                 aer(icl_a,jtotal,ibin)  +   &
1448                 aer(imsa_a,jtotal,ibin) +   &
1449              2.*aer(ico3_a,jtotal,ibin))-   &
1450             (2.*aer(ica_a,jtotal,ibin)  +   &
1451                 aer(ina_a,jtotal,ibin)  +   &
1452                 aer(inh4_a,jtotal,ibin))
1453     aer_H = max(aer_H, 0.0d0)
1454     do iaer = 1, naer
1455        vol_dry_a(ibin) = vol_dry_a(ibin) +   &
1456             aer(iaer,jtotal,ibin)*mw_aer_mac(iaer)/dens_aer_mac(iaer)       ! ng/m^3(air)
1457     enddo
1458     vol_dry_a(ibin) = vol_dry_a(ibin) + aer_H
1459     vol_dry_a(ibin) = vol_dry_a(ibin)*1.e-15                                ! cc(aer)/cc(air)
1460
1461
1462     ! update size
1463     !
1464     ! Box-model only
1465
1466     mosaic_vars_aa%f_mos_fail = -1
1467     if(vol_dry_a(ibin)<0.0_r8) then
1468        write(202,*)'EXITING due to negative vol_dry_a(',ibin,')=', &
1469           vol_dry_a(ibin), mosaic_vars_aa%it_mosaic, mosaic_vars_aa%hostgridinfo(1:3)
1470        mosaic_vars_aa%f_mos_fail = 1
1471        return
1472     endif
1473     Dp_dry_a(ibin) = (vol_dry_a(ibin)/(piover6*num_a(ibin)))**third
1474
1475     return
1476   end subroutine conform_aerosol_size
1477
1478
1479
1480   !***********************************************************************
1481   ! computes MTEM ternary parameters only once per transport time-step
1482   ! for a given aH2O (= RH)
1483   !
1484   ! author: Rahul A. Zaveri
1485   ! update: jan 2005
1486   ! reference: Zaveri, R.A., R.C. Easter, and A.S. Wexler,
1487   ! A new method for multicomponent activity coefficients of electrolytes
1488   ! in aqueous atmospheric aerosols, J. Geophys. Res., 2005.
1489   !-----------------------------------------------------------------------
1490   subroutine MTEM_compute_log_gamZ(aH2O,log_gamZ,b_mtem,aw_min)
1491     use module_data_mosaic_aero, only: nelectrolyte,                            &!Parameters
1492          jhno3,jnh4so4,jnh4no3,jnh4cl,jna2so4,jnano3,jnacl,jcano3,jcacl2,jhcl,  &
1493          jh2so4,jnh4hso4,jlvcite,jnahso4,jna3hso4,jhhso4                            !TBD
1494
1495     implicit none
1496
1497     !Sub args
1498     real(r8), intent(in) :: aH2O
1499     real(r8), intent(in), dimension(nelectrolyte) :: aw_min
1500     real(r8), intent(inout), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1501     real(r8), intent(in), dimension(6,nelectrolyte,nelectrolyte) :: b_mtem
1502     ! local variables
1503     integer jA
1504     ! functions
1505     !real(r8) :: fnlog_gamZ, bin_molality
1506
1507
1508     ! sulfate-poor species
1509     jA = jhno3
1510     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1511     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1512     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1513     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1514     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1515     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1516     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1517     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1518     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1519     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1520     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1521     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1522     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1523     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1524     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1525
1526
1527     jA = jhcl
1528     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1529     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1530     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1531     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1532     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1533     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1534     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1535     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1536     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1537     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1538     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1539     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1540     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1541     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1542     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1543
1544
1545     jA = jnh4so4
1546     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1547     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1548     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1549     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1550     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1551     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1552     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1553     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1554     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1555     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1556     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1557     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1558     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1559     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1560     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1561
1562
1563     jA = jnh4no3
1564     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1565     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1566     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1567     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1568     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1569     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1570     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1571     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1572     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1573     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1574
1575
1576     jA = jnh4cl
1577     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1578     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1579     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1580     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1581     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1582     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1583     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1584     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1585     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1586     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1587
1588
1589     jA = jna2so4
1590     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1591     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1592     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1593     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1594     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1595     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1596     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1597     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1598     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1599     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1600     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1601     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1602     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1603     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1604     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1605
1606
1607     jA = jnano3
1608     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1609     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1610     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1611     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1612     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1613     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1614     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1615     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1616     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1617     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1618
1619
1620     jA = jnacl
1621     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1622     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1623     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1624     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1625     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1626     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1627     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1628     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1629     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1630     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1631
1632
1633     jA = jcano3
1634     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1635     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1636     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1637     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1638     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1639     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1640     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1641     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1642     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1643     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1644
1645
1646     jA = jcacl2
1647     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1648     log_gamZ(jA,jnh4no3) = fnlog_gamZ(jA,jnh4no3,aH2O,b_mtem,aw_min)
1649     log_gamZ(jA,jnh4cl)  = fnlog_gamZ(jA,jnh4cl,aH2O,b_mtem,aw_min)
1650     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1651     log_gamZ(jA,jnano3)  = fnlog_gamZ(jA,jnano3,aH2O,b_mtem,aw_min)
1652     log_gamZ(jA,jnacl)   = fnlog_gamZ(jA,jnacl,aH2O,b_mtem,aw_min)
1653     log_gamZ(jA,jcano3)  = fnlog_gamZ(jA,jcano3,aH2O,b_mtem,aw_min)
1654     log_gamZ(jA,jcacl2)  = fnlog_gamZ(jA,jcacl2,aH2O,b_mtem,aw_min)
1655     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1656     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1657
1658
1659     ! sulfate-rich species
1660     jA = jh2so4
1661     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1662     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1663     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1664     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1665     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1666     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1667     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1668     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1669     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1670
1671
1672     jA = jhhso4
1673     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1674     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1675     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1676     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1677     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1678     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1679     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1680     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1681     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1682
1683
1684     jA = jnh4hso4
1685     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1686     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1687     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1688     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1689     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1690     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1691     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1692     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1693     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1694
1695
1696     jA = jlvcite
1697     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1698     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1699     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1700     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1701     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1702     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1703     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1704     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1705     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1706
1707
1708     jA = jnahso4
1709     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1710     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1711     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1712     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1713     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1714     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1715     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1716     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1717     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1718
1719
1720     jA = jna3hso4
1721     log_gamZ(jA,jh2so4)  = fnlog_gamZ(jA,jh2so4,aH2O,b_mtem,aw_min)
1722     log_gamZ(jA,jnh4hso4)= fnlog_gamZ(jA,jnh4hso4,aH2O,b_mtem,aw_min)
1723     log_gamZ(jA,jlvcite) = fnlog_gamZ(jA,jlvcite,aH2O,b_mtem,aw_min)
1724     log_gamZ(jA,jnh4so4) = fnlog_gamZ(jA,jnh4so4,aH2O,b_mtem,aw_min)
1725     log_gamZ(jA,jnahso4) = fnlog_gamZ(jA,jnahso4,aH2O,b_mtem,aw_min)
1726     log_gamZ(jA,jna3hso4)= fnlog_gamZ(jA,jna3hso4,aH2O,b_mtem,aw_min)
1727     log_gamZ(jA,jna2so4) = fnlog_gamZ(jA,jna2so4,aH2O,b_mtem,aw_min)
1728     log_gamZ(jA,jhno3)   = fnlog_gamZ(jA,jhno3,aH2O,b_mtem,aw_min)
1729     log_gamZ(jA,jhcl)    = fnlog_gamZ(jA,jhcl,aH2O,b_mtem,aw_min)
1730
1731     return
1732   end subroutine MTEM_compute_log_gamZ
1733
1734
1735
1736   subroutine degas_acids(jp,ibin,XT,aer,gas,electrolyte)
1737     use module_data_mosaic_aero, only: naer,nelectrolyte,nbin_a_max,            &
1738          ngas_aerchtot,jliquid,jsolid,jtotal,                                      &
1739          jhno3,jhcl,ihno3_g,ihcl_g,ino3_a,icl_a
1740
1741     implicit none
1742
1743     ! subr arguments
1744     integer, intent(in) :: jp, ibin
1745     real(r8), intent(in) :: XT
1746     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
1747     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
1748     real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
1749     ! local variables
1750     real(r8) :: ehno3, ehcl
1751
1752
1753
1754     if(jp .ne. jliquid)then
1755        write(6,*)'Error in degas_acids'
1756        write(6,*)'wrong jp'
1757     endif
1758
1759     ehno3 = electrolyte(jhno3,jp,ibin)
1760     ehcl  = electrolyte(jhcl,jp,ibin)
1761
1762     ! add to gas
1763     gas(ihno3_g) = gas(ihno3_g) + ehno3
1764     gas(ihcl_g)  = gas(ihcl_g)  + ehcl
1765
1766     ! remove from aer
1767     aer(ino3_a,jp,ibin) = aer(ino3_a,jp,ibin) - ehno3
1768     aer(icl_a, jp,ibin) = aer(icl_a, jp,ibin) - ehcl
1769
1770     ! update jtotal
1771     aer(ino3_a,jtotal,ibin) = aer(ino3_a,jliquid,ibin) +   &
1772          aer(ino3_a,jsolid, ibin)
1773
1774     aer(icl_a,jtotal,ibin)  = aer(icl_a,jliquid,ibin) +   &
1775          aer(icl_a,jsolid, ibin)
1776
1777     electrolyte(jhno3,jp,ibin) = 0.0
1778     electrolyte(jhcl,jp,ibin)  = 0.0
1779
1780     return
1781   end subroutine degas_acids
1782
1783
1784
1785
1786
1787   !***********************************************************************
1788   ! updates all temperature dependent thermodynamic parameters
1789   !
1790   ! author: Rahul A. Zaveri
1791   ! update: jan 2005
1792   !-----------------------------------------------------------------------
1793   subroutine update_thermodynamic_constants(     aH2O,         T_K_in,               & !intent-ins
1794        sat_soa,    aH2O_a,   log_gamZ,  Keq_sl,  sigma_water,  Kp_nh4cl,             & !intent-outs
1795        Kp_nh4no3,  Kp_nh3,   Keq_ll,    Keq_gl,  Keq_sg,       MDRH_T,               &
1796        molality0                                                                     )
1797
1798     use module_data_mosaic_aero, only: r8, nbin_a_max, ngas_volatile,  nelectrolyte,         &
1799          nrxn_aer_sg, nrxn_aer_gl, nrxn_aer_sl, nrxn_aer_ll, MDRH_T_NUM, d_mdrh_DIM2,        &
1800          nbin_a, b_mtem, b_zsr, a_zsr, aw_min, d_mdrh,                                       &
1801          jnh4so4, jlvcite, jnh4hso4, jnh4msa, jnh4no3, jnh4cl, jna2so4, jnahso4, jna3hso4,   &
1802          jnamsa, jnano3, jnacl, jcacl2, jcano3, jcamsa2, iaro1_g, iaro2_g, ialk1_g, iole1_g, &
1803          iapi1_g, iapi2_g, ilim1_g, ilim2_g, isoa_first, msoa_flag1,                         &
1804          use_cam5mam_soa_params
1805
1806     use module_mosaic_ext, only: bin_molality
1807
1808     use module_mosaic_soa_vbs, only: soa_vbs_update_thermcons
1809
1810     implicit none
1811
1812     !Subroutine Arguments
1813     real(r8), intent(in) :: aH2O, T_K_in
1814
1815     real(r8), intent(out) :: sigma_water,Kp_nh4cl,Kp_nh4no3,Kp_nh3
1816     real(r8), intent(out), dimension(nbin_a_max)    :: aH2O_a
1817     real(r8), intent(out), dimension(ngas_volatile) :: sat_soa
1818     real(r8), intent(out), dimension(nrxn_aer_ll)   :: Keq_ll
1819     real(r8), intent(out), dimension(nrxn_aer_sl)   :: Keq_sl
1820     real(r8), intent(out), dimension(nrxn_aer_gl)   :: Keq_gl
1821     real(r8), intent(out), dimension(nrxn_aer_sg)   :: Keq_sg
1822     real(r8), intent(out), dimension(MDRH_T_NUM)    :: MDRH_T
1823     real(r8), intent(out), dimension(nelectrolyte,nbin_a_max)  :: molality0 !BSINGH(05/23/2014) - Added dimension nbin_a_max
1824     real(r8), intent(out), dimension(nelectrolyte,nelectrolyte) :: log_gamZ
1825
1826     ! local variables
1827     integer :: ibin, iv, j_index, je
1828     logical :: use_sorgam_soa_species
1829
1830 !   real(r8) :: MWsoa
1831     real(r8) :: Po_soa(ngas_volatile)
1832     real(r8) :: rt
1833     real(r8) :: sat_factor
1834     real(r8) :: T_K             ! 6/3/2015 RAZ : bound temperature between 220 and 330 K
1835     real(r8) :: tr, term
1836     ! function
1837     !real(r8) :: fn_Keq, fn_Po, drh_mutual, bin_molality
1838
1839
1840 ! bound temperature between 220 and 330 K  ! 6/3/2015 RAZ
1841     T_K = max( 220.0_r8, T_K_in )
1842     T_K = min( 330.0_r8, T_K )
1843
1844     tr = 298.15                   ! reference temperature
1845     rt = 82.056*T_K/(1.e9*1.e6)   ! [m^3 atm/nmol]
1846
1847     ! gas-liquid
1848     Keq_gl(1)= 1.0                                        ! Kelvin Effect (default)
1849     Keq_gl(2)= fn_Keq(57.64d0, 13.79d0, -5.39d0,T_K)*rt     ! NH3(g)  <=> NH3(l)
1850     Keq_gl(3)= fn_Keq(2.63d6,  29.17d0, 16.83d0,T_K)*rt     ! HNO3(g) <=> NO3- + H+
1851     Keq_gl(4)= fn_Keq(2.00d6,  30.20d0, 19.91d0,T_K)*rt     ! HCl(g)  <=> Cl- + H+
1852
1853     ! liquid-liquid
1854     Keq_ll(1)= fn_Keq(1.0502d-2, 8.85d0, 25.14d0,T_K)      ! HSO4- <=> SO4= + H+
1855     Keq_ll(2)= fn_Keq(1.805d-5, -1.50d0, 26.92d0,T_K)      ! NH3(l) + H2O = NH4+ + OH-
1856     Keq_ll(3)= fn_Keq(1.01d-14,-22.52d0, 26.92d0,T_K)      ! H2O(l) <=> H+ + OH-
1857
1858
1859     Kp_nh3   = Keq_ll(3)/(Keq_ll(2)*Keq_gl(2))
1860     Kp_nh4no3= Kp_nh3/Keq_gl(3)
1861     Kp_nh4cl = Kp_nh3/Keq_gl(4)
1862
1863
1864     ! solid-gas
1865     Keq_sg(1)= fn_Keq(4.72d-17,-74.38d0,6.12d0,T_K)/rt**2  ! NH4NO3<=>NH3(g)+HNO3(g)
1866     Keq_sg(2)= fn_Keq(8.43d-17,-71.00d0,2.40d0,T_K)/rt**2  ! NH4Cl <=>NH3(g)+HCl(g)
1867
1868
1869     ! solid-liquid
1870     Keq_sl(jnh4so4) = fn_Keq(1.040d0,-2.65d0, 38.57d0, T_K)  ! amSO4(s) = 2NH4+ + SO4=
1871     Keq_sl(jlvcite) = fn_Keq(11.8d0, -5.19d0, 54.40d0, T_K)  ! lvcite(s)= 3NH4+ + HSO4- + SO4=
1872     Keq_sl(jnh4hso4)= fn_Keq(117.0d0,-2.87d0, 15.83d0, T_K)  ! amHSO4(s)= NH4+ + HSO4-
1873     Keq_sl(jnh4msa) = 1.e15                                      ! NH4MSA(s)= NH4+ + MSA-
1874     Keq_sl(jnh4no3) = fn_Keq(12.21d0,-10.4d0, 17.56d0, T_K)  ! NH4NO3(s)= NH4+ + NO3-
1875     Keq_sl(jnh4cl)  = fn_Keq(17.37d0,-6.03d0, 16.92d0, T_K)  ! NH4Cl(s) = NH4+ + Cl-
1876     Keq_sl(jna2so4) = fn_Keq(0.491d0, 0.98d0, 39.75d0, T_K)  ! Na2SO4(s)= 2Na+ + SO4=
1877     Keq_sl(jnahso4) = fn_Keq(313.0d0, 0.8d0,  14.79d0, T_K)  ! NaHSO4(s)= Na+ + HSO4-
1878     Keq_sl(jna3hso4)= 1.e15                                      ! Na3H(SO4)2(s) = 2Na+ + HSO4- + SO4=
1879     Keq_sl(jnamsa)  = 1.e15                                      ! NaMSA(s) = Na+ + MSA-
1880     Keq_sl(jnano3)  = fn_Keq(11.95d0,-8.22d0, 16.01d0, T_K)  ! NaNO3(s) = Na+ + NO3-
1881     Keq_sl(jnacl)   = fn_Keq(38.28d0,-1.52d0, 16.89d0, T_K)  ! NaCl(s)  = Na+ + Cl-
1882     Keq_sl(jcacl2)  = fn_Keq(8.0d11,  32.84d0,44.79d0, T_K)  ! CaCl2(s) = Ca++ + 2Cl-
1883     Keq_sl(jcano3)  = fn_Keq(4.31d5,   7.83d0,42.01d0, T_K)  ! Ca(NO3)2(s) = Ca++ + 2NO3-
1884     Keq_sl(jcamsa2) = 1.e15                                ! CaMSA2(s)= Ca+ + 2MSA-
1885
1886
1887
1888     ! soa parameters
1889     if (msoa_flag1 == 1) then
1890        use_sorgam_soa_species = .true.
1891     else
1892        use_sorgam_soa_species = .false.
1893     end if
1894
1895     ! vapor pressures of soa species
1896     Po_soa(:) = 1.0e5_r8  ! this default value gives csat=6.1e9 ug/m^3 (~1e9 ppbv) for 298K and MW=150
1897     if ( use_sorgam_soa_species ) then
1898     Po_soa(iaro1_g) = fn_Po(5.7d-5, 156.0d0, T_K) ! [Pascal]
1899     Po_soa(iaro2_g) = fn_Po(1.6d-3, 156.0d0, T_K) ! [Pascal]
1900     Po_soa(ialk1_g) = fn_Po(5.0d-6, 156.0d0, T_K) ! [Pascal]
1901     Po_soa(iole1_g) = fn_Po(5.0d-6, 156.0d0, T_K) ! [Pascal]
1902     Po_soa(iapi1_g) = fn_Po(4.0d-6, 156.0d0, T_K) ! [Pascal]
1903     Po_soa(iapi2_g) = fn_Po(1.7d-4, 156.0d0, T_K) ! [Pascal]
1904     Po_soa(ilim1_g) = fn_Po(2.5d-5, 156.0d0, T_K) ! [Pascal]
1905     Po_soa(ilim2_g) = fn_Po(1.2d-4, 156.0d0, T_K) ! [Pascal]
1906     end if
1907
1908     sat_soa(:) = 0.0_r8
1909     sat_factor = 0.5  ! = 1.0 for original SORGAM parameters
1910     do iv = isoa_first, ngas_volatile
1911        !       sat_soa(iv) = 1.e9*Po_soa(iv)/(8.314*T_K)  ! [nmol/m^3(air)]
1912        sat_soa(iv) = sat_factor * 1.e9*Po_soa(iv)/(8.314*T_K)  ! [nmol/m^3(air)]
1913     enddo
1914
1915     if ( ( use_cam5mam_soa_params > 0 ) .and. &
1916          ( 1 <= ilim2_g .and. ilim2_g <= ngas_volatile ) ) then
1917        Po_soa(ilim2_g) = fn_Po(1.0d-10, 156.0d0, T_K) ! [Pascal]
1918        sat_soa(ilim2_g) = 1.e9*Po_soa(ilim2_g)/(8.314*T_K)  ! [nmol/m^3(air)]
1919     end if
1920
1921     ! MWsoa = 120.0
1922     ! sat_soa(iapi1_g) = 1000.*2564.1/MWsoa ! [nmol/m^3(air)]
1923     ! sat_soa(iapi2_g) = 1000.*11.803/MWsoa ! [nmol/m^3(air)]
1924
1925     if (msoa_flag1 >= 1000) call soa_vbs_update_thermcons( t_k, po_soa, sat_soa )
1926
1927
1928
1929     ! water surface tension
1930     term = (647.15 - T_K)/647.15
1931     sigma_water = 0.2358*term**1.256 * (1. - 0.625*term) ! surface tension of pure water in N/m
1932
1933     ! MDRH(T)
1934     do j_index = 1, 63
1935        MDRH_T(j_index) = drh_mutual(j_index,T_K)
1936     enddo
1937
1938
1939
1940     ! RH dependent parameters
1941     do ibin = 1, nbin_a
1942        aH2O_a(ibin) = aH2O                        ! initialize
1943
1944       do je = 1, nelectrolyte
1945         molality0(je,ibin) = bin_molality(je,ibin,aH2O_a,b_zsr,a_zsr,aw_min)  ! compute aH2O dependent binary molalities. RAZ 5/20/2014
1946       enddo
1947
1948     enddo
1949
1950     call MTEM_compute_log_gamZ(aH2O,log_gamZ,b_mtem,aw_min)              ! function of aH2O and T
1951
1952
1953     return
1954   end subroutine update_thermodynamic_constants
1955
1956
1957
1958
1959
1960   !***********************************************************************
1961   ! functions used in MOSAIC
1962   !
1963   ! author: Rahul A. Zaveri
1964   ! update: jan 2005
1965   !-----------------------------------------------------------------------
1966
1967
1968
1969   !----------------------------------------------------------
1970   function fn_Keq(Keq_298, a, b, T)
1971     implicit none
1972     real(r8) :: fn_Keq
1973     ! subr. arguments
1974     real(r8) :: Keq_298, a, b, T
1975     ! local variables
1976     real(r8) :: tt
1977
1978
1979     tt = 298.15/T
1980     fn_Keq = Keq_298*exp(a*(tt-1.)+b*(1.+log(tt)-tt))
1981
1982     return
1983   end function fn_Keq
1984   !----------------------------------------------------------
1985
1986
1987
1988   !----------------------------------------------------------
1989   function fn_Po(Po_298, DH, T)   ! TOUCH
1990     implicit none
1991     real(r8) :: fn_Po
1992     ! subr. arguments
1993     real(r8) :: Po_298, DH, T
1994     ! local variables
1995
1996     fn_Po = Po_298*exp(-(DH/8.314e-3)*(1./T - 3.354016435e-3))
1997
1998     return
1999   end function fn_Po
2000   !----------------------------------------------------------
2001
2002
2003
2004   !----------------------------------------------------------
2005   function drh_mutual(j_index,T_K)            ! TOUCH
2006     use module_data_mosaic_aero, only: d_mdrh
2007
2008     implicit none
2009
2010
2011     !Subr. arguments
2012     integer,  intent(in) ::  j_index
2013     real(r8), intent(in) :: T_K
2014
2015     !Local variables
2016     integer j
2017     real(r8) :: drh_mutual
2018
2019     j = j_index
2020
2021     if(j_index .eq. 7 .or. j_index .eq. 8 .or.   &
2022          (j_index.ge. 34 .and. j_index .le. 51))then
2023
2024        drh_mutual = 10.0  ! cano3 or cacl2 containing mixtures
2025
2026     else
2027
2028        drh_mutual =  d_mdrh(j,1) + T_K*   &
2029             (d_mdrh(j,2) + T_K*   &
2030             (d_mdrh(j,3) + T_K*   &
2031             d_mdrh(j,4) )) + 1.0
2032
2033 ! bound drh_mutual between 0% and 100%  ! RAZ 6/3/2015
2034        drh_mutual = max(   0.0_r8, drh_mutual )
2035        drh_mutual = min( 100.0_r8, drh_mutual )
2036
2037     endif
2038
2039
2040     return
2041   end function drh_mutual
2042   !----------------------------------------------------------
2043
2044
2045 ! RAZ
2046 ! Moved the following code to module_mosaic_ext.f90
2047 ! function bin_molality
2048
2049
2050   !----------------------------------------------------------
2051   function fnlog_gamZ(jA,jE,aH2O,b_mtem,aw_min) ! jA in jE
2052     use module_data_mosaic_aero, only: nelectrolyte
2053
2054     implicit none
2055
2056     real(r8) :: fnlog_gamZ
2057     ! subr. arguments
2058     integer, intent(in) :: jA, jE
2059     real(r8), intent(in) :: aH2O
2060     real(r8), intent(in),dimension(nelectrolyte) :: aw_min
2061     real(r8), intent(in), dimension(6,nelectrolyte,nelectrolyte) :: b_mtem
2062     ! local variables
2063     real(r8) :: aw
2064
2065
2066     aw = max(aH2O, aw_min(jE))
2067
2068     fnlog_gamZ = b_mtem(1,jA,jE) + aw*   &
2069          (b_mtem(2,jA,jE) + aw*   &
2070          (b_mtem(3,jA,jE) + aw*   &
2071          (b_mtem(4,jA,jE) + aw*   &
2072          (b_mtem(5,jA,jE) + aw*   &
2073          b_mtem(6,jA,jE) ))))
2074
2075     return
2076   end function fnlog_gamZ
2077   !----------------------------------------------------------
2078
2079
2080
2081   !----------------------------------------------------------
2082   ! currently not used
2083   !
2084   ! two roots of a quadratic equation
2085   !
2086   subroutine quadratix(a,b,c, qx1,qx2)
2087     implicit none
2088     ! subr. arguments
2089     real(r8) :: a, b, c, qx1, qx2
2090     ! local variables
2091     real(r8) :: x, dum
2092
2093
2094     if(b .ne. 0.0)then
2095        x = 4.*(a/b)*(c/b)
2096     else
2097        x = 1.e+6
2098     endif
2099
2100     if(abs(x) .lt. 1.e-6)then
2101        dum = ( (0.5*x) +   &
2102             (0.125*x**2) +   &
2103             (0.0625*x**3) )
2104
2105        qx1 = (-0.5*b/a)*dum
2106        qx2 = -b/a - qx1
2107
2108     else
2109
2110        qx1 = ((-b)+sqrt((b*b)-(4.*a*c)))/   &
2111             (2.*a)
2112        qx2 = ((-b)-sqrt((b*b)-(4.*a*c)))/   &
2113             (2.*a)
2114
2115     endif
2116
2117     return
2118   end subroutine quadratix
2119
2120
2121
2122   !***********************************************************************
2123   ! computes aerosol optical properties
2124   !
2125   ! author: Rahul A. Zaveri
2126   ! update: jan 2005
2127   !-----------------------------------------------------------------------
2128   subroutine aerosol_optical_properties(                                 &
2129           gas, aer, num_a, water_a,                                      &
2130           dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &
2131           Dp_dry_a, Dp_wet_a, dp_core_a,                                 &
2132           ri_shell_a, ri_core_a, ri_avg_a, jaerosolstate, jphase,        &
2133           tot_cl_in, tot_nh4_in, tot_no3_in, XT, area_dry_a, area_wet_a, &
2134           dens_dry_a, dens_wet_a, mass_dry_a, mass_wet_a, vol_dry_a,     &
2135           vol_wet_a, total_species, electrolyte     )
2136
2137     use module_data_mosaic_aero, only: &
2138        icl_a, inh4_a, ino3_a, ihcl_g, inh3_g, ihno3_g, jtotal, &
2139        naer, naercomp, nbin_a, nbin_a_max, nelectrolyte,       &
2140        ngas_aerchtot, ngas_volatile, no_aerosol
2141     use module_mosaic_ext, only: calc_dry_n_wet_aerosol_props, &
2142          conform_electrolytes
2143
2144     implicit none
2145
2146     ! subr arguments
2147     integer,  intent(inout), dimension(nbin_a_max) :: jaerosolstate, jphase
2148
2149     real(r8), intent(in), dimension(naer)       :: dens_aer_mac, mw_aer_mac
2150     real(r8), intent(in), dimension(naercomp)   :: dens_comp_a,mw_comp_a
2151
2152     real(r8), intent(inout) :: tot_cl_in, tot_nh4_in, tot_no3_in, XT
2153     real(r8), intent(inout), dimension(ngas_aerchtot) :: gas
2154     real(r8), intent(inout), dimension(naer,3,nbin_a_max) :: aer
2155     real(r8), intent(inout), dimension(nbin_a_max) :: water_a
2156     real(r8), intent(inout), dimension(nbin_a_max) :: num_a
2157     real(r8), intent(inout), dimension(nbin_a_max) :: Dp_dry_a
2158     real(r8), intent(inout), dimension(nbin_a_max) :: Dp_wet_a, dp_core_a
2159     real(r8), intent(inout), dimension(nbin_a_max) :: area_dry_a, area_wet_a
2160     real(r8), intent(inout), dimension(nbin_a_max) :: dens_dry_a, dens_wet_a
2161     real(r8), intent(inout), dimension(nbin_a_max) :: mass_dry_a, mass_wet_a
2162     real(r8), intent(inout), dimension(nbin_a_max) :: vol_dry_a, vol_wet_a
2163     real(r8), intent(inout), dimension(ngas_volatile) :: total_species
2164     real(r8), intent(inout), dimension(nelectrolyte,3,nbin_a_max) :: electrolyte
2165
2166
2167     complex,  intent(in), dimension(naercomp)   :: ref_index_a
2168     complex,  intent(inout), dimension(nbin_a_max) :: ri_shell_a, ri_avg_a, ri_core_a
2169
2170     ! local variables
2171     integer iaer, ibin, je, k
2172
2173     ! initialize to zero
2174     do ibin = 1, nbin_a
2175        do je = 1, nelectrolyte
2176           electrolyte(je,jtotal,ibin)  = 0.0
2177        enddo
2178        jaerosolstate(ibin) = -1   ! initialize to default value
2179     enddo
2180
2181     ! calc total_species for conform_electrolytes call
2182     total_species(:) = 0.0_r8
2183     tot_no3_in = gas(ihno3_g)
2184     tot_cl_in  = gas(ihcl_g)
2185     tot_nh4_in = gas(inh3_g)
2186     do ibin = 1, nbin_a
2187        tot_no3_in = tot_no3_in + aer(ino3_a,jtotal,ibin)
2188        tot_cl_in  = tot_cl_in  + aer(icl_a, jtotal,ibin)
2189        tot_nh4_in = tot_nh4_in + aer(inh4_a,jtotal,ibin)
2190     enddo
2191     total_species(inh3_g) = tot_nh4_in
2192     total_species(ihno3_g)= tot_no3_in
2193     total_species(ihcl_g) = tot_cl_in
2194
2195
2196     ! calc properties for each bin
2197     do  ibin = 1, nbin_a
2198
2199        call check_aerosol_mass( ibin, jaerosolstate, jphase, aer, num_a, mass_dry_a )
2200
2201        if(jaerosolstate(ibin) .ne. no_aerosol) then
2202
2203           ! conforms aer(jtotal) to a valid aerosol
2204           call conform_electrolytes( jtotal, ibin, XT, aer, gas, electrolyte, total_species, tot_cl_in )
2205
2206           ! check mass again after conform_electrolytes
2207           call check_aerosol_mass( ibin, jaerosolstate, jphase, aer, num_a, mass_dry_a )
2208
2209           if(jaerosolstate(ibin) .ne. no_aerosol) then
2210              ! adjusts number conc so that it conforms with bin mass and diameter
2211              call conform_aerosol_number( ibin, jaerosolstate, aer, num_a, vol_dry_a, Dp_dry_a)
2212
2213              ! calc Dp_wet, ref index
2214              call calc_dry_n_wet_aerosol_props(                                &
2215                   ibin, jaerosolstate, aer, electrolyte, water_a, num_a,         &  ! input
2216                   dens_comp_a, mw_comp_a, dens_aer_mac, mw_aer_mac, ref_index_a, &  ! input
2217                   Dp_dry_a, Dp_wet_a, dp_core_a,                                 &  ! output
2218                   area_dry_a, area_wet_a, mass_dry_a, mass_wet_a,                &  ! output
2219                   vol_dry_a, vol_wet_a, dens_dry_a, dens_wet_a,                  &  ! output
2220                   ri_shell_a, ri_core_a, ri_avg_a                                )  ! output
2221           endif
2222        endif
2223
2224     enddo
2225
2226     return
2227   end subroutine aerosol_optical_properties
2228
2229
2230 end module module_mosaic_box_aerchem
2231