subroutine uEMEP_annual_mean_pdf_correction_NO2_O3(bin_min, bin_max, delta_log10_bin, run_all, no2_in, nox_in, o3_in, J_photo_in, temperature_in, ox_sigma_ratio_in, nox_sigma_ratio_in, lon_in, lat_in, no2_out, o3_out)
real, intent(in) :: bin_min, bin_max, delta_log10_bin
logical, intent(in) :: run_all
real, intent(in) :: J_photo_in, temperature_in
real, intent(in) :: no2_in, nox_in, o3_in
real, intent(in) :: ox_sigma_ratio_in, nox_sigma_ratio_in
real, intent(in) :: lon_in, lat_in
real, intent(out) :: no2_out, o3_out
! Loca variables
real :: mol_no2_out, mol_o3_out
real :: Na, Na_fac, k1
integer, parameter :: no2_i = 1
!integer, parameter :: no_i = 2
integer, parameter :: nox_i = 3
integer, parameter :: o3_i = 4
!integer, parameter :: ox_i = 5
integer, parameter :: n_i = 5
real :: mmass(n_i) = [46.0, 30.0, 46.0, 48.0, 47.0]
real :: log10_bin_max, log10_bin_min
real, allocatable :: log10_bin(:), bin(:), delta_bin(:)
integer :: n_bin
real :: nox_sigma_ratio, ox_sigma_ratio
real :: ox_sigma, ox_sig_sqr, ox_mu, nox_sigma, nox_sig_sqr, nox_mu
real, allocatable :: y_nox(:), y_ox(:), y_Jd(:)
integer, parameter :: nbin_Jd = 8760 ! Hours in a year
integer :: bin_temp
real, save :: y_Jd_acc(nbin_Jd)
real :: mean_y_Jd_acc, y_Jd_acc_temp, y_Jd_acc_temp_log10
real :: y_all_val, y_all_prob, y_all_sum, y_all_prob_sum, y_all, y_annual, y_scale
real :: azimuth_ang, zenith_ang
real :: mol_nox, mol_no2, mol_ox, mol_o3, Jd
double precision :: date_num
integer :: i, j, k, l
integer :: date_array(6)
real :: min_sigma_ratio = 0.01
real :: bin_temp2
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_in) ! (cm^3/s) and temperature in Kelvin
mol_nox = nox_in*Na_fac/mmass(nox_i)
mol_no2 = no2_in*Na_fac/mmass(no2_i)
mol_o3 = o3_in*Na_fac/mmass(o3_i)
mol_ox = mol_no2 + mol_o3
Jd = J_photo_in/k1
no2_out = no2_in
o3_out = o3_in
!Create the bins for the pdf in (mol/cm3). ox, nox and J
log10_bin_max = log10(bin_max)
log10_bin_min = log10(bin_min)
n_bin = int((log10_bin_max - log10_bin_min)/delta_log10_bin) + 1
! Creates 80 bins with these settings
allocate (log10_bin(n_bin))
allocate (bin(n_bin))
allocate (delta_bin(n_bin))
do i = 1, n_bin
log10_bin(i) = log10_bin_min + (i - 1)*delta_log10_bin
end do
do i = 1, n_bin
bin(i) = (10.0**log10_bin(i))*Na_fac/mmass(nox_i)
delta_bin(i) = (10.0**(log10_bin(i) + delta_log10_bin/2.0) - 10.0**(log10_bin(i) - delta_log10_bin/2.0))*Na_fac/mmass(nox_i)
end do
! Distribute concentrations into the pdf bins
! Minimum needed to avoid NaNs in the calculation
nox_sigma_ratio = 1.14
ox_sigma_ratio = 0.21
if (abs(nox_sigma_ratio_in) > epsilon0) nox_sigma_ratio = max(nox_sigma_ratio_in, min_sigma_ratio)
if (abs(ox_sigma_ratio_in) > epsilon0) ox_sigma_ratio = max(ox_sigma_ratio_in, min_sigma_ratio)
ox_sigma = mol_ox*ox_sigma_ratio
ox_sig_sqr = log(1.0 + ox_sigma**2.0/mol_ox**2.0)
ox_mu = log(mol_ox**2.0/sqrt(mol_ox**2.0 + ox_sigma**2.0))
nox_sigma = mol_nox * nox_sigma_ratio
nox_sig_sqr = log(1.0 + nox_sigma**2.0/mol_nox**2.0)
nox_mu = log(mol_nox**2.0/sqrt(mol_nox**2.0 + nox_sigma**2.0))
allocate (y_ox(n_bin))
allocate (y_nox(n_bin))
if ( .not. allocated(y_Jd)) allocate (y_Jd(n_bin))
y_Jd = 0.0
y_ox = 0.0
y_nox = 0.0
do i = 1, n_bin
y_ox(i) = 1.0/bin(i)/sqrt(ox_sig_sqr)/sqrt(2.0*pi)*exp(-((log(bin(i)) - ox_mu)**2)/2.0/ox_sig_sqr)*delta_bin(i)
y_nox(i) = 1.0/bin(i)/sqrt(nox_sig_sqr)/sqrt(2.0*pi)*exp(-((log(bin(i)) - nox_mu)**2)/2.0/nox_sig_sqr)*delta_bin(i)
end do
! Create the Jd_acc distribution by looping through every hour in the year and extracting the zenith angle
! Only do this if requested for the first time
if (run_all) then
y_Jd_acc = 0
do i = 1, nbin_Jd
date_num = 1.0 + i/24.0
date_array = 0
zenith_ang = 0.0
call get_sun_angles(lat_in, lon_in, date_array, date_num, 0.0, azimuth_ang, zenith_ang)
if (zenith_ang .ge. 90.0) then
y_Jd_acc(i) = 0
else
y_Jd_acc(i) = ((cosd(zenith_ang))**0.28)
end if
end do
end if
mean_y_Jd_acc = sum(y_Jd_acc)/nbin_Jd
do i = 1, nbin_Jd
y_Jd_acc_temp = Jd*y_Jd_acc(i)/mean_y_Jd_acc
if (y_Jd_acc_temp/Na_fac*mmass(nox_i) < bin_min) then
bin_temp = 1
else
y_Jd_acc_temp_log10 = log10(y_Jd_acc_temp/Na_fac*mmass(nox_i))
bin_temp = floor((y_Jd_acc_temp_log10 - log10_bin_min)/delta_log10_bin + 0.5) + 1
bin_temp = min(max(bin_temp, 1), n_bin)
end if
y_Jd(bin_temp) = y_Jd(bin_temp) + 1
end do
! Normalise all distributions
y_Jd = y_Jd/sum(y_Jd)
y_ox = y_ox/sum(y_ox)
y_nox = y_nox/sum(y_nox)
!Calculate scaling based on photostationary assumption
y_all_sum = 0
y_all_prob_sum = 0
do j = 1, n_bin
do k = 1, n_bin
do l = 1, n_bin
! Calculate weighting
y_all_prob = y_nox(j)*y_Jd(k)*y_ox(l)
if (y_all_prob .gt. 0.0) then
! Calculate NO2 value
bin_temp2 = bin(j) + bin(k) + bin(l)
y_all_val = (bin_temp2 - sqrt(bin_temp2*bin_temp2 - 4.0*bin(j)*bin(l)))/2.0
! Add weighted value
y_all_sum = y_all_sum + y_all_val*y_all_prob
! Calculate sum of weights for normalisation later
y_all_prob_sum = y_all_prob_sum + y_all_prob
end if
end do
end do
end do
y_all = y_all_sum/y_all_prob_sum
! Calculate the mean value
y_annual = ((mol_nox + Jd + mol_ox) - sqrt((mol_nox + Jd + mol_ox)**2 - 4.0*mol_nox*mol_ox))/2.0
y_scale = y_all/y_annual
mol_no2_out = y_scale*mol_no2
mol_o3_out = mol_ox - mol_no2_out
no2_out = mol_no2_out/Na_fac*mmass(no2_i)
o3_out = mol_o3_out/Na_fac*mmass(o3_i)
deallocate (log10_bin)
deallocate (bin)
deallocate (delta_bin)
deallocate (y_ox)
deallocate (y_nox)
deallocate (y_Jd)
end subroutine uEMEP_annual_mean_pdf_correction_NO2_O3