Skip to content

Commit

Permalink
updates with slow version fj
Browse files Browse the repository at this point in the history
  • Loading branch information
@jialinl6
jialinl6 committed Nov 11, 2024
1 parent 6d0e7da commit fe2c13e
Showing 1 changed file with 67 additions and 114 deletions.
181 changes: 67 additions & 114 deletions src/Fast-JX.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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
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

0 comments on commit fe2c13e

Please sign in to comment.