uEMEP_Romberg_NO2 Subroutine

    subroutine uEMEP_Romberg_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_out, no2_out, o3_out, romberg_params)
        ! From Norwegian obs fit
        ! real :: a_rom=20
        ! real :: b_rom=30
        ! From model fit
        ! real :: a_rom=30
        ! real :: b_rom=35
        ! real :: c_rom=0.20
        ! Gral values 30 35 0.18
        ! B�chlin and B�singer (2008) 29 35 0.217
        real, intent(in) :: nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc
        real, intent(in) :: romberg_params(3)
        real, intent(out) :: nox_out, no2_out, o3_out

        ! Local variables
        real :: a_rom = 30.0
        real :: b_rom = 35.0
        real :: c_rom = 0.20
        real :: ox_init, no2_init, no2_equ
        
        ! If available, use custom parameter values
        if (abs(romberg_params(1)) > epsilon0) then
            a_rom = romberg_params(1)
            b_rom = romberg_params(2)
            c_rom = romberg_params(3)
        end if

        nox_out = nox_bg + nox_loc
        no2_equ = a_rom*nox_bg/(nox_bg+b_rom) + nox_bg*c_rom
        no2_out = a_rom*nox_out/(nox_out+b_rom) + nox_out*c_rom
        no2_out = no2_out - no2_equ + no2_bg
        no2_out = max(no2_bg, no2_out)
        no2_init = no2_bg + f_no2_loc*nox_loc

        ! Small adjustments for molecular weights
        ox_init = no2_init*47.0/46.0 + o3_bg*47.0/48.0
        o3_out = ox_init*48.0/47.0 - no2_out*48.0/46.0
    end subroutine uEMEP_Romberg_NO2