uEMEP_photostationary_NO2 Subroutine

    subroutine uEMEP_photostationary_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, 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
        real, intent(out) :: nox_out, no2_out, o3_out, p_bg_out, p_out
        integer no2_i, no_i, nox_i, o3_i, ox_i, nox_bg_i, no2_bg_i
        integer, parameter :: n_i = 7
        real :: Na, Na_fac, k1
        real :: mass(n_i)
        real :: mmass(n_i) = [46.0, 30.0, 46.0, 48.0, 47.0, 46.0, 46.0]
        real :: mol(n_i)
        real :: f_no2, f_ox, Jd, fac_sqrt
        real :: min_nox = 1.0e-6


        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))
        
        ! Test the photostationary state for the bg input data
        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)

        ! Test the photostationary state for the input data. Will not be in equilibrium
        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

        ! Set the photolysis rate
        Jd = J_photo/k1/mol(nox_i)

        ! Calculate fraction of NO2 in photostationary state
        fac_sqrt = max(0.0, (1.0 + f_ox + Jd)**2 - 4.0*f_ox)
        f_no2 = 0.5*((1.0 + f_ox + Jd) - sqrt(fac_sqrt))

        ! Convert back to mass
        mol(no2_i) = 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)

        ! 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_photostationary_NO2