!---------------------------------------------------------------- ! SOA formation from glyoxal with complex formulation including ! * reversible formation (Kampf et al., ES&T, submitted) ! * dark/ammonium-catalyzed formation (Noziere, J. Phys. Chem., 2009) ! * OH chemistry (Ervens and Volkamer, ACP, 2010) ! * surface uptake (Ervens and Volkamer, ACP, 2010) ! Christoph Knote, ACD, NCAR, 20130326 !---------------------------------------------------------------- MODULE module_mosaic_gly IMPLICIT NONE INTEGER, PARAMETER :: nspecs = 13, & igly_g = 1, & igly_r1 = 2, & igly_r2 = 3, & igly_nh4 = 4, & igly_sfc = 5, & igly_oh = 6, & ic_as = 7, & ic_an = 8, & ia_nh4 = 9, & ioh_g = 10, & iph = 11, & iwater = 12, & iarea = 13 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> LOGICAL, PARAMETER :: lfast_tau1 = .FALSE. ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CONTAINS ! adapted from Numerical Recipes, Second Edition, p. 706-707 ! changes: added "h" (timestep) argument to DERIVS calls SUBROUTINE rk4(y, dydx, n, x, h, yout, derivs) INTEGER n REAL(kind=8) :: h, x, dydx(n), y(n), yout(n) EXTERNAL derivs INTEGER i REAL(kind=8) :: h6, hh, xh, dym(nspecs), dyt(nspecs), yt(nspecs) hh=h*0.5 h6=h/6. xh=x+hh DO i=1, n yt(i) = y(i) + hh * dydx(i) ENDDO CALL derivs(xh, yt, h, dyt) DO i=1, n yt(i) = y(i) + hh * dyt(i) ENDDO CALL derivs(xh, yt, h, dym) DO i=1, n yt(i) = y(i) + h * dym(i) dym(i) = dyt(i) + dym(i) ENDDO CALL derivs(x+h, yt, h, dyt) DO i=1, n yout(i) = y(i) + h6 * ( dydx(i) + dyt(i) + 2. * dym(i)) ENDDO RETURN END SUBROUTINE rk4 ! Simple SOA formation from glyoxal as presented in ! Washenfelder et al, JGR, 2011 SUBROUTINE glysoa_simple(dtchem) USE module_data_mosaic_therm, ONLY: t_k, area_wet_a, gas, aer, & jtotal, igly, iglysoa_sfc_a, nbin_a IMPLICIT NONE REAL(kind=8), INTENT(IN) :: dtchem REAL(kind=8) :: omega, gamma_gly, A, delta_gly, frac_A INTEGER :: ibin ! mean molecular velocity of glyoxal (cm/s) omega = 1.455e4 * sqrt(t_k / 58.0_8) ! aerosol uptake coefficient for glyoxal (-) ! Washenfelder et al., 2011: ! 0 - 8 x 10^-4 ! 2 x 10^-4 (+/- 1 x 10^-4) ! Volkamer et al., 2007 ! 3.7 x 10^-3 ! B. Ervens, pers. comm., 2010 gamma_gly = 3.3E-3 ! get total aerosol surface area (cm^2 / cm^3) A = 0.0 DO ibin = 1, nbin_a A = A + area_wet_a(ibin) ENDDO ! no aerosol surface area - no uptake IF (A > 0.0) THEN ! first order uptake, Fuchs and Sutugin, 1971 ! dCg = 1/4 * gamma * A * |v_mol| * Cg * dt delta_gly = 0.25 * gamma_gly * A * omega * gas(igly) * dtchem ! avoid negative concentrations delta_gly = MIN(gas(igly), delta_gly) ! update partitioning gas(igly) = gas(igly) - delta_gly ! distribute onto bins according to fraction of surface area DO ibin = 1, nbin_a frac_A = area_wet_a(ibin) / A ! we take the "photochemical" glysoa aerosol as surrogate aer(iglysoa_sfc_a, jtotal, ibin) = aer(iglysoa_sfc_a, jtotal, ibin) & + frac_A * delta_gly ENDDO ENDIF END SUBROUTINE glysoa_simple SUBROUTINE glysoa_complex_derivs(x, y, dt, dydx) USE module_data_mosaic_therm, ONLY: conv1a, & ! converts q/mol(air) to nq/m^3 (q = mol or g) p_atm, & ! pressure (atm) t_k ! temperature (K) REAL(kind=8), INTENT(IN) :: x, y(nspecs), dt REAL(kind=8), INTENT(OUT) :: dydx(nspecs) REAL(kind=8), PARAMETER :: eps = 1.e-16 , & ! minimum allowed concentration in reservoirs Kh_water = 4.19e5 , & ! effective Henry's law constant of glyoxal in pure water (M atm-1) (Ip et al., GRL, 2011) Kh_oh = 25.0, & ! Henry's law constant of OH in pure water (M atm-1) (Klaening et al., 1985) k_oh = 1.1e9 ! OH reaction rate (mol L-1 s-1) (Ervens and Volkamer, ACP, 2010) REAL(kind=8) :: gly_g_atm, & ! gas-phase concentration in atm f_A1, & ! fraction of glyoxal in reservoir 1 tau1, & ! characteristical timescale reservoir 1 tau2, & ! characteristical timescale reservoir 2 oh_g_atm, & ! gas-phase OH concentration in atm oh_a, & ! liquid-phase OH concentration c_tot, & ! total concentration of dissolved salts (M) Kh_eq, & ! Henry's law constant at equilibrium (M atm-1) gly_ptot_eq,& ! total glyoxal concentration (reservoirs 1 and 2) at equilibrium gly_r1_eq, & ! glyoxal concentration (reservoir 1) at equilibrium gly_r2_eq, & ! glyoxal concentration (reservoir 2) at equilibrium anh4, & ! ammonium-ion activity (constrained) kII, kI, & ! second and first order dark rate constants omega ! mean molecular velocity of glyoxal (cm/s) ! tendencies REAL(kind=8) :: dg_r1, & ! from gas-phase to reservoir 1 dr1_r2, & ! reservoir 1 to reservoir 2 (or vice versa) dr1_nh4, & ! reservoir 1 to nh4 dr1_oh, & ! reservoir 1 to oh dg_sfc ! gas-phase to surface uptake REAL(kind=8) :: accloss, & ! acc. loss scaling ! tendency scaling dg_r1 = 0.0 dr1_r2 = 0.0 dr1_nh4 = 0.0 dr1_oh = 0.0 dg_sfc = 0.0 ! convert gas-phase glyoxal concentration ! ------------------------------------------------------------- gly_g_atm = y(igly_g) / conv1a ! mole / mole gly_g_atm = gly_g_atm * p_atm ! atm ! with A2/A1 = Kolig, and Kolig is 1 f_A1 = 0.5 tau1 = 2.5e2 ! s tau2 = 5.5e3 ! s IF ( y(ic_as) + y(ic_an) .GT. 12.0 ) THEN ! with A2/A1 = Kolig, and Kolig is 0.5 f_A1 = 0.6667 tau1 = 4.4e4 ! s tau2 = 4.7e4 ! s ENDIF ! Kampff et al., ES&T, 2013, submitted: kinetic limitation of ! salting-in for SO4 concentrations > 12 M. c_tot = MIN( 12.0, y(ic_as) + y(ic_an) ) ! effective Henry's law constant including salting-in effect (Kampff et al.) ! at equilibrium ! derived from eqn. 3 in Kampff et al., -0.24 is "salting-in" constant Kh_eq = Kh_water / 10**(-0.24D0 * c_tot) ! mol L-1 atm-1 == mol kg-1 atm-1 ! gly_g_atm in atm, Kh_eq in mol kg-1 atm-1, water in kg m-3 air ! total glyoxal concentration (reservoirs 1 and 2) at equilibrium gly_ptot_eq = gly_g_atm * Kh_eq * y(iwater) * 1e9 ! nmol / m3 air gly_r1_eq = gly_ptot_eq * f_A1 ! in reservoir 1 at equilibrium gly_r2_eq = gly_ptot_eq * (1.0 - f_A1) ! in reservoir 2 at equilibrium ! process tendencies in nmol m-3 s-1 ! from gas-phase to reservoir 1 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (.NOT. lfast_tau1) THEN ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< dg_r1 = (1.0/tau1) * (gly_r1_eq - y(igly_r1)) ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ENDIF ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! reservoir 1 to reservoir 2 dr1_r2 = (1.0/tau2) * (gly_r2_eq - y(igly_r2)) ! OH reaction (Ervens and Volkamer, Atm. Chem. Phys., 2010) ! ------------------------------------------------------------- oh_g_atm = y(ioh_g) / conv1a ! now its in mole / mole oh_g_atm = oh_g_atm * p_atm ! partial pressure (atm) oh_a = oh_g_atm * Kh_oh * y(iwater) * 1e9 ! nmol / m3 air dr1_oh = k_oh * y(igly_r1) * oh_a ! nmole m-3 s-1 ! dark pathway (Noziere et al., J. Phys. Chem., 2009) ! ------------------------------------------------------------- ! second-order rate constant (mol-1 kg s-1): anh4 = MAX( 0.0, MIN( 4.0, y(ia_nh4)) ) ! restrict to measured range kII = 2.e-10 * exp(1.5 * anh4) * exp(2.5 * y(iph)) ! dark process uses reservoir 1 kI = kII * y(igly_r1) / y(iwater) * 1e-9 ! s-1 dr1_nh4 = kI * y(igly_r1) ! nmol m-3 s-1 ! surface uptake (Ervens and Volkamer, ACP, 2010) ! ------------------------------------------------------------- ! mean molecular velocity of glyoxal (cm/s) omega = 1.455e4 * sqrt(t_k / 58.0D0) ! first order uptake, Fuchs and Sutugin, 1971 ! dCg = 1/4 * gamma * A * |v_mol| * Cg * dt, ! gamma downscaled to 1.0e-3 according to Waxman et al., 2013 dg_sfc = 0.25D0 * 1.e-3 * y(iarea) * omega * y(igly_g) ! Numerical integration ! ------------------------------------------------------------- ! check for undershoots, avoid negative concentrations ! while ensuring we don't loose mass IF ( y(igly_g) < eps ) THEN dg_r1 = 0.0 dg_sfc = 0.0 ELSE accloss = (dg_r1 + dg_sfc) * dt IF ( ( y(igly_g) - accloss ) < eps ) THEN scaling = y(igly_g) / (accloss + eps) dg_r1 = dg_r1 * scaling dg_sfc = dg_sfc * scaling ENDIF ENDIF IF ( y(igly_r1) < eps ) THEN dr1_r2 = 0.0 dr1_nh4 = 0.0 dr1_oh = 0.0 ELSE accloss = (-dg_r1 + dr1_r2 + dr1_nh4 + dr1_oh) * dt IF ( ( y(igly_r1) - accloss ) < eps) THEN scaling = y(igly_r1) / (accloss + eps) dr1_r2 = dr1_r2 * scaling dr1_nh4 = dr1_nh4 * scaling dr1_oh = dr1_oh * scaling ENDIF ENDIF IF ( y(igly_r2) < eps ) THEN dr1_r2 = MAX( 0.0, dr1_r2 ) ELSE accloss = -dr1_r2 * dt IF ( ( y(igly_r2) - accloss ) < eps ) THEN scaling = y(igly_r2) / (accloss + eps) dr1_r2 = dr1_r2 * scaling ENDIF ENDIF ! sum tendencies dydx(igly_g) = -dg_r1 -dg_sfc dydx(igly_r1) = dg_r1 -dr1_r2 -dr1_nh4 -dr1_oh dydx(igly_r2) = dr1_r2 dydx(igly_nh4) = +dr1_nh4 dydx(igly_oh) = dr1_oh dydx(igly_sfc) = dg_sfc END SUBROUTINE glysoa_complex_derivs SUBROUTINE glysoa_complex(dtchem) USE module_data_mosaic_therm, ONLY : jaerosolstate, all_liquid, mixed, & jtotal, jliquid, nbin_a, & area_wet_a, gas, water_a, aer, mc, & ph, a_nh4, c_as, c_an, & igly, iho, & iglysoa_r1_a, iglysoa_r2_a, & iglysoa_oh_a, & iglysoa_nh4_a, iglysoa_sfc_a, & iso4_a, ino3_a, inh4_a, jc_h ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> USE module_data_mosaic_therm, ONLY: conv1a, & ! converts q/mol(air) to nq/m^3 (q = mol or g) p_atm ! pressure (atm) ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REAL(kind=8), INTENT(IN) :: dtchem REAL(kind=8) :: A, conv, y(nspecs), yout(nspecs), & dydx(nspecs), gly_g INTEGER :: i, ii, nbin_proc, bin_proc(nbin_a) ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> REAL(kind=8), PARAMETER :: Kh_water = 4.19e5 ! effective Henry's law constant of glyoxal in pure water (M atm-1) (Ip et al., GRL, 2011) REAL(kind=8) :: gly_g_atm, & ! gas-phase concentration in atm f_A1, & ! fraction of glyoxal in reservoir 1 c_tot, & ! total concentration of dissolved salts (M) Kh_eq, & ! Henry's law constant at equilibrium (M atm-1) gly_ptot_eq,& ! total glyoxal concentration (reservoirs 1 and 2) at equilibrium gly_r1_eq, & ! glyoxal concentration (reservoir 1) at equilibrium deltagly ! delta to bring r1 in equilibrium ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! liquid / mixed phase bins available? ! --------------------------------------------------------------- nbin_proc = 0 bin_proc(:) = -1 ! see which bins are either mixed, ! or all_liquid (only these are to be processed) DO i = 1, nbin_a IF (jaerosolstate(i) == all_liquid .OR. & jaerosolstate(i) == mixed) THEN nbin_proc = nbin_proc + 1 bin_proc(nbin_proc) = i ENDIF ENDDO IF (nbin_proc == 0) RETURN ! aerosol surface area available? ! --------------------------------------------------------------- A = 0.0 ! total aerosol surface area (cm^2 / cm^3) DO i = 1, nbin_proc ii = bin_proc(i) A = A + area_wet_a(ii) ENDDO IF (A <= 0) RETURN ! clean diagnostic arrays ! --------------------------------------------------------------- ph(:) = -9999.0 ! aerosol pH a_nh4(:) = 0.0 ! ammonium ion activity (M, mol/m^3) c_as(:) = 0.0 ! ammonium sulfate concentration (M, mol/kg) c_an(:) = 0.0 ! ammonium nitrate concentration (M, mol/kg) ! get gas-phase ! --------------------------------------------------------------- ! gly_g will be re-used for all bins gly_g = gas(igly) ! nmol / m3 DO i = 1, nbin_proc ii = bin_proc(i) ! load concentrations array ! ------------------------------------------------------------- conv = 1.e-9 / water_a(ii) ! nmol/m^3 (air) -> mol/kg (water) y(:) = 0.0 ! nmol/m^3 y(igly_g) = gly_g y(igly_r1) = aer(iglysoa_r1_a,jtotal,ii) y(igly_r2) = aer(iglysoa_r2_a,jtotal,ii) y(igly_nh4) = aer(iglysoa_nh4_a,jtotal,ii) y(igly_sfc) = aer(iglysoa_sfc_a,jtotal,ii) y(igly_oh) = aer(iglysoa_oh_a,jtotal,ii) y(ic_as) = aer(iso4_a,jliquid,ii) * conv ! assume we can only form (NH4)2SO4 y(ic_an) = aer(ino3_a,jliquid,ii) * conv ! assume we can only form NH4NO3 y(ia_nh4) = aer(inh4_a,jliquid,ii) * conv ! set activity == concentration y(iph) = MIN(14.0D0, MAX(0.0D0, -log10(mc(jc_h,ii)) )) y(iwater) = water_a(ii) ! kg m-3 y(ioh_g) = gas(iho) ! nmol m-3 y(iarea) = area_wet_a(ii) ! cm^2 / cm^3 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (lfast_tau1) THEN ! instantaneous Henry's law equilibrium ! for fast partitioning (in case characteristic timescale ! of Kampf et al. is only artefact...) gly_g_atm = y(igly_g) / conv1a ! mole / mole gly_g_atm = gly_g_atm * p_atm ! atm ! with A2/A1 = Kolig, and Kolig is 1 f_A1 = 0.5 IF ( y(ic_as) + y(ic_an) .GT. 12.0 ) THEN f_A1 = 0.6667 ENDIF ! Kampff et al., ES&T, 2013, submitted: kinetic limitation of ! salting-in for SO4 concentrations > 12 M. c_tot = MIN( 12.0, y(ic_as) + y(ic_an) ) ! effective Henry's law constant including salting-in effect (Kampff et al.) ! at equilibrium ! derived from eqn. 3 in Kampff et al., -0.24 is "salting-in" constant Kh_eq = Kh_water / 10**(-0.24D0 * c_tot) ! mol L-1 atm-1 == mol kg-1 atm-1 ! gly_g_atm in atm, Kh_eq in mol kg-1 atm-1, water in kg m-3 air ! total glyoxal concentration (reservoirs 1 and 2) at equilibrium gly_ptot_eq = gly_g_atm * Kh_eq * y(iwater) * 1e9 ! nmol / m3 air gly_r1_eq = gly_ptot_eq * f_A1 ! in reservoir 1 at equilibrium deltagly = gly_r1_eq - y(igly_r1) y(igly_g) = y(igly_g) - deltagly y(igly_r1) = y(igly_r1) + deltagly IF (y(igly_g) < 0.0 .OR. y(igly_r1) < 0.0) THEN WRITE(*,*) "THIS IS NOT RIGHT: ",y(igly_g), y(igly_r1),deltagly ENDIF ENDIF ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! integrate system ! ------------------------------------------------------------- CALL glysoa_complex_derivs(1._8, y, dtchem, dydx) CALL rk4(y, dydx, nspecs, 1._8, dtchem, yout, glysoa_complex_derivs) ! update transported fields aer(iglysoa_r1_a,jtotal,ii) = yout(igly_r1) aer(iglysoa_r2_a,jtotal,ii) = yout(igly_r2) aer(iglysoa_nh4_a,jtotal,ii) = yout(igly_nh4) aer(iglysoa_sfc_a,jtotal,ii) = yout(igly_sfc) aer(iglysoa_oh_a,jtotal,ii) = yout(igly_oh) ! do not put gas-phase glyoxal back yet, as we will use it ! for the next bin as well... gly_g = yout(igly_g) ! save diagnostics c_as(ii) = yout(ic_as) c_an(ii) = yout(ic_an) ph(ii) = yout(iph) a_nh4(ii) = yout(ia_nh4) ENDDO ! update gas-phase reservoir, after all bins have been treated. gas(igly) = gly_g END SUBROUTINE glysoa_complex END MODULE module_mosaic_gly