subroutine uEMEP_nonlocal_NO2_O3(nox_bg, nox_emep, no2_emep, o3_emep, J_photo, temperature, f_no2, no2_out, o3_out)
real, intent(in) :: J_photo, temperature, f_no2
real, intent(in) :: nox_bg
real, intent(in) :: nox_emep, no2_emep, o3_emep
real, intent(out) :: no2_out, o3_out
! Local variables
real :: mol_nox_bg
real :: mol_no2_out, mol_o3_out
real :: mol_nox_emep, mol_no2_emep, mol_o3_emep, mol_ox_emep, mol_no_emep
real :: b, d, r, c
real :: Na, Na_fac, k1
real :: p_phot, r_phot
integer no2_i, no_i, nox_i, o3_i, ox_i
integer, parameter :: n_i = 5
real :: mmass(n_i) = [46.0, 30.0, 46.0, 48.0, 47.0]
no2_i = 1; no_i = 2; nox_i = 3; o3_i = 4; ox_i = 5
Na = 6.022e23 ! (molecules/mol)
Na_fac = Na/1.0e12 ! Conversion from ug/m3 to molecules/cm3 included
k1 = 1.4e-12*exp(-1310.0/temperature) ! (cm^3/s) and temperature in Kelvin
! Normally multiplied by *Na_fac but not necessary as it is just a scaling
mol_o3_emep = o3_emep/mmass(o3_i)*Na_fac
mol_no2_emep = no2_emep/mmass(no2_i)*Na_fac
mol_nox_emep = nox_emep/mmass(nox_i)*Na_fac
mol_ox_emep = mol_o3_emep + mol_no2_emep - f_no2*mol_nox_emep
mol_nox_bg = nox_bg/mmass(nox_i)*Na_fac
mol_no_emep = max(0.0, mol_nox_emep - mol_no2_emep)
if (mol_no2_emep .gt. 0) then
r = mol_o3_emep*mol_no_emep/mol_no2_emep
p_phot = J_photo/k1*mol_no2_emep/mol_o3_emep/mol_no_emep
r_phot = J_photo/k1/Na_fac
r = r_phot*Na_fac
b = mol_ox_emep + mol_nox_bg + r
c = max(0.0, b**2 - 4.0*mol_ox_emep*mol_nox_bg) ! Should never be less than 0 but can be -0.0
d = sqrt(c)
mol_no2_out = (b - d)/2.0
mol_o3_out = mol_ox_emep - mol_no2_out
no2_out = max(0.0, mol_no2_out*mmass(no2_i)/Na_fac)
o3_out = max(0.0, mol_o3_out*mmass(o3_i)/Na_fac)
else
no2_out = 0.0
o3_out = o3_emep
end if
end subroutine uEMEP_nonlocal_NO2_O3