From fe2c13ebf15bd5e7ca3646b2bd579d97616c8b9f Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Mon, 11 Nov 2024 13:28:00 -0600 Subject: [PATCH] updates with slow version fj --- src/Fast-JX.jl | 181 ++++++++++++++++++------------------------------- 1 file changed, 67 insertions(+), 114 deletions(-) diff --git a/src/Fast-JX.jl b/src/Fast-JX.jl index 6ce1d85c..068b9c29 100644 --- a/src/Fast-JX.jl +++ b/src/Fast-JX.jl @@ -3,50 +3,16 @@ export FastJX # Effective wavelength in 18 bins covering 177–850 nm const WL = SA_F32[187, 191, 193, 196, 202, 208, 211, 214, 261, 267, 277, 295, 303, 310, 316, 333, 380, 574] -function interp2_func(x1, x2, T1, T2) - if x1 ≈ x2 - return (T) -> x1 - end - function interp2(T) - T = clamp(T, T1, T2) - x1 + (T - T1) * (x2 - x1) / (T2 - T1) - end -end - -function interp3_func(x1, x2, x3, T1, T2, T3) - if x1 ≈ x2 && x2 ≈ x3 - return (T) -> x1 - end - function interp3(T) - T = clamp(T, T1, T3) - ifelse(T < T2, - x1 + (T - T1) * (x2 - x1) / (T2 - T1), - x2 + (T - T2) * (x3 - x2) / (T3 - T2), - ) - end -end - -function interp_func(temperatures, cross_sections) - if length(temperatures) == 2 - return interp2_func(cross_sections[1], cross_sections[2], temperatures[1], temperatures[2]) - elseif length(temperatures) == 3 - return interp3_func(cross_sections[1], cross_sections[2], cross_sections[3], - temperatures[1], temperatures[2], temperatures[3]) - end - LinearInterpolation(temperatures, cross_sections, extrapolation_bc=Flat()) -end - """ Create a vector of interpolators to interpolate the cross sections σ (TODO: What are the units?) for different wavelengths (in nm) and temperatures (in K). We use use linear interpolation with flat extrapolation. """ function create_fjx_interp(temperatures::Vector{Float32}, cross_sections::Vector{SVector{18,Float32}}) - [interp_func(temperatures, [x[i] for x ∈ cross_sections]) for i ∈ 1:length(WL)] - #[linear_interpolation(temperatures, [x[i] for x ∈ cross_sections], extrapolation_bc=Flat()) for i ∈ 1:length(WL)] + [linear_interpolation(temperatures, [x[i] for x ∈ cross_sections], extrapolation_bc=Flat()) for i ∈ 1:length(WL)] end -# Cross sections and quantum yield from GEOS-CHEM "FJX_spec.dat" for photo-reactive species included in SuperFast: +# Cross sections and quantum yield from GEOS-CHEM "FJX_spec.dat" for photo-reactive species mentioned in Superfast: # Cross sections σ for O3(1D) for different wavelengths(18 bins), temperatures # O3 -> O2 + O(1D) (superfast: O3 -> 2OH including O(1D) + H2O -> 2OH) @@ -105,27 +71,27 @@ This function is to compute the cosine of the solar zenith angle, given the unix The input variables: lat=latitude(°), long=longitude(°), t=unixtime(s) the cosine of the solar zenith angle (SZA) is given by: . cos(SZA) = sin(LAT)*sin(DEC) + cos(LAT)*cos(DEC)*cos(AHR) - + where LAT = the latitude angle, DEC = the solar declination angle, AHR = the hour angle, all in radians. All in radians """ -function cos_solar_zenith_angle(t, lat, long) +function cos_solar_zenith_angle(lat::T, t::T2, long::T)::T where {T,T2} ut = Dates.unix2datetime(t) DOY = dayofyear(ut) hours = Dates.hour(ut) + Dates.minute(ut) / 60 + Dates.second(ut) / 3600 y = Dates.year(ut) - yearday_temp = 2 * pi * (DOY - 1 + (hours - 12) / 24) - γ = ifelse(mod(y, 4) == 0, # the fraction year in radians - yearday_temp / 366, - yearday_temp / 365, - ) + if mod(y, 4) == 0 + γ = 2 * pi / 366 * (DOY - 1 + (hours - 12) / 24) # the fraction year in radians + else + γ = 2 * pi / 365 * (DOY - 1 + (hours - 12) / 24) # the fraction year in radians + end Eot = 229.18 * (0.000075 + 0.001868 * cos(γ) - 0.032077 * sin(γ) - 0.014615 * cos(γ * 2) - 0.040849 * sin(γ * 2)) timezone = floor(long / 15) - time_offset = Eot + 4 * (long - 15 * timezone) # in minutes - dt = floor(long / 15) # in hours - t_local = t + dt * 3600 # in seconds + time_offset = Eot + 4 * (long - 15 * timezone) # in minute + dt = floor(long / 15) # in hour + t_local = t + dt * 3600 # in second ut_local = Dates.unix2datetime(t_local) tst = Dates.hour(ut_local) + Dates.minute(ut_local) / 60 + Dates.second(ut_local) / 3600 + time_offset / 60 @@ -137,60 +103,54 @@ function cos_solar_zenith_angle(t, lat, long) return CSZA end -@register_symbolic cos_solar_zenith_angle(t, lat, long) - -# Dummy function for unit validation. Basically ModelingToolkit -# will call the function with a DynamicQuantities.Quantity or an integer to -# get information about the type and units of the output. -cos_solar_zenith_angle(t::DynamicQuantities.Quantity, lat, long) = 1.0 - - -# calculate actinic flux at the given cosine of the solar zenith angle `csa` and -# maximium actinic flux `max_actinic_flux` -function calc_flux(csa, max_actinic_flux) +""" +calculate actinic flux at the given cosine of the solar zenith angle `csa` +maximium actinic flux `max_actinic_flux` +""" +function calc_flux(csa::T, max_actinic_flux::T)::T where {T} max(zero(max_actinic_flux), max_actinic_flux * csa) end -# calculate all actinic fluxes -function calc_fluxes(csa) - [calc_flux(csa, actinic_flux[i]) for i in 1:18] -end - -# Symbolic equations for actinic flux -function flux_eqs(csa) - flux_vals = calc_fluxes(csa) - flux_vars = [] - @constants c_flux = 1.0 #[unit = u"s^-1", description = "Constant actinic flux (for unit conversion)"] - for i in 1:18 - wl = WL[i] - n = Symbol("F_", Int(round(wl))) - v = @variables $n(t) #[unit = u"s^-1", description = "Actinic flux at $wl nm"] - push!(flux_vars, only(v)) - end - flux_vars, (flux_vars .~ flux_vals .* c_flux) -end - - -""" +""" Get mean photolysis rates at different times """ -function j_mean(σ_interp, ϕ, Temperature, fluxes) +function j_mean(σ_interp, ϕ::Float32, time::T2, lat::T, long::T, Temperature::T)::T where {T,T2} + csa = cos_solar_zenith_angle(lat, time, long) j = zero(Temperature) - for i in 1:18 - j += fluxes[i] * σ_interp[i](Temperature) * ϕ + @inbounds for i in 1:18 + j += calc_flux(csa, T(actinic_flux[i])) * σ_interp[i](Temperature) * T(ϕ) end j end +function j_mean(σ_interp, ϕ::Float32, time::T2, lat::T, long::T, Temperature::T)::T where {T<:Integer,T2} # Dummy function for symbolic tracing. + round(j_mean(σ_interp, ϕ, Float32(time), Float32(lat), Float32(long), Float32(Temperature))) +end -j_mean_H2O2(T, fluxes) = j_mean(σ_H2O2_interp, ϕ_H2O2_jx, T, fluxes) -j_mean_CH2Oa(T, fluxes) = j_mean(σ_CH2Oa_interp, ϕ_CH2Oa_jx, T, fluxes) -j_mean_CH2Ob(T, fluxes) = j_mean(σ_CH2Ob_interp, ϕ_CH2Ob_jx, T, fluxes) -j_mean_CH3OOH(T, fluxes) = j_mean(σ_CH3OOH_interp, ϕ_CH3OOH_jx, T, fluxes) -j_mean_NO2(T, fluxes) = j_mean(σ_NO2_interp, ϕ_NO2_jx, T, fluxes) -j_mean_o31D(T, fluxes) = j_mean(σ_o31D_interp, ϕ_o31D_jx, T, fluxes) +# Dummy functions for unit validation. Basically ModelingToolkit +# will call the function with a DynamicQuantities.Quantity or an integer to +# get information about the type and units of the output. +j_mean_H2O2(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" +j_mean_CH2Oa(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" +j_mean_CH2Ob(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" +j_mean_CH3OOH(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" +j_mean_NO2(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" +j_mean_o31D(t, lat, long, T::DynamicQuantities.Quantity) = 0.0#u"s^-1" + +j_mean_H2O2(t, lat, long, T) = j_mean(σ_H2O2_interp, ϕ_H2O2_jx, t, lat, long, T) +@register_symbolic j_mean_H2O2(t, lat, long, T) +j_mean_CH2Oa(t, lat, long, T) = j_mean(σ_CH2Oa_interp, ϕ_CH2Oa_jx, t, lat, long, T) +@register_symbolic j_mean_CH2Oa(t, lat, long, T) +j_mean_CH2Ob(t, lat, long, T) = j_mean(σ_CH2Ob_interp, ϕ_CH2Ob_jx, t, lat, long, T) +@register_symbolic j_mean_CH2Ob(t, lat, long, T) +j_mean_CH3OOH(t, lat, long, T) = j_mean(σ_CH3OOH_interp, ϕ_CH3OOH_jx, t, lat, long, T) +@register_symbolic j_mean_CH3OOH(t, lat, long, T) +j_mean_NO2(t, lat, long, T) = j_mean(σ_NO2_interp, ϕ_NO2_jx, t, lat, long, T) +@register_symbolic j_mean_NO2(t, lat, long, T) +j_mean_o31D(t, lat, long, T) = j_mean(σ_o31D_interp, ϕ_o31D_jx, t, lat, long, T) +@register_symbolic j_mean_o31D(t, lat, long, T) """ - adjust_j_o31D(T, P, H2O) + adjust_j_o31D(j_o31D::T, C_H2O::T, C_H2::T, C_N2::T, C_O2::T, T::T) -> T where {T} Adjust the photolysis rate of O3 -> O2 + O(1D) to represent the effective rate for O3 -> 2OH. This adjustment is based on the fraction of O(1D) that reacts with H2O to produce 2 OH. @@ -230,47 +190,40 @@ struct FastJXCoupler end """ -Description: This is a box model used to calculate the photolysis reaction rate constant using the Fast-JX scheme +Description: This is a box model used to calculate the photolysis reaction rate constant using the Fast-JX scheme (Neu, J. L., Prather, M. J., and Penner, J. E. (2007), Global atmospheric chemistry: Integrating over fractional cloud cover, J. Geophys. Res., 112, D11306, doi:10.1029/2006JD008007.) +Build Fast-JX model # Example -Build Fast-JX model: ``` julia fj = FastJX() ``` """ -function FastJX(; name=:FastJX) - @constants T_unit = 1.0 [unit = u"K", description = "Unit temperature (for unit conversion)"] +function FastJX(;name=:FastJX) @parameters T = 298.0 [unit = u"K", description = "Temperature"] + @parameters P = 101325 [unit = u"Pa", description = "Pressure"] @parameters lat = 40.0 [description = "Latitude (Degrees)"] @parameters long = -97.0 [description = "Longitude (Degrees)"] - @parameters P = 101325 [unit = u"Pa", description = "Pressure"] - @parameters H2O = 450 [unit = u"ppb"] + @parameters H2O = 450.0 [unit = u"ppb"] - @variables j_h2o2(t) = 1.0097 * 10.0^-5 #[unit = u"s^-1"] + @variables j_h2o2(t) = 1.0097e-5 #[unit = u"s^-1"] @variables j_CH2Oa(t) = 0.00014 #[unit = u"s^-1"] - @variables j_CH2Ob(t) = 0.00014 #[unit = u"s^-1"] - @variables j_o31D(t) = 4.0 * 10.0^-3 #[unit = u"s^-1"] + @variables j_o31D(t) = 4.0e-3 #[unit = u"s^-1"] @variables j_o32OH(t) = 2.27e-4 #[unit = u"s^-1"] - @variables j_CH3OOH(t) = 8.9573 * 10.0^-6 #[unit = u"s^-1"] + @variables j_CH3OOH(t) = 8.9573e-6 #[unit = u"s^-1"] @variables j_NO2(t) = 0.0149 #[unit = u"s^-1"] - @variables cosSZA(t) [description = "Cosine of the solar zenith angle"] - # TODO(JL): What's difference between the two photolysis reactions of CH2O, do we really need both? - # (@variables j_CH2Ob(t) = 0.00014 [unit = u"s^-1"]) (j_CH2Ob ~ mean_J_CH2Ob(t,lat,T)*j_unit) + @variables j_CH2Ob(t) = 0.00014 #[unit = u"s^-1"] - flux_vars, fluxeqs = flux_eqs(cosSZA) eqs = [ - cosSZA ~ cos_solar_zenith_angle(t, lat, long); - fluxeqs; - j_h2o2 ~ j_mean_H2O2(T/T_unit, flux_vars); - j_CH2Oa ~ j_mean_CH2Oa(T/T_unit, flux_vars); - j_CH2Ob ~ j_mean_CH2Ob(T/T_unit, flux_vars); - j_o31D ~ j_mean_o31D(T/T_unit, flux_vars); - j_o32OH ~ j_o31D*adjust_j_o31D(T, P, H2O); - j_CH3OOH ~ j_mean_CH3OOH(T/T_unit, flux_vars); - j_NO2 ~ j_mean_NO2(T/T_unit, flux_vars) + j_h2o2 ~ j_mean_H2O2(t, lat, long, T) + j_CH2Oa ~ j_mean_CH2Oa(t, lat, long, T) + j_CH2Ob ~ j_mean_CH2Ob(t, lat, long, T) + j_o31D ~ j_mean_o31D(t, lat, long, T)*1e-17 + j_o32OH ~ j_o31D*adjust_j_o31D(T, P, H2O) + j_CH3OOH ~ j_mean_CH3OOH(t, lat, long, T) + j_NO2 ~ j_mean_NO2(t, lat, long, T) ] - ODESystem(eqs, t, [j_h2o2, j_CH2Oa, j_CH2Ob, j_o32OH, j_o31D, j_CH3OOH, j_NO2, cosSZA, flux_vars...], - [lat, long, T, P, H2O]; name=name, metadata=Dict(:coupletype => FastJXCoupler)) -end \ No newline at end of file + ODESystem(eqs, t, [j_h2o2, j_CH2Oa, j_CH2Ob, j_CH3OOH, j_NO2, j_o32OH], [lat, long, T, P, H2O]; name=name, + metadata=Dict(:coupletype => FastJXCoupler)) +end