[aerodyne-pam-users] revised PAM_chem version
Aerodyne's PAM Users list
aerodyne-pam-users at aerodyne.com
Thu May 18 11:37:20 EDT 2017
We have revised the model code to fix the following mistake:
Previously, 'k_oh_x' on lines 70, 101, or 132:
Xi = [0] ; k_oh_x = 2.7e-11; (k_oh_x specified by user)
was overwritten by line 322:
k_oh_x = 1.7e-12*exp(940./T); (k_oh_x assumes value of alpha-pinene)
The revised model code prompts the user to specify both 'k_oh_x' and
'k_o3_x' on lines 70, 101, or 132. These values should be retained
throughout the code.
A copy of "PAM_chem_v8_rev20170518" is attached and has also been uploaded
to the manual.
Please contact Bill Brune and me with any questions or issues related to
the model implementation.
Andy
Total Control Panel
https://asp.reflexion.net/history-detail?hID=21305674384
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.aerodyne.com/pipermail/aerodyne-pam-users/attachments/20170518/f8d0808c/attachment-0001.html>
-------------- next part --------------
% PAM_chem_v8.m models the oxidation chemistry in the PAM chamber. It does
% not model the detailed organic chemistry nor does it model the aerosol
% chemistry. Its main use is to check the oxidant levels for different
% scenarios that are being used or considered for use. WHB September 2014
%
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%
% The chemistry in this version includes the following:
% fairly complete HOx-Ox chemistry
% HOx-NOx chemistry
% highly simplified RO2 chemistry, starting with species "X"
% SO2 oxidation chemistry
%
% The chemistry in this version does NOT include the following:
% Halogen chemistry
% NO3 chemistry - will be in the next version
% aerosol formation or aging
%
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%
clear all % gets rid of previous runs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Input initial conditions
%
% Variables that are typically changed to run different scenarios are
% initialized in this section.
%
% Set up the appropriate conditions for your experiment in this section.
% Give the study a name (i.e., 'nox') after the word "case". Then write
% this name in "study = ' '.
%
% some example studies :
% no_nox - VOC chemistry with no NOx.
% nox - add NOx chemistry
% 254nm_only - no 185nm UV; add ozone and water vapor.
%
study = '254nm_only'
%
switch study % enables different scenarios to be saved. To set up a different
% scenario, copy and paste an existing scenario,
% rename the copied case, and then change the variables to
% the desired values.
%
%
case 'no_nox'
time = 120; % total time in seconds
tres_mean = 65; % mean residence time in seconds
dt = 0.0005 ; % time step for calculating changes
p = 974; % pressure in hPa
T = 300; % temperature in Kelvin
M = 2.69e19.*(p/1013).*(273.1./T); % calculate the number density of air
n2 = 0.78*M; %
o2 = 0.21*M;
Flux254 = [2e14; 2.0e15]; % 254nm flux in photons cm^-2 sec^-1.
Flux185 = 0.005*Flux254 ; % 185nm light as a fraction of the 254nm light
Flux365 = 1e14; % black lights or LEDs
ozone = [0]; % in ppmv
wv = [1; 2] ; % water vapor in percent
sulfur_dioxide = 100; % in ppbv
OHi = 0; % in ppbv
HO2i = 0; % in ppbv
RO2i = 0; % in ppbv
NOi = 0; % in ppbv
NO2i = 0; % in ppbv
N2Oi = 0; % in ppbv
H2O2i = 0 ; % in ppbv
carbon_monoxide = 0 ; % in ppbv
Xi = [0; 100; 500] ; k_oh_x = 2.7e-11; k_o3_x = 5e-17 ; % X is some species of interest in ppbv
oh_source = 0; % add a separate constant production to doh
ho2_source = 0; % add a separate constant production to dho2
kw_oh = 0; % OH wall loss in units = 1/s
kw_o3 = 0; % O3 wall loss in units = 1/s
ro2_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 8e-13 cm^3 s^-1
sci_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 6e-13 cm^3 s^-1
%
case 'nox'
time = 120; % total time in seconds
tres_mean = 65; % mean residence time in seconds
dt = 0.0005 ; % time step for calculating changes
p = 974; % pressure in hPa
T = 300; % temperature in Kelvin
M = 2.69e19.*(p/1013).*(273.1./T); % calculate the number density of air
n2 = 0.78*M; %
o2 = 0.21*M;
Flux254 = [1e14; 5e14; 2.0e15; 5e15; 1e16 ]; % 254nm flux in photons cm^-2 sec^-1.
Flux185 = 0.005*Flux254 ; % 185nm light as a fraction of the 254nm light
Flux365 = 0; % black lights or LEDs
ozone = [0]; % in ppmv
wv = [1] ; % water vapor in percent
sulfur_dioxide = 100; % in ppbv
OHi = 0; % in ppbv
HO2i = 0; % in ppbv
RO2i = 0; % in ppbv
NOi = 100; % in ppbv
NO2i = 0; % in ppbv
N2Oi = 0; % in ppbv
H2O2i = 0 ; % in ppbv
carbon_monoxide = 0 ; % in ppbv
Xi = [1; 10; 20; 50; 100; 200; 300; 400; 600; 800; 1000] ; k_oh_x = 2.7e-11; k_o3_x = 5e-17; % X is some species of interest in ppbv
oh_source = 0; % add a separate constant production to doh
ho2_source = 0; % add a separate constant production to dho2
kw_oh = 0; % OH wall loss in units = 1/s
kw_o3 = 0; % O3 wall loss in units = 1/s
ro2_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 8e-13 cm^3 s^-1
sci_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 6e-13 cm^3 s^-1
%
case '254nm_only'
time = 120; % total time in seconds
tres_mean = 65; % mean residence time in seconds
dt = 0.0005 ; % time step for calculating changes
p = 974; % pressure in hPa
T = 300; % temperature in Kelvin
M = 2.69e19.*(p/1013).*(273.1./T); % calculate the number density of air
n2 = 0.78*M; %
o2 = 0.21*M;
Flux254 = [1e14; 5e14; 2.0e15; 5e15; 1e16 ]; % 254nm flux in photons cm^-2 sec^-1.
Flux185 = 0*Flux254 ; % 185nm light as a fraction of the 254nm light
Flux365 = 0; % black lights or LEDs
ozone = [6]; % in ppmv
wv = [1] ; % water vapor in percent
sulfur_dioxide = 100; % in ppbv
OHi = 0; % in ppbv
HO2i = 0; % in ppbv
RO2i = 0; % in ppbv
NOi = 100; % in ppbv
NO2i = 0; % in ppbv
N2Oi = 0; % in ppbv
H2O2i = 0 ; % in ppbv
carbon_monoxide = 0 ; % in ppbv
Xi = [1; 10; 20; 50; 100; 200; 300; 400; 600; 800; 1000] ; k_oh_x = 2.7e-11; k_o3_x = 5e-17; % X is some species of interest in ppbv
oh_source = 0; % add a separate constant production to doh
ho2_source = 0; % add a separate constant production to dho2
kw_oh = 0; % OH wall loss in units = 1/s
kw_o3 = 0; % O3 wall loss in units = 1/s
ro2_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 8e-13 cm^3 s^-1
sci_so2_rxn_multiplier = 0; % 0 means no reaction; 1 means 6e-13 cm^3 s^-1
%
end
%
% The sampled air has a distribution of residence times. The parameters for the
% distribution function are as follows using the ideas in Lambe et al. AMT, 2011
% on the PAM and TPOT distributions.
%
pdf_disp(1) = 0;
D1 = 0.06;
tres1 = tres_mean; % From Lambe et al., tres1 = 65 seconds
D2 = 0.21;
tres2 = 1.5*tres_mean; % From Lambe et al., tres2 = 93 seconds
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The concentration calculations start in this section. This model is a
% simple time-step model; it is possible to do better integrations,
% but it's not clear that it is necessary.
%
% Pay attention to the time step "dt". Choose a "dt" to make sure that the
% model results are not dependent on "dt".
%
num_calc = ceil(time/dt); % maximum index number, which is just the time divivded
% by the size of the time steps
num_mean = tres_mean/dt; % the index number for the mean time
%
% convert mixing ratios to concentrations. upper case is used for mixing
% ratios; lower case is used for concentration (molecules cm^{-3})
so2i = sulfur_dioxide*1e-9.*M;
ohi = OHi*1e-12.*M;
ho2i = HO2i*1e-12.*M;
ro2i = RO2i*1e-12.*M;
noi = NOi*1e-9.*M;
no2i = NO2i*1e-9.*M;
n2oi = N2Oi*1e-9.*M;
h2o2i = H2O2i*1e-9.*M;
coi = carbon_monoxide*1e-9.*M;
xi = Xi.*1e-9.*M;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate reaction rate coefficients
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% bimolecular rate coefficients
%
% O(1D) + H2O -> 2 OH
k_o1d_h2o=1.63e-10.*exp(60./T);
% O(1D) + N2 -> O(3P)
k_o1d_n2=2.15e-11.*exp(110./T);
% O(1D) + O2 -> O(3P)
k_o1d_o2=3.3e-11.*exp(55./T);
% O(1D) + N2O -> 2 NO
k_o1d_n2o_1=6.7e-11*exp(20./T);
% O(1D) + N2O -> N2 + O2
k_o1d_n2o_2=4.7e-11*exp(20./T);
% O(1D) + O3 -> O2 + 2 O
k_o1d_o3_1 = 1.2e-10;
% O(1D) + O3 -> 2O2
k_o1d_o3_2 = 1.2e-10;
% O(1D) + H2 -> OH + H
k_o1d_h2 = 1.2e-10;
% O + OH -> O2 + H
k_o_oh = 2.2e-11*exp(120./T);
% O + HO2 -> OH + O2
k_o_ho2 = 3.02e-11*exp(200./T);
% O + H2O2 -> OH + HO2
k_o_h2o2 = 1.4e-12*exp(-2000./T);
% O + O3 -> 2O2
k_o_o3 =8.0e-12*exp(-2060./T);
% NO2 + O -> NO + O2
k_o_no2 = 5.1e-12*exp(210./T);
% H + O3 -> OH + O2
k_h_o3 =1.4e-10*exp(-470./T);
% OH + O3 -> HO2 + O2
k_oh_o3 = 1.7e-12*exp(-940./T);
%k_oh_o3 = 1.28e-13 ;
% OH + HO2 -> H2O + O2
k_oh_ho2 = 4.8e-11*exp(250./T);
% OH + HONO -> H2O + NO2
k_oh_hono = 1.8e-11*exp(-390./T);
% OH + H2O2 -> H2O + HO2
k_oh_h2o2 = 2.9e-12*exp(-160./T);
% OH + H2 -> H2O + H
k_oh_h2 = 2.8e-12*exp(-1800./T);
% OH + OH -> H2O + O
k_oh_oh = 1.8e-12;
% HO2 + O3 -> OH + O2
k_ho2_o3 = 1.0e-14*exp(-490./T);
% HO2 + H -> 2 OH
k_ho2_h_1 = 7.2e-11 ;
% HO2 + H -> O + H2O
k_ho2_h_2 = 1.6e-12 ;
% HO2 + H -> O2 + H2
k_ho2_h_3 = 6.9e-12 ;
% HO2 + NO -> OH + NO2
k_ho2_no = 3.5e-12*exp(270./T);
% NO + O3 --> NO2 + O2
k_no_o3 = 2.0e-12.*exp(-1400./T);
% NO2 + O3 -> NO3 + O2
k_no2_o3 = 1.2e-13.*exp(-2450./T);
%-----------------------------------------------------------------
% termolecular reactions
%
% O + O2 + M -> O3 + M
% k12o and k12h are the low- and high-pressure limit rate constants
k_o_o2_m = 6.0e-34*M.*(300./T).^2.4; % from JPL evaluation 13
%
% H + O2 + M -> HO2 + M
% ko and kh are the low- and high-pressure limit rate constants
ko=4.4e-32*M.*(300./T).^1.3; % from JPL evaluation 13
kh=7.5e-11*(300./T).^0.2;
k_h_o2_m =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
% OH + OH + M -> H2O2 + M
% ko and kh are the low- and high-pressure limit rate constants
ko=6.9e-31*M.*(300./T).^1.0; % from JPL evaluation 13
kh=2.6e-11*(300./T).^0.0;
k_oh_oh_m =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
% OH + NO + M -> HONO + M
% ko and kh are the low- and high-pressure limit rate constants
ko=7.0e-31*M.*(300./T).^2.6; % from JPL evaluation 15
kh=3.6e-11*(300./T).^0.1;
k_oh_no =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
% OH + NO2 + M -> HNO3 + M
% ko and kh are the low- and high-pressure limit rate constants
ko=1.8e-30*M.*(300./T).^3.0; % from JPL evaluation 13
kh=2.8e-11*(300./T).^0.0;
k_oh_no2 =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
% OH + HNO3 --> H2O + NO3
%
k00 = 2.4e-14*exp(460./T);
k01 = 6.5e-34*exp(1335./T);
k02 = 2.7e-17*exp(2199./T);
k_oh_hno3 = k00 + (k01.*M)./(1 + (k01.*M)./k02);
%
% HO2 + NO2 + M -> HO2NO2 + M
% k10o is the low-pressure limit rate constant
% k10h is the high-pressure limit rate constant
% k10r is the reverse thermal decomposition rate
%
ko=1.8e-31*M.*(300./T).^3.2;
kh=4.7e-12*(300./T).^1.4;
k_ho2_no2=(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
k_ho2_no2_r = k_ho2_no2./(2.1e-27*exp(10900./T));
%
% OH + HNO4 --> products
%
k_oh_hno4 = 1.3e-12.*exp(250./T) ;
%
% OH + SO2 + M -> SO3 + M
% ko and kh are the low- and high-pressure limit rate constants
ko=3.3e-31*M.*(300./T).^4.3; % from JPL evaluation 15
kh=1.6e-12*(300./T).^0.0;
k_oh_so2 =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
% OH + CO + M -> H + CO2 + M
% OH + CO --> H + CO2
% ko and kh are the low- and high-pressure limit rate constants
ko=1.5e-13.*(300./T).^(-0.6); % from JPL evaluation 15
kh=(2.1e9*(300./T).^(-6.1))./M;
k_oh_co_1 =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
% OH + CO --> HOCO
% ko and kh are the low- and high-pressure limit rate constants
ko=5.9e-33*M.*(300./T).^1.4; % from JPL evaluation 15
kh=1.1e-12*(300./T).^(-1.3);
k_oh_co_2 =(ko./(1+(ko./kh))).* 0.6.^((1+(log10(ko./kh)).^2).^(-1));
%
k_oh_co = k_oh_co_1 + k_oh_co_2;
%
%-----------------------
% simplified R chemistry, with peroxide and nitrate products.
% Species X is the VOC that initiates the RO2.
%
% typical rate coeeficient values obtained from looking at IUPAC and then
% using the typical value. Rough range of vales are given in comments.
%
% sCI + SO2 -> products
k_sci_so2 = 6e-13.*sci_so2_rxn_multiplier;
% sCI + H2O -> products
k_sci_h2o = 4e-15;
% RO2 + NO -> RO + NO2
k_ro2_no = 8e-12; %(6-10)e-12
% RO2 + HO2 -> ROOH + O2
k_ro2_ho2 = 1.2e-11; %(0.5-1.3)e-11
% ROOH + OH -> RO2 + H2O
k_rooh_oh_1 = 5.3e-12*exp(190./T)*0.6;
% ROOH + OH -> RPHO + OH + H2O
k_rooh_oh_2 = 5.3e-12*exp(190./T)*0.4;
% RO2 + OH -> RPO2 + H2O
k_ro2_oh = 2.3e-10; % Rb is a second organic radical
% RO2 + RO2 -> RP
k_ro2_ro2 = 5e-12; %(0.012-15)e-12
% RO + O2 -> RPO + HO2
k_ro_o2 = 6e-15;
% RO2 + NO + M --> RNO3 + M set equal to 0.02 of k_ro2_no bimolecular reaction (from IUPAC)
k_ro2_no_m = 0.02*k_ro2_no;
% RO + NO + M -> RONO + M
k_ro_no = 3e-11;
% RO + NO2 + M -> RONO2 + M
k_ro_no2 = 3e-11;
% RO2 + SO2 --> products
k_ro2_so2 = 1.4e-14*ro2_so2_rxn_multiplier;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Three nested loops for water vapor, Flux 254nm, and carbon monoxide
% The rate equations are time-stepped in this section.
%
% The program can be run in loops so that results can easily be compared
% amonst different scenarios of CO, lamp flux, and water vapor.
%
for m = 1:length(wv)
%
for k = 1:length(Flux254)
%
for j = 1:length(Xi)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the photolysis rates - some will be recalculated once water
% vapor is known for each time step
%
flux185(1) = Flux185(k);
flux254(1) = Flux254(k);
flux365(1) = Flux365;
%
JO2(1) = 1.1e-20*flux185(1);
JO1D(1) = 6.2e-19*flux185(1) +1.03e-17*flux254(1); % photolysis frequency in s^{-1)
JH2O = 6.78e-21*flux185(1); %H2O+hv-->H+OH
JH2O2_HO2 = 1e-19*flux185(1) ; %HOOH+hv-->HO2+H
JH2O2_OH = 6.7e-20*flux254(1); %HOOH+hv-->2OH
JHO2 = 3.68e-18*flux185(1) + 2.63e-19*flux254(1); %HO2+hv -> OH+O
JNO2 = 1.0e-20*flux185(1) + 1.0e-20*flux254(1) + 5.8e-19*flux365(1) ; %NO2+hv-->O+NO
JHONO = 8.5e-19*flux185(1) + 1.4e-19*flux254(1); %HONO+hv-->OH+NO
JHNO3_1 = 0.6*1.7e-17*flux185(1) + 1.95e-20*flux254(1); %HNO3+hv-->OH+NO2
JHNO3_2 = 0.4*1.7e-17*flux185(1); %HNO3+hv-->O+HONO
JHNO4 = 1.2e-17*flux185(1) + 3.6e-19*flux254(1); %HNO4+hv-->HO2+NO2
JN2O = 1.43e-19*flux185(1); %N2O+hv-->N2+O(1D)
%
% intialize time and the concentrations of gases, for which the index i = 1.
%
%
o1d_ss(1)=0;
o_ss(1) = 0;
h_ss(1) = 0;
o3(1) = ozone*M/1e6;
h2o(1) = wv(m)*M/100;
h2(1) = 0.5e-6.*M ;
ho2(1) = ho2i;
oh(1) = ohi;
h2o2(1) = h2o2i;
so2(1) = so2i;
co(1) = coi;
no(1) = noi;
no2(1) = no2i;
hono(1) = 0;
hno3(1) = 0;
hno4(1) = 0;
n2o(1) = n2oi;
x(1) = xi(j) ;
ro_ss(1) = 0;
ro2(1) = ro2i;
rooh(1) = 0;
rpo(1) = 0;
rp(1) = 0;
rno3(1) = 0;
rono(1) = 0;
rono2(1) = 0;
sci(1) = 0;
t(1) = 0;
%
% "ohr" is the OH reactivity. 'ohr-int' is the OH reactivity due to
% reactions that are internal to the HOx chemistry, such as OH+O3, OH+HO2,
% ... . 'ohr_ext' is the OH reactivity due to the addition of other
% chemical species, such as CO, VOCs, NOx, or SO2.
% Please see Peng et al., AMT, 8, 2015 for more detail.
%
ohr_int(1) = k_o_oh.*o_ss(1) + k_oh_o3.*o3(1) + k_oh_ho2.*ho2(1) + k_oh_hono.*hono(1) + ...
k_oh_h2o2.*h2o2(1) + k_oh_h2.*h2(1) + 2*k_oh_oh.*oh(1) + 2*k_oh_oh_m.*oh(1);
%
ohr_ext(1) = k_oh_no.*no(1) + k_oh_no2.*no2(1) + k_oh_hno3.*hno3(1) + ...
k_oh_x.*x(1) + k_oh_so2.*so2(1) + k_oh_co.*co(1);
%
ohr(1) = ohr_int(1) + ohr_ext(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the variables for all time steps, starting at the second step.
%
for i = 2:num_calc
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% adjust the lamps flux for water vapor, oxygen, and ozone absorption (from
% University of Colorado group)
%
flux185(i) = flux185(1).*exp(-7*(1.1e-20*o2+6.75e-20*h2o(i-1)));
flux254(i) = flux254(1).*exp(-7*(1.15e-17*o3(i-1)));
%
%
% J-values must be recalculated for every loop of the j index because the
% amount of ozone can change.
%
JO2(1) = 1.1e-20*flux185(1);
JO1D(1) = 6.2e-19*flux185(1) +1.03e-17*flux254(1); % photolysis frequency in s^{-1)
JH2O = 6.78e-21*flux185(1); %H2O+hv-->H+OH
JH2O2_HO2 = 1e-19*flux185(1) ; %HOOH+hv-->HO2+H
JH2O2_OH = 6.7e-20*flux254(1); %HOOH+hv-->2OH
JHO2 = 3.68e-18*flux185(1) + 2.63e-19*flux254(1); %HO2+hv -> OH+O
JNO2 = 1.0e-20*flux185(1) + 1.0e-20*flux254(1) + 5.8e-19*flux365(1) ; %NO2+hv-->O+NO
JHONO = 8.5e-19*flux185(1) + 1.4e-19*flux254(1); %HONO+hv-->OH+NO
JHNO3_1 = 0.6*1.7e-17*flux185(1) + 1.95e-20*flux254(1); %HNO3+hv-->OH+NO2
JHNO3_2 = 0.4*1.7e-17*flux185(1); %HNO3+hv-->O+HONO
JHNO4 = 1.2e-17*flux185(1) + 3.6e-19*flux254(1); %HNO4+hv-->HO2+NO2
JN2O = 1.43e-19*flux185(1); %N2O+hv-->N2+O(1D)
%
% Calculating the rate cefficient for HO2+HO2 requires water vapor,
% which changes with water vapor. All other termolecular reactions
% depend only on T or p.
%
% HO2 + HO2 -> H2O2 + O2
k_ho2_ho2=(2.3e-13*exp(600./T)+ 1.7e-33.*M.*exp(1000./T)).*...
(1+1.4e-21.*h2o(i-1).*exp(2200./T));
% k_ho2_ho2 = 2.08e-12
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculate steady-state value for O(1D), O, H, and RO2. Their
% lifetimes are very short compared to the other speices and so we need to
% calculate their steady-state values.
% Steady-state values are the production divided by the loss frequencies.
%
o1d_ss(i) = (JO1D.*o3(i-1) + JN2O.*n2o(i-1))./ ...
(k_o1d_n2.*n2 + k_o1d_o2.*o2 + k_o1d_h2o.*h2o(i-1) + ...
k_o1d_n2o_1.*n2o(i-1) + k_o1d_n2o_2.*n2o(i-1) + ...
k_o1d_o3_1.*o3(i-1) + k_o1d_o3_2.*o3(i-1) + k_o1d_h2.*h2(i-1));
%
o_ss(i) = (k_o1d_n2.*n2.*o1d_ss(i-1) + k_o1d_o2.*o2.*o1d_ss(i-1) + ...
2*k_o1d_o3_1.*o3(i-1).*o1d_ss(i-1) + k_oh_oh.*oh(i-1).*oh(i-1) + ...
JN2O.*n2o(i-1) + 2*JO2.*o2 + JHO2.*ho2(i-1) + JNO2.*no2(i-1) + ...
JHNO3_2.*hno3(i-1) + k_ho2_h_2.*h_ss(i-1).*ho2(i-1))./ ...
(k_o_o2_m.*o2 + k_o_oh.*oh(i-1) + k_o_ho2.*ho2(i-1) + ...
k_o_h2o2.*h2o2(i-1) + k_o_o3.*o3(i-1) + k_o_no2.*no2(i-1));
%
h_ss(i) = (JH2O2_HO2.*h2o2(i-1) + k_o1d_h2.*h2(i-1).*o1d_ss(i-1) + ...
k_o_oh.*o_ss(i-1).*oh(i-1) + k_oh_h2.*oh(i-1).*h2(i-1))./...
(k_h_o3.*o3(i-1) + k_h_o2_m.*o2 + ...
(k_ho2_h_1+k_ho2_h_2+k_ho2_h_3).*ho2(i-1)) ;
%
ro_ss(i) = (k_ro2_no.*no(i-1).*ro2(i-1))./ ...
(k_ro_o2.*o2 + k_ro_no.*no(i-1) + k_ro_no2.*no2(i-1));
%
% Do not consider NO3 and HNO4 chemistry for now. Assume any NO3 formed
% becomes NO2 + O3 due to the photolysis.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up the differential changes of all the chemicals. Each term is
% multiplied by the time step, dt, and then is the production of that
% chemical followed by minus the loss terms that are in parentheses.
%
do3(i) = dt*((k_o_o2_m.*o2.*o_ss(i) + k_oh_hno3.*oh(i-1).*hno3(i-1))...
-(JO1D.*o3(i-1) + k_o_o3.*o_ss(i-1).*o3(i-1) + ...
k_oh_o3.*oh(i-1).*o3(i-1) + k_ho2_o3.*ho2(i-1).*o3(i-1) + ...
k_no_o3.*no(i-1).*o3(i-1) + k_o1d_o3_1.*o1d_ss(i-1).*o3(i-1) + ...
k_o1d_o3_2.*o1d_ss(i-1).*o3(i-1) + kw_o3.*o3(i-1)));
%
doh(i) = dt*((oh_source + JH2O.*h2o(i-1) + 2*k_o1d_h2o.*h2o(i-1).*o1d_ss(i-1) + ...
2*JH2O2_OH.*h2o2(i-1) + JHO2.*ho2(i-1) + JHNO3_1.*hno3(i-1) + ...
JHONO.*hono(i-1) + k_o_ho2.*o_ss(i-1).*ho2(i-1) + ...
k_o_h2o2.*o_ss(i-1).*h2o2(i-1) + k_ho2_o3.*ho2(i-1).*o3(i-1) + ...
k_ho2_no.*ho2(i-1).*no(i-1) + k_o1d_h2.*h2(i-1).*o1d_ss(i-1) + ...
k_h_o3.*h_ss(i-1).*o3(i-1) + 2*k_ho2_h_1.*h_ss(i-1).*ho2(i-1)) ...
- (k_o_oh.*o_ss(i-1).*oh(i-1) + ...
k_oh_o3.*o3(i-1).*oh(i-1) + k_oh_ho2.*ho2(i-1).*oh(i-1) + ...
k_oh_h2o2.*h2o2(i-1).*oh(i-1) + 2*k_oh_oh_m.*oh(i-1).*oh(i-1) + ...
k_oh_no.*no(i-1).*oh(i-1) + k_oh_no2.*no2(i-1).*oh(i-1) + ...
k_oh_hno3.*hno3(i-1).*oh(i-1) + k_oh_hono.*hono(i-1).*oh(i-1) +...
k_oh_so2.*so2(i-1).*oh(i-1) + k_oh_co.*co(i-1).*oh(i-1)+ ...
k_oh_h2.*h2(i-1).*oh(i-1) + k_oh_x.*x(i-1).*oh(i-1) + ...
k_ro2_oh.*ro2(i-1).*oh(i-1) + (k_rooh_oh_1 + k_rooh_oh_2).*rooh(i-1).*oh(i-1) +...
2*k_oh_oh.*oh(i-1).*oh(i-1) + kw_oh.*oh(i-1)));
%
dho2(i) = dt*((ho2_source + JH2O.*h2o(i-1) + JH2O2_HO2.*h2o2(i-1) + k_oh_o3.*o3(i-1).*oh(i-1) + ...
k_h_o2_m.*h_ss(i-1).*o2 + k_o_h2o2.*o_ss(i-1).*h2o2(i-1) + ...
k_oh_h2o2.*h2o2(i-1).*oh(i-1) + k_oh_so2.*so2(i-1).*oh(i-1) + ...
k_oh_co.*co(i-1).*oh(i-1) + k_ro_o2*ro_ss(i-1).*o2 + ...
k_o1d_h2.*h2(i-1).*o1d_ss(i-1)) ...
- (JHO2.*ho2(i-1) + k_o_ho2.*o_ss(i-1).*ho2(i-1) + ...
k_ho2_o3.*o3(i-1).*ho2(i-1) + k_oh_ho2.*oh(i-1).*ho2(i-1) + ...
2*k_ho2_ho2.*ho2(i-1).*ho2(i-1) + k_ho2_no.*no(i-1).*ho2(i-1) + ...
k_ro2_ho2.*ro2(i-1).*ho2(i-1)) + ...
(k_ho2_h_1+k_ho2_h_2+k_ho2_h_3).*h_ss(i-1).*ho2(i-1)) ;
%
dh2o2(i) = dt*((k_ho2_ho2.*ho2(i-1).*ho2(i-1)+ k_oh_oh_m.*oh(i-1).*oh(i-1)) ...
- (k_o_h2o2.*o_ss(i-1).*h2o2(i-1) + k_oh_h2o2.*oh(i-1).*h2o2(i-1) + ...
(JH2O2_OH + JH2O2_HO2).*h2o2(i-1)));
%
dh2o(i) = dt*((k_oh_ho2.*oh(i-1).*ho2(i-1) + k_oh_h2o2.*h2o2(i-1).*oh(i-1) + ...
k_o_h2o2.*h2o2(i-1).*o_ss(i-1) + k_oh_hono.*hono(i-1).*oh(i-1) + ...
k_oh_hno3.*oh(i-1).*hno3(i-1) + k_oh_x.*x(i-1).*oh(i-1) + ...
k_ho2_h_2.*ho2(i-1).*h_ss(i-1) + k_oh_h2.*oh(i-1).*h2(i-1) + ...
k_oh_oh.*oh(i-1).*oh(i-1)) + k_ro2_oh.*ro2(i-1).*oh(i-1) + ...
- (k_o1d_h2o.*h2o(i-1).*o1d_ss(i-1))) ;
%
dh2(i) = dt*((k_ho2_h_3.*ho2(i-1).*h_ss(i-1))...
- (k_o1d_h2.*h2(i-1).*o1d_ss(i-1) - k_oh_h2.*oh(i-1).*h2(i-1)));
%
dno(i) = dt*((JNO2.*no2(i-1) + JHONO.*hono(i-1) + ...
2*k_o1d_n2o_1.*o1d_ss(i-1).*n2o(i-1) + k_o_no2.*no2(i-1).*o_ss(i-1))...
- (k_ho2_no.*ho2(i-1).*no(i-1) + k_no_o3.*o3(i-1).*no(i-1) + ...
k_oh_no.*oh(i-1).*no(i-1) + k_ro2_no.*ro2(i-1).*no(i-1) + ...
k_ro2_no_m.*ro2(i-1).*no(i-1) + k_ro_no.*ro_ss(i-1).*no(i-1)));
%
dno2(i) = dt*((JHNO3_1.*hno3(i-1) + k_ho2_no.*ho2(i-1).*no(i-1) + ...
k_no_o3.*no(i-1).*o3(i-1) + k_oh_hno3.*oh(i-1).*hno3(i-1) + ...
k_oh_hono.*oh(i-1).*hono(i-1) + k_ro2_no.*ro2(i-1).*no(i-1)) ...
- (JNO2.*no2(i-1) + k_oh_no2.*oh(i-1).*no2(i-1) + ...
k_ro_no2.*ro_ss(i-1).*no2(i-1)));
%
dhono(i) = dt*((k_oh_no.*oh(i-1).*no(i-1) + JHNO3_2.*hno3(i-1)) ...
- (JHONO.*hono(i-1) + k_oh_hono.*oh(i-1).*hono(i-1)));
%
dhno3(i) = dt*((k_oh_no2.*oh(i-1).*no2(i-1)) ...
- ((JHNO3_1+JHNO3_2).*hno3(i-1) + k_oh_hno3.*oh(i-1).*hno3(i-1)));
%
dn2o(i) = dt*(-(JN2O.*n2o(i-1) + (k_o1d_n2o_1+k_o1d_n2o_2).*o1d_ss(i-1).*n2o(i-1)));
%
% dhno4 = dt* Assume HNO4 is in steady state.
%
dso2(i) = -dt*(k_oh_so2.*so2(i-1).*oh(i-1) + k_ro2_so2.*so2(i-1).*ro2(i-1));
%
dco(i) = dt*(-k_oh_co.*co(i-1).*oh(i-1));
%
dx(i) = dt*(-k_oh_x.*x(i-1).*oh(i-1));
%
dco(i) = dt*(-(k_oh_co.*co(i-1).*oh(i-1)));
%
dx(i) = dt*(-k_oh_x.*x(i-1).*oh(i-1));
%
dro2(i) = dt*((k_oh_x.*x(i-1).*oh(i-1) + k_rooh_oh_1.*rooh(i-1).*oh(i-1))...
- (k_ro2_no.*no(i-1).*ro2(i-1) + k_ro2_no_m.*no(i-1).*ro2(i-1) + ...
k_ro2_oh.*ro2(i-1).*oh(i-1) + k_ro2_so2.*so2(i-1).*ro2(i-1) + ...
k_ro2_ho2.*ho2(i-1).*ro2(i-1) + k_ro2_ro2.*ro2(i-1).*ro2(i-1)));
%
drp(i) = dt*(k_ro2_ro2.*ro2(i-1).*ro2(i-1));
%
drooh(i) = dt*((k_ro2_ho2.*ho2(i-1).*ro2(i-1)) ...
- ((k_rooh_oh_1 + k_rooh_oh_2).*rooh(i-1).*oh(i-1)));
%
drno3(i) = dt*(k_ro2_no_m.*no(i-1).*ro2(i-1));
%
drono(i) = dt*(k_ro_no.*no(i-1).*ro_ss(i-1));
%
drono2(i) = dt*(k_ro_no2.*no2(i-1).*ro_ss(i-1));
%
drpo(i) = dt*(k_ro_o2.*o2.*ro_ss(i-1));
%
dsci(i) = dt*((0.3.*k_o3_x.*o3(i-1).*x(i-1)) ...
- (k_sci_h2o.*sci(i-1).*h2o(i-1) + k_sci_so2.*sci(i-1).*so2(i-1)));
%
%
% determine new concentration values by adding the calculated changes
% and adding them to the values at a previous time step
%
o3(i) = o3(i-1) + do3(i);
oh(i) = oh(i-1) + doh(i);
ho2(i) = ho2(i-1) + dho2(i);
h2o2(i) = h2o2(i-1) + dh2o2(i);
h2o(i) = h2o(i-1) + dh2o(i) ;
h2(i) = h2(i-1) + dh2(i);
no(i) = no(i-1) + dno(i);
no2(i) = no2(i-1) + dno2(i);
hono(i) = hono(i-1) + dhono(i);
hno3(i) = hno3(i-1) + dhno3(i);
n2o(i) = n2o(i-1) + dn2o(i);
so2(i) = so2(i-1) + dso2(i);
co(i) = co(i-1) + dco(i);
x(i) = x(i-1) + dx(i) ;
ro2(i) = ro2(i-1) + dro2(i);
rooh(i) = rooh(i-1) + drooh(i);
rpo(i) = rpo(i-1) + drpo(i);
rp(i) = rp(i-1) + drp(i);
rno3(i) = rno3(i-1) + drno3(i);
rono(i) = rono(i-1) + drono(i);
rono2(i) = rono2(i-1) + drono2(i);
sci(i) = sci(i-1) + dsci(i);
%
% increase the time by adding a time step
%
t(i) = t(i-1) + dt;
%
% calculate values for the internal OHR (HOx chemistry only) and the
% exterbal OHR (everything but the HOx chemistry
%
ohr_int(i) = k_o_oh.*o_ss(i) + k_oh_o3.*o3(i) + k_oh_ho2.*ho2(i) + k_oh_hono.*hono(i) + ...
k_oh_h2o2.*h2o2(i) + k_oh_h2.*h2(i) + 2*k_oh_oh.*oh(i) + 2*k_oh_oh_m.*oh(i);
%
ohr_ext(i) = k_oh_no.*no(i) + k_oh_no2.*no2(i) + k_oh_hno3.*hno3(i) + ...
k_oh_x.*x(i) + k_oh_so2.*so2(i) + k_oh_co.*co(i) + ...
(k_rooh_oh_1 + k_rooh_oh_2).*rooh(i) + k_ro2_oh.*ro2(i) ;
%
ohr(i) = ohr_int(i) + ohr_ext(i);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculate the distribution function for the sampled air using the
% equations determined in Lambe et al. 2011 (flow tube comparison)
%
tau1 = t(i)/tres1;
tau2 = t(i)/tres2;
pdf_disp(i) = exp(-(1-tau1).^2./(D1.*tau1))./(2*sqrt(3.14156*D1.*tau1))...
+ exp(-(1-tau2).^2./(D2.*tau2))./(2*sqrt(3.14156*D2.*tau2));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
end %for loop with i - time steps
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Convert concentrations to mixing ratios
%
O3 = 1e6.*o3./M; % ozone in ppmv
OH = 1e9.*oh./M; % OH in ppbv
HO2 = 1e9.*ho2./M; % HO2 in ppbv
RO2 = 1e9.*ro2./M; % RO2 in ppbv
H2O2 = 1e9.*h2o2./M; % HOOH in ppbv
ROOH = 1e9.*rooh./M; % ROOH in ppbv
H2O = 100*h2o./M; % H2O in percent
NO = 1e9.*no./M; % NO in ppbv
NO2 = 1e9.*no2./M; % NO2 in ppbv
HONO = 1e9.*hono./M; % HONO in ppbv
HNO3 = 1e9.*hno3./M; % HNO3 in ppbv
N2O = 1e9.*n2o./M; % N2O in ppbv
SO2 = 1e9.*so2./M; % SO2 in ppbv
CO = 1e9.*co./M; % CO in ppbv
X = 1e9.*x./M; % X in ppbv
sCI = 1e12.*sci./M; % sCI in pptv
%
% Calculate the OH exposure
%
oh_exp = cumsum(oh)*dt;
oh_exp_so2 = log(SO2(1)./SO2)./k_oh_so2;
ho2_exp = cumsum(ho2)*dt;
o3_exp = cumsum(o3)*dt;
%
% set up the matrices for this species of interest. These values are the
% values for the index i associated with the mean time and are followed
% with an "f".
%
Flux254f(j,k,m) = flux254(1);
Flux185f(j,k,m) = flux185(1);
H2Os = H2O(1);
H2Of(j,k,m) = H2O(num_mean);
O3s(j,k,m) = O3(1);
O3f(j,k,m) = O3(num_mean);
OHf(j,k,m) = OH(num_mean);
HO2f(j,k,m) = HO2(num_mean);
H2O2f(j,k,m) = H2O2(num_mean);
ROOHf(j,k,m) = ROOH(num_mean);
deltaO3f(j,k,m) = O3(1)-O3(num_mean);
fracO3f(j,k,m) = (O3(1)-O3(num_mean))./O3(1);
NOf(j,k,m) = NO(num_mean);
NO2f(j,k,m) = NO2(num_mean);
HONOf(j,k,m) = HONO(num_mean);
HNO3f(j,k,m) = HNO3(num_mean);
N2Os = N2O(1);
N2Of(j,k,m) = N2O(num_mean);
SO2s(j,k,m) = SO2(1);
SO2f(j,k,m) = SO2(num_mean);
COs(j,k,m) = CO(1);
COf(j,k,m) = CO(num_mean);
deltaCOf(j,k,m) = COs(j,k,m)-COf(j,k,m);
Xs(j,k,m) = X(1);
Xf(j,k,m) = X(num_mean);
oh_expf(j,k,m) = oh_exp(floor(tres_mean./dt));
oh_exp_so2f(j,k,m) = oh_exp_so2(floor(tres_mean./dt));
OHRs_int(j,k,m) = ohr_int(1);
OHRs_ext(j,k,m) = ohr_ext(1);
OHRs(j,k,m) = ohr(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% sum the product of the normalized distribution function for all the
% variables of interest. These will be disignated with the letter "d"
%
pdf_sum = sum(pdf_disp); % needed to normalize the pdf function
%
H2Od(j,k,m) = sum(H2O.*pdf_disp)./pdf_sum;
O3d(j,k,m) = sum(O3.*pdf_disp)./pdf_sum;
deltaO3d(j,k,m) = O3(1)-O3d(j,k,m);
fracO3d(j,k,m) = (O3(1)-O3d(j,k,m))./O3(1);
OHd(j,k,m) = sum(OH.*pdf_disp)./pdf_sum;
HO2d(j,k,m) = sum(HO2.*pdf_disp)./pdf_sum;
H2O2d(j,k,m) = sum(H2O2.*pdf_disp)./pdf_sum;
ROOHd(j,k,m) = sum(ROOH.*pdf_disp)./pdf_sum;
NOd(j,k,m) = sum(NO.*pdf_disp)./pdf_sum;
NO2d(j,k,m) = sum(NO2.*pdf_disp)./pdf_sum;
HONOd(j,k,m) = sum(HONO.*pdf_disp)./pdf_sum;
HNO3d(j,k,m) = sum(HNO3.*pdf_disp)./pdf_sum;
N2Od(j,k,m) = sum(N2O.*pdf_disp)./pdf_sum;
SO2d(j,k,m) = sum(SO2.*pdf_disp)./pdf_sum;
COd(j,k,m) = sum(CO.*pdf_disp)./pdf_sum;
deltaCOd(j,k,m) = COs(j,k,m)-COd(j,k,m);
Xd(j,k,m) = sum(X.*pdf_disp)./pdf_sum;
oh_expd(j,k,m) = sum(oh_exp.*pdf_disp)./pdf_sum;
oh_exp_so2d(j,k,m) = sum(oh_exp_so2.*pdf_disp)./pdf_sum;
OHRd_int(j,k,m) = sum(ohr_int.*pdf_disp)./pdf_sum;
OHRd_ext(j,k,m) = sum(ohr_ext.*pdf_disp)./pdf_sum;
OHRd(j,k,m) = sum(ohr.*pdf_disp)./pdf_sum;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot the results for the major chemical species as a function of time.
% Figure of OH, HO2, O3, and SO2 for individual values of water vapor,
% photolysis, and CO
%
figure
subplot(3,1,1)
semilogy(t,O3,'r-',t,OH,'b-',t,HO2,'g-',t,X,'k-',t,SO2,'c-','linewidth',2)
axis([0 max(t) .1 1e4])
xlabel('time (sec)')
ylabel('O_3ppm OHppb HO_2ppb Xppb SO_2ppb')
legend('O_3','OH','HO_2','X','SO_2','Location','Best')
title(['Flux=',num2str(Flux254(k),'%4.1e'),'ph cm^{-2}s^{-1}; O_3(initial)=',num2str(ozone),'ppbv; H_2O=',num2str(wv(m)),'%; OHexp =',num2str(oh_expf(j,k,m),'%4.2e'),'; OHexpSO2 =',num2str(oh_exp_so2f(j,k,m),'%4.2e\n' )])
grid
%
subplot(3,1,2)
semilogy(t,O3,'r-',t,NO,'k-',t,NO2,'b-',t,HNO3,'g-','linewidth',2)
axis([0 max(t) .1 1e4])
xlabel('time (sec)')
ylabel('O_3ppm NOppb NO_2ppb HNO_3ppb')
legend('O_3','NO','NO_2','HNO_3','Location','Best')
grid
%
subplot(3,1,3)
semilogy(t,O3,'r-',t,ohr_ext,'k-',t,ohr_int,'b-',t,oh_exp/1e10,'c-','linewidth',2)
axis([0 max(t) .1 1e4])
xlabel('time (sec)')
ylabel('O_3ppm OHRext s^{-1} OHRint s^{-1} OHexp/1e8')
legend('O_3','OHR_{ext}','OHR_{int}','OHexp/1e10','Location','Best')
grid
%
%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%
end % for j loop - X (first index)
%
end % for k loop - lamp flux (second index)
%
end % for m loop - water vapor (third index)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot some results as a function of valuse in the nested loops.
%
%
figure;
for m = 1:length(wv)
for k=1:length(Flux254)
scatter(OHRs_ext(:,k,m),(oh_expd(:,k,m)),60,Flux254f(:,k,m),'filled');
hold on
end
end
grid;colorbar;
ylabel('OH exp (molecules cm^{-3} s)'); xlabel('initial external OHR (s^{-1})')
title('PAM chemistry model')
hold off
%
% Print out the OH exposure.
%
fprintf('OH exposure (weighted) = %4.2e cm^{-3} s \r', oh_expd);
fprintf('OH exposure SO2 (weighted) = %4.2e cm^{-3} s \r', oh_exp_so2d);
More information about the Aerodyne-PAM-users
mailing list