subroutine uEMEP_phototimescale_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, time_scale, nox_out, no2_out, o3_out, p_bg_out, p_out)
real, intent(in) :: nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, time_scale
real, intent(out) :: nox_out, no2_out, o3_out, p_bg_out, p_out
! Local variables
integer :: no2_i,no_i,nox_i,o3_i,ox_i,nox_bg_i,no2_bg_i
integer, parameter :: n_i = 7
double precision :: Na, Na_fac, k1
double precision :: mass(n_i)
double precision :: mmass(n_i) = [46.0, 30.0, 46.0, 48.0, 47.0, 46.0, 46.0]
double precision :: mol(n_i)
double precision :: fac_sqrt
double precision :: f_no2,f_ox,Jd,Jd_bg
double precision :: min_nox=1.0e-6
double precision :: c, b, BB, td, f_no2_0, f_no2_ps
complex(8) :: AA
double precision :: p_tot_out, f_ox_bg, f_no2_bg_ps, f_no2_bg
no2_i = 1; no_i = 2; nox_i = 3; o3_i = 4; ox_i = 5; nox_bg_i = 6; no2_bg_i = 7
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
mass(1:n_i) = 0.0
mol(1:n_i) = 0.0
! Test for 0 NOx. If so leave the routine
mass(nox_i) = nox_loc + nox_bg
if (mass(nox_i) .le. min_nox) then
nox_out = 0.0
no2_out = 0.0
o3_out = o3_bg
return
end if
! Check the photostationary assumption for the input data
mass(nox_i) = nox_bg
mass(no2_i) = no2_bg
mass(o3_i) = o3_bg
mol = mass/mmass*Na_fac ! (molecules per cm3)
mol(ox_i) = mol(o3_i) + mol(no2_i)
mol(no_i) = max(0.0, mol(nox_i) - mol(no2_i))
f_ox_bg = mol(ox_i)/mol(nox_i)
Jd_bg = J_photo/k1/mol(nox_i)
f_no2_bg_ps = 0.5*((1 + f_ox_bg + Jd_bg) - sqrt((1 + f_ox_bg+Jd_bg)**2 - 4.0*f_ox_bg))
f_no2_bg = mol(no2_i)/mol(nox_i)
p_bg_out = f_no2_bg/f_no2_bg_ps
! Check input
if (abs(J_photo) > epsilon0 .and. abs(mol(no_i)) > epsilon0 .and. abs(mol(o3_i)) > epsilon0) then
p_bg_out = J_photo*mol(no2_i)/k1/mol(o3_i)/mol(no_i)
else
p_bg_out = mol(no2_i)/(mol(ox_i) + mol(nox_i) - abs(mol(ox_i) - mol(nox_i)))*2.0
end if
! Add the local contribution for calculation
mass(nox_i) = nox_loc + nox_bg
mass(no2_i) = f_no2_loc*nox_loc + no2_bg
mass(o3_i) = o3_bg
mol = mass/mmass*Na_fac ! (molecules per cm3)
mol(ox_i) = mol(o3_i) + mol(no2_i)
mol(no_i) = max(0.0, mol(nox_i) - mol(no2_i))
f_no2 = mol(no2_i)/mol(nox_i)
f_ox = mol(ox_i)/mol(nox_i)
! Set the photolysis rate
Jd = J_photo/k1/mol(nox_i)
! Calculate photostationary for total nox, ox
fac_sqrt = max(0.0, (1 + f_ox + Jd)**2 - 4.0*f_ox)
f_no2_ps = 0.5*((1 + f_ox + Jd) - sqrt(fac_sqrt))
p_tot_out = f_no2/f_no2_ps
! Calculate fraction of NO2 based on the time scale
c = f_ox
b = 1 + f_ox + Jd
BB = sqrt(max(0.0, b**2 - 4.0*c))! max avoids roundoff errors
td = time_scale*k1*mol(nox_i)
f_no2_0 = f_no2
AA = clog(cmplx((BB + b - 2.0*f_no2_0)/(BB - b + 2.0*f_no2_0)))
f_no2 = real(-BB/2.0*((exp(AA + BB*td) - 1.0)/(exp(AA + BB*td) + 1.0)) + b/2.0)
if (isnan(f_no2)) f_no2 = -BB/2.0 + b/2.0
fac_sqrt = max(0.0, (1 + f_ox + Jd)**2 - 4.0*f_ox)
f_no2_ps = 0.5*((1.0 + f_ox + Jd) - sqrt(fac_sqrt))
p_out = f_no2/f_no2_ps
! Convert back to mass
mol(no2_i) = max(0.0, f_no2*mol(nox_i))
mol(o3_i) = max(0.0, mol(ox_i) - mol(no2_i)) ! Rounding errors possible
mol(no_i) = max(0.0, mol(nox_i) - mol(no2_i)) ! Rounding errors possible
mass = mol*mmass/Na_fac ! (ug/m3)
no2_out = mass(no2_i)
nox_out = mass(nox_i)
o3_out = mass(o3_i)
if (isnan(no2_out)) then
write(*,'(8a12)') 'nox_bg', 'no2_bg', 'o3_bg', 'nox_loc', 'f_no2_loc', 'J_photo', 'temperature', 'time_scale'
write(*,'(8es12.2)') nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, time_scale
write(*,'(4a12)') 'f_no2', 'BB', 'b', 'b**2-4.*c'
write(*,'(4es12.2)') f_no2, BB, b, b**2 - 4.0*c
stop 1
end if
! Check output
if (abs(J_photo) > epsilon0 .and. abs(mol(no_i)) > epsilon0 .and. abs(mol(o3_i)) > epsilon0) then
p_out = J_photo*mol(no2_i)/k1/mol(o3_i)/mol(no_i)
else
p_out = mol(no2_i)/(mol(ox_i) + mol(nox_i) - abs(mol(ox_i) - mol(nox_i)))*2.0
end if
end subroutine uEMEP_phototimescale_NO2