uEMEP_During_NO2 Subroutine

    subroutine uEMEP_During_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_emep, no2_emep, o3_emep, J_photo, temperature, nox_out, no2_out, o3_out, p_bg_out, p_out)
        ! Reference
        ! D�ring, I., B�chlin, W., Ketzel, M., Baum, A., Friedrich, U., Wurzler, S., 2011.
        ! A new simplified NO/NO2 conversion model under consideration of direct NO2-emissions.
        ! Meteorol. Zeitschrift 20, 67�73. doi:10.1127/0941-2948/2011/0491
        !
        ! Improved Methodologies for NO2 Exposure Assessment in the EU, page 53
        ! https://ec.europa.eu/environment/air/pdf/NO2_Exposure_Final_Report.pdf
        real, intent(in) :: nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, J_photo, temperature
        real, intent(in) :: nox_emep, no2_emep, o3_emep
        real, intent(out) :: nox_out, no2_out, o3_out
        real, intent(out) :: p_bg_out, p_out

        ! Local variables
        real :: mol_nox_bg, mol_no2_bg, mol_nox_loc, mol_o3_bg, mol_no2_loc, mol_ox_loc, mol_no_bg, mol_ox_bg
        real :: mol_no2_out, mol_o3_out, mol_no_out
        real :: mol_nox_emep, mol_no2_emep, mol_o3_emep, mol_no_emep, mol_ox_emep, p_emep_out
        real :: b, d, r, c, k1
        real :: Na, Na_fac
        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
        k1 = 1.4e-12*exp(-1310.0/temperature) ! (cm^3/s) and temperature in Kelvin
        Na = 6.022e23 ! (molecules/mol)
        Na_fac = Na/1.0e12 ! Conversion from ug/m3 to molecules/cm3 included
        
        ! Normally multiplied by *Na_fac but not necessary as it is just a scaling
        mol_no2_bg = no2_bg/mmass(no2_i)
        mol_no2_loc = f_no2_loc*nox_loc/mmass(no2_i)
        mol_nox_bg = nox_bg/mmass(nox_i)
        mol_nox_loc = nox_loc/mmass(nox_i)
        mol_o3_bg = o3_bg/mmass(o3_i)
        mol_no_bg = (nox_bg - no2_bg)/mmass(nox_i)
        mol_ox_bg = mol_o3_bg + mol_no2_bg

        mol_o3_emep = o3_emep/mmass(o3_i)
        mol_nox_emep = nox_emep/mmass(nox_i)
        mol_no2_emep = no2_emep/mmass(no2_i)
        mol_no_emep = max(0.0, mol_nox_emep - mol_no2_emep)
        mol_ox_emep = mol_o3_emep + mol_no2_emep
        mol_ox_loc = mol_o3_bg + mol_no2_bg + mol_no2_loc

        if (mol_no2_emep .gt. 0) then
            r = mol_o3_emep*mol_no_emep/mol_no2_emep
        else
            r = 0.0
        end if

        b = mol_ox_loc + mol_nox_bg + mol_nox_loc + r
        c = max(0.0, b**2 - 4.0*mol_ox_loc*(mol_nox_bg + mol_nox_loc)) ! Should never be less than 0 but can be -0
        d = sqrt(c)
        mol_no2_out = (b - d)/2.0
        mol_o3_out = mol_ox_loc - mol_no2_out
        mol_no_out = mol_nox_bg + mol_nox_loc - mol_no2_out

        nox_out = nox_bg + nox_loc
        no2_out = mol_no2_out*mmass(no2_i)
        o3_out = mol_o3_out*mmass(o3_i)

        ! Not correct as it does not calculate the actual photostationary equation
        p_out = r
        p_bg_out = mol_o3_bg*mol_no_bg/mol_no2_bg

        ! Check output
        if (abs(J_photo) > epsilon0 .and. abs(mol_no_out) > epsilon0 .and. abs(mol_o3_out) > epsilon0) then
            p_out = J_photo*mol_no2_out/k1/mol_o3_out/mol_no_out/Na_fac
            p_emep_out = J_photo*mol_no2_emep/k1/mol_o3_emep/mol_no_emep/Na_fac
            p_bg_out = J_photo*mol_no2_bg/k1/mol_o3_bg/mol_no_bg/Na_fac
        else
            p_out = mol_no2_out/(mol_ox_loc + mol_nox_bg + mol_nox_loc - abs(mol_ox_loc - mol_nox_bg - mol_nox_loc))*2.0
            p_emep_out = mol_no2_emep/(mol_ox_emep + mol_nox_emep - abs(mol_ox_emep - mol_nox_emep))*2.0
            p_bg_out = mol_no2_bg/(mol_ox_bg + mol_nox_bg - abs(mol_ox_bg - mol_nox_bg))*2.0
        end if
    end subroutine uEMEP_During_NO2