From a9ec7dc4367e29dc2f672429376ad7579415baa7 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 15:32:12 -0500 Subject: [PATCH 01/32] Dry deposition --- src/DepositionMTK.jl | 12 +-- src/dry_deposition.jl | 233 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 237 insertions(+), 8 deletions(-) create mode 100644 src/dry_deposition.jl diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index 80ed2cfa..b1d4ba66 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -1,6 +1,3 @@ -module DepositionMTK -import Base.show, StaticArrays -export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, inf, r_i, r_lu, r_ac, r_gsS, r_gsO, r_clS, r_clO, So2Data, O3Data, No2Data, NoData, Hno3Data, H2o2Data, GasData, AldData, HchoData, OpData, PaaData, OraData, Nh3Data, PanData, Hno2Data using StaticArrays @@ -89,7 +86,7 @@ export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. - function r_s(G::T, Ts::T, iSeason::Int, iLandUse::Int, rainOrDew::Bool) where T<:AbstractFloat + function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) rs = 0.0 if Ts >= 39.9 || Ts <= 0.1 rs = inf @@ -105,7 +102,7 @@ export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, end # Calculate the resistance from the effects of mixing forced by buoyant convection when sunlight heats the ground or lower canopy and by penetration of wind into canopies on the sides of hills [s m-1] when given the solar irradiation (G [W m-2]) and the slope of the local terrain (θ [radians]). From Wesely (1989) equation 5. - function r_dc(G::T, θ::T) where T<:AbstractFloat + function r_dc(G, θ) return 100.0 * (1.0 + 1000.0/(G+10.0)) / (1.0 + 1000.0*θ) end @@ -172,8 +169,8 @@ export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, # Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3). # From Wesely (1989) regarding rain and dew inputs: "A direct computation of the surface wetness would be most desirable, e.g. by estimating the amount of free surface water accumulated and then evaporated. Alternatively, surface relative humidity might be a useful index. After dewfall and rainfall events are completed, surface wetness often disappears as a result of evaporation after approximately 2 hours of good atmospheric mixing, the period of time recommended earlier (Sheih et al., 1986)". - function SurfaceResistance(gasData::GasData, G::T, Ts::T, θ::T, - iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where T<:AbstractFloat + function SurfaceResistance(gasData::GasData, G, Ts, θ, + iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) rs = r_s(G, Ts, iSeason, iLandUse, rain || dew) rmx = r_mx(gasData.Hstar, gasData.Fo) rsmx = r_smx(rs, gasData.Dh2oPerDx, rmx) @@ -208,4 +205,3 @@ export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, r_c = min(r_c, 9999.0) return r_c end -end diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl new file mode 100644 index 00000000..3437fb95 --- /dev/null +++ b/src/dry_deposition.jl @@ -0,0 +1,233 @@ +using Unitful +using ModelingToolkit +using StaticArrays +using Test +include("DepositionMTK.jl") + +g = 9.81u"m*s^-2" # gravitational acceleration [m/s2] +κ = 0.4 # von Karman constant +k = 1.3806488e-23u"m^2*kg*s^-2/K" # Boltzmann constant +M_air = 28.97e-3u"kg/mol" # molecular weight of air +R = 8.3144621u"J/K/mol" # Gas constant + +""" +Function Ra calculates aerodynamic resistance to dry deposition +where z is the top of the surface layer [m], z₀ is the roughness length [m], u_star is friction velocity [m/s], and L is Monin-Obukhov length [m] +Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 & 19.14. +""" + +function ra(z, z₀, u_star, L) + ζ = z/L + ζ₀= z₀/L + if 0 < ζ < 1 + rₐ = 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)) + elseif ζ == 0 + rₐ = 1/(κ*u_star)*log(z/z₀) + elseif -1 < ζ < 0 + η₀=(1-15*ζ₀)^1/4 + η =(1-15*ζ)^1/4 + rₐ = 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))] + else + print("wrong rₐ") + end + return rₐ +end + +""" +Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. +""" +function mu(T) + return (1.8e-5*(T/298u"K")^0.85)u"kg/m/s" +end + +""" +Function mfp calculates the mean free path of air [m] +where T is temperature [K] P is pressure [Pa], and Mu is dynamic viscosity [kg/(m s)]. +From Seinfeld and Pandis (2006) equation 9.6 +""" +function mfp(T,P,μ) + return 2*μ/(P*(8*M_air/(pi*R*T))^0.5) +end + +""" +Function cc calculates the Cunnningham slip correction factor +where Dp is particle diameter [m], T is temperature [K], and P is pressure [Pa]. +From Seinfeld and Pandis (2006) equation 9.34. +""" +function cc(Dₚ,T,P,μ) + λ = mfp(T,P,μ) + return 1+2*λ/Dₚ*(1.257+0.4*exp(-1.1*Dₚ/(2*λ))) +end + +""" +Function vs calculates the terminal setting velocity of a +particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. +From equation 9.42 in Seinfeld and Pandis (2006). +""" +function vs(Dₚ, ρₚ, Cc, μ) + if Dₚ > 20.e-6u"m" + print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.") + else + return Dₚ^2*ρₚ*g*Cc/(18*μ) + end +end + +""" +Function dParticle calculates the brownian diffusivity of a particle [m2/s] using the Stokes-Einstein-Sutherland relation +(Seinfeld and Pandis eq. 9.73) +where T is air temperature [K], P is pressure [Pa], Dp is particle diameter [m], and μ is air dynamic viscosity [kg/(s m)] +""" +function dParticle(T,P,Dₚ,Cc,μ) + return k*T*Cc/(3*pi*μ*Dₚ) +end + +""" +Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K] +using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf +""" +function dH2O(T) + T_unitless = T*1u"K^-1" + return (-2.775e-6 + 4.479e-8*T_unitless + 1.656e-10*T_unitless^2)u"m^2/s" +end + +""" +Function sc computes the dimensionless Schmidt number, +where μ is dynamic viscosity of air [kg/(s m)], ρ is air density [kg/m3], and D is the molecular diffusivity of the gas speciesof interest [m2/s] +""" +function sc(μ, ρ, D) + return μ/(ρ*D) +end + +""" +Function stSmooth computes the dimensionless Stokes number for dry deposition of particles on smooth surfaces or surfaces with bluff roughness elements, +where vs is settling velocity [m/s], u_star is friction velocity [m/s], μ is dynamic viscosity of air [kg/(s m)], and ρ is air density [kg/m3], +based on Seinfeld and Pandis (2006) equation 19.23. +""" +function stSmooth(vₛ, u_star, μ, ρ) + return vₛ*u_star^2/(g*ρ*μ) +end + +""" +Function stVeg computes the dimensionless Stokes number for dry deposition of particles on vegetated surfaces, +where vs is settling velocity [m/s], u_star is friction velocity [m/s], and A is the characteristic collector radius [m], +based on Seinfeld and Pandis (2006) equation 19.24. +""" +function stVeg(vₛ, u_star, A) + return vₛ*u_star/(g*A) +end + +""" +Function RbGas calculates the quasi-laminar sublayer resistance to dry deposition for a gas species [s/m], +where Sc is the dimensionless Schmidt number and u_star is the friction velocity [m/s]. +From Seinfeld and Pandis (2006) equation 19.17. +""" +function RbGas(Sc, u_star) + return 5*Sc^(2/3)/u_star +end + +""" +Values for the characteristic radii of collectors [m] +where the columns are land use categories and the rows are seasonal categories. +Land-use categories (LUCs) +1. Evergreen–needleleaf trees +2. Deciduous broadleaf trees +3. Grass +4. Desert +5. Shrubs and interrupted woodlands +Seasonal categories (SC) +1. Midsummer with lush vegetation +2. Autumn with cropland not harvested +3. Late autumn after frost, no snow +4. Winter, snow on ground +5. Transitional +given in Seinfeld and Pandis Table 19.2 +""" +z₀_table = SA_F32[ + 0.8 1.05 0.1 0.04 0.1 + 0.9 1.05 0.1 0.04 0.1 + 0.9 0.95 0.05 0.04 0.1 + 0.9 0.55 0.02 0.04 0.1 + 0.8 0.75 0.05 0.04 0.1 +] # unit:[m] + +A_table = SA_F32[ + 2.0 5.0 2.0 Inf 10.0 + 2.0 5.0 2.0 Inf 10.0 + 2.0 10.0 5.0 Inf 10.0 + 2.0 10.0 2.0 Inf 10.0 + 2.0 5.0 2.0 Inf 10.0 +] # unit:[mm] + +α_table = SA_F32[ + 1.0 0.8 1.2 50.0 1.3 +] + +γ_table = SA_F32[ + 0.56 0.56 0.54 0.54 0.54 +] + +""" +Function RbParticle calculates the quasi-laminar sublayer resistance to dry deposition for a particles [s/m], +where Sc is the dimensionless Schmidt number, u_star is the friction velocity [m/s], St is the dimensionless Stokes number, +Dp is particle diameter [m], and iSeason and iLandUse are season and land use indexes, respectively. +From Seinfeld and Pandis (2006) equation 19.27. +""" + +function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) + α = α_table[iLandUse] + γ = γ_table[iLandUse] + A = (A_table[iSeason,iLandUse]*10^(-3))u"m" + R1 = exp(-St^0.5) + term_1 = Sc^(-γ) + term_2 = (St/(α+St))^2 + term_3 = 1/2*(Dₚ/A)^2 + return 1/(3*u_star*(term_1+term_2+term_3)*R1) +end + +""" +Function DryDepGas calculates dry deposition velocity [m/s] for a gas species, +where z is the height of the surface layer [m], zo is roughness length [m], u_star is friction velocity [m/s], +L is Monin-Obukhov length [m], T is surface air temperature [K], ρA is air density [kg/m3] +gasData is data about the gas species for surface resistance calculations, G is solar +irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and iLandUse are indexes for the season and land use, +dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. +Based on Seinfeld and Pandis (2006) equation 19.2. +""" +function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) + Ra = ra(z, z₀, u_star, L) + μ = mu(Ts) + Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s] + Sc = sc(μ,ρA, Dg) + Rb = RbGas(Sc, u_star) + Rc = SurfaceResistance(gasData, G/u"W*m^-2", Ts/u"K", θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m" + return 1/(Ra+Rb+Rc) +end + +""" +Function DryDepParticle calculates particle dry deposition velocity [m/s] +where z is the height of the surface layer [m], zo is roughness length [m], u_star is friction velocity [m/s], L is Monin-Obukhov length [m], +Dp is particle diameter [m], Ts is surface air temperature [K], P is pressure [Pa], ρParticle is particle density [kg/m3], ρAir is air density [kg/m3], +and iSeason and iLandUse are indexes for the season and land use. +Based on Seinfeld and Pandis (2006) equation 19.7. +""" +function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason::Int, iLandUse::Int) + Ra = ra(z, z₀, u_star, L) + μ = mu(Ts) + Cc = cc(Dp, Ts, P, μ) + Vs = vs(Dp, ρParticle, Cc, μ) + if iLandUse == 4 # dessert + St = stSmooth(Vs, u_star, μ, ρA) + else + St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))u"m") + end + D = dParticle(Ts, P, Dp, Cc, μ) + Sc = sc(μ, ρA, D) + Rb = RbParticle(Sc, u_star, St, Dp, iSeason, iLandUse) + return 1/(Ra+Rb+Ra*Rb*Vs)+Vs +end + +@test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" +@test unit(dH2O(300u"K"))==u"m^2/s" +@test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" +@test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" + From d5dcb7fb1dad2b99b1f44a7238b2517e66b0578e Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 18:41:16 -0500 Subject: [PATCH 02/32] Finish Unit Tests --- src/DepositionMTK.jl | 211 ++-------------------------------------- src/dry_deposition.jl | 23 ++--- src/wesley1989.jl | 207 +++++++++++++++++++++++++++++++++++++++ test/drydep_test.jl | 79 +++++++++++++++ test/runtests.jl | 165 ++----------------------------- test/wesely1989_test.jl | 163 +++++++++++++++++++++++++++++++ 6 files changed, 476 insertions(+), 372 deletions(-) create mode 100644 src/wesley1989.jl create mode 100644 test/drydep_test.jl create mode 100644 test/wesely1989_test.jl diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index b1d4ba66..a45e0f33 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -1,207 +1,10 @@ +module DepositionMTK - using StaticArrays +using Unitful +using ModelingToolkit +using StaticArrays - const inf = 1.e25 +include("wesley1989.jl") +include("dry_deposition.jl") - # r_i represents the minimum bulk canopy resistances for water vapor. - const r_i = SA_F32[ - inf 60 120 70 130 100 inf inf 80 100 150 - inf inf inf inf 250 500 inf inf inf inf inf - inf inf inf inf 250 500 inf inf inf inf inf - inf inf inf inf 400 800 inf inf inf inf inf - inf 120 240 140 250 190 inf inf 160 200 300] - - # r_lu signifies leaf cuticles in healthy vegetation and otherwise the outer surfaces in the upper canopy. - const r_lu = SA_F32[ - inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000 - inf 9000 9000 9000 4000 8000 inf inf 9000 9000 9000 - inf inf 9000 9000 4000 8000 inf inf 9000 9000 9000 - inf inf inf inf 6000 9000 inf inf 9000 9000 9000 - inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000] - - # r_ac signifies transfer that depends only on canopy height and density. - const r_ac = SA_F32[ - 100.0 200 100 2000 2000 2000 0 0 300 150 200 - 100 150 100 1500 2000 1700 0 0 200 120 140 - 100 10 100 1000 2000 1500 0 0 100 50 120 - 100 10 10 1000 2000 1500 0 0 50 10 50 - 100 50 80 1200 2000 1500 0 0 200 60 120] - - # r_gs signifies uptake at the "ground" by soil, leaf litter, snow, water etc. 'S' and 'O' stand for SO2 and O3 respectively. - const r_gsS = SA_F32[ - 400.0 150 350 500 500 100 0 1000 0 220 40 - 400 200 350 500 500 100 0 1000 0 300 400 - 400 150 350 500 500 200 0 1000 0 200 400 - 100 100 100 100 100 100 0 1000 100 100 50 - 500 150 350 500 500 200 0 1000 0 250 40] - - const r_gsO = SA_F32[ - 300.0 150 200 200 200 300 2000 400 1000 180 200 - 300 150 200 200 200 300 2000 400 800 180 200 - 300 150 200 200 200 300 2000 400 1000 180 20 - 600 3500 3500 3500 3500 3500 2000 400 3500 3500 3500 - 300 150 200 200 200 300 2000 400 1000 180 200] - - # r_cl is meant to account for uptake pathways at the leaves, bark, etc. 'S' and 'O' stand for SO2 and O3 respectively. - const r_clS = SA_F32[ - inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000 - inf 9000 9000 9000 2000 4000 inf inf 9000 9000 9000 - inf inf 9000 9000 3000 6000 inf inf 9000 9000 9000 - inf inf inf 9000 200 400 inf inf 9000 inf 9000 - inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000] - - const r_clO = SA_F32[ - inf 1000 1000 1000 1000 1000 inf inf 1000 1000 1000 - inf 400 400 400 1000 600 inf inf 400 400 400 - inf 1000 400 400 1000 600 inf inf 800 600 600 - inf 1000 1000 400 1500 600 inf inf 800 1000 800 - inf 1000 500 500 1500 700 inf inf 600 800 800] - - # Holder for gas properties from Wesely (1989) Table 2.' - struct GasData - Dh2oPerDx::AbstractFloat - Hstar::AbstractFloat - Fo::AbstractFloat - end - - const So2Data = GasData(1.9, 1.e5, 0) - const O3Data = GasData(1.6, 0.01, 1) - const No2Data = GasData(1.6, 0.01, 0.1) # Wesely (1989) suggests that, - # in general, the sum of NO and NO2 should be considered rather - # than NO2 alone because rapid in-air chemical reactions can cause - # a significant change of NO and NO2 vertical fluxes between the - # surface and the point at which deposition velocities are applied, - # but the sum of NO and NO2 fluxes should be practically unchanged. - const NoData = GasData(1.3, 3.e-3, 0) # Changed according to Walmsley (1996) - const Hno3Data = GasData(1.9, 1.e14, 0) - const H2o2Data = GasData(1.4, 1.e5, 1) - const AldData = GasData(1.6, 15, 0) # Acetaldehyde (aldehyde class) - const HchoData = GasData(1.3, 6.e3, 0) # Formaldehyde - const OpData = GasData(1.6, 240, 0.1) - # Methyl hydroperoxide (organic peroxide class) - const PaaData = GasData(2.0, 540, 0.1) # Peroxyacetyl nitrate - const OraData = GasData(1.6, 4.e6, 0) # Formic acid (organic acid class) - const Nh3Data = GasData(0.97, 2.e4, 0) # Changed according to Walmsley (1996) - const PanData = GasData(2.6, 3.6, 0.1) # Peroxyacetyl nitrate - const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid - - # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. - function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) - rs = 0.0 - if Ts >= 39.9 || Ts <= 0.1 - rs = inf - else - rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * - (400.0 * 1.0 / (Ts * (40.0 - Ts))) - end - # Adjust for dew and rain (from "Effects of dew and rain" section). - if rainOrDew - rs *= 3 - end - return rs - end - - # Calculate the resistance from the effects of mixing forced by buoyant convection when sunlight heats the ground or lower canopy and by penetration of wind into canopies on the sides of hills [s m-1] when given the solar irradiation (G [W m-2]) and the slope of the local terrain (θ [radians]). From Wesely (1989) equation 5. - function r_dc(G, θ) - return 100.0 * (1.0 + 1000.0/(G+10.0)) / (1.0 + 1000.0*θ) - end - - #Calculate mesophyll resistance [s m-1] based on Wesely (1989) equation 6 when given the effective Henry's law coefficient (Hstar [M atm-1]) and the reactivity factor (fo [-]), both available in Wesely (1989) table 2. - function r_mx(Hstar::T, fo::T) where T<:AbstractFloat - return 1.0 / (Hstar/3000.0 + 100.0*fo) - end - - #Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]). - function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat - return r_s * Dh2oPerDx + r_mx - end - - # Calculate the resistance of the outer surfaces in the upper canopy (leaf cuticular resistance in healthy vegetation) based on Wesely (1989) equations 7 and 10-14 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 or O3. - function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where T<:AbstractFloat - rlux = 0.0 - if dew && (iSeason != 4) # Dew doesn't have any effect in the winter - if isSO2 - if iLandUse == 1 - rlux = 50.0 # equation 13 and a half - else - rlux = 100 # equation 10. - end - elseif isO3 - # equation 11 - rlux = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) - else - rluO = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) #equation 11 - rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.0e-5*Hstar+fo)) - + 1.0e-7*Hstar + fo/rluO) # equation 14, modified to match Walmsley eq. 5g - end - elseif rain && (iSeason != 4) - if isSO2 - if iLandUse == 1 - rlux = 50 #equation 13 and a half - else - # equation 12 - rlux = 1.0 / (1.0/5000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) - end - elseif isO3 - # equation 13 - rlux = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) - else - rluO = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) # equation 13 - rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.e-5*Hstar+fo)) + 1.0e-7*Hstar +fo/rluO) # equation 14, modified to match Walmsley eq. 5g - end - else - rlux = r_lu[iSeason, iLandUse] / (1.0e-5*Hstar + fo) - end - return rlux - end - - # Calculate the resistance of the exposed surfaces in the lower portions of structures (canopies, buildings) above the ground based on Wesely (1989) equation 8 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). - function r_clx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat - return 1.0 / (Hstar/(1.0e5*r_clS[iSeason, iLandUse]) + - fo/r_clO[iSeason, iLandUse]) - end - - # Calculate the resistance to uptake at the 'ground' surface based on Wesely (1989) equation 9 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). - function r_gsx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat - return 1.0 / (Hstar/(1.0e5*r_gsS[iSeason, iLandUse]) + - fo/r_gsO[iSeason, iLandUse]) - end - - # Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3). - # From Wesely (1989) regarding rain and dew inputs: "A direct computation of the surface wetness would be most desirable, e.g. by estimating the amount of free surface water accumulated and then evaporated. Alternatively, surface relative humidity might be a useful index. After dewfall and rainfall events are completed, surface wetness often disappears as a result of evaporation after approximately 2 hours of good atmospheric mixing, the period of time recommended earlier (Sheih et al., 1986)". - function SurfaceResistance(gasData::GasData, G, Ts, θ, - iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) - rs = r_s(G, Ts, iSeason, iLandUse, rain || dew) - rmx = r_mx(gasData.Hstar, gasData.Fo) - rsmx = r_smx(rs, gasData.Dh2oPerDx, rmx) - rdc = r_dc(G, θ) - rlux = r_lux(gasData.Hstar, gasData.Fo, iSeason, iLandUse, - rain, dew, isSO2, isO3) - rclx = 0.0 - rgsx = 0.0 - if isSO2 - rclx = r_clS[iSeason, iLandUse] - rgsx = r_gsS[iSeason, iLandUse] - elseif isO3 - rclx = r_clO[iSeason, iLandUse] - rgsx = r_gsO[iSeason, iLandUse] - else - rclx = r_clx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) - rgsx = r_gsx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) - end - - rac = r_ac[iSeason, iLandUse] - - # Correction for cold temperatures from page 4 column 1. - if Ts < 0.0 - correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark - rlux += correction - rclx += correction - rgsx += correction - end - - r_c = 1.0 / (1.0/(rsmx) + 1.0/rlux + 1.0/(rdc+rclx) + 1.0/(rac+rgsx)) - r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces. - r_c = min(r_c, 9999.0) - return r_c - end +end \ No newline at end of file diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 3437fb95..5b76abb5 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,8 +1,8 @@ using Unitful using ModelingToolkit using StaticArrays -using Test -include("DepositionMTK.jl") +include("wesley1989.jl") +export ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle g = 9.81u"m*s^-2" # gravitational acceleration [m/s2] κ = 0.4 # von Karman constant @@ -17,8 +17,14 @@ Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 """ function ra(z, z₀, u_star, L) - ζ = z/L - ζ₀= z₀/L + if L == 0u"m" + ζ = 0 + ζ₀= 0 + else + ζ = z/L + ζ₀= z₀/L + end + print(ζ) if 0 < ζ < 1 rₐ = 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)) elseif ζ == 0 @@ -104,7 +110,7 @@ where vs is settling velocity [m/s], u_star is friction velocity [m/s], μ is dy based on Seinfeld and Pandis (2006) equation 19.23. """ function stSmooth(vₛ, u_star, μ, ρ) - return vₛ*u_star^2/(g*ρ*μ) + return vₛ*u_star^2*ρ/(g*μ) end """ @@ -199,7 +205,7 @@ function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s] Sc = sc(μ,ρA, Dg) Rb = RbGas(Sc, u_star) - Rc = SurfaceResistance(gasData, G/u"W*m^-2", Ts/u"K", θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m" + Rc = SurfaceResistance(gasData, G/u"W*m^-2", (Ts/u"K"-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m" return 1/(Ra+Rb+Rc) end @@ -226,8 +232,3 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: return 1/(Ra+Rb+Ra*Rb*Vs)+Vs end -@test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" -@test unit(dH2O(300u"K"))==u"m^2/s" -@test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" -@test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" - diff --git a/src/wesley1989.jl b/src/wesley1989.jl new file mode 100644 index 00000000..8835ea43 --- /dev/null +++ b/src/wesley1989.jl @@ -0,0 +1,207 @@ + export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, inf, r_i, r_lu, r_ac, r_gsS, r_gsO, r_clS, r_clO, So2Data, O3Data, No2Data, NoData, Hno3Data, H2o2Data, GasData, AldData, HchoData, OpData, PaaData, OraData, Nh3Data, PanData, Hno2Data + + + const inf = 1.e25 + + # r_i represents the minimum bulk canopy resistances for water vapor. + const r_i = SA_F32[ + inf 60 120 70 130 100 inf inf 80 100 150 + inf inf inf inf 250 500 inf inf inf inf inf + inf inf inf inf 250 500 inf inf inf inf inf + inf inf inf inf 400 800 inf inf inf inf inf + inf 120 240 140 250 190 inf inf 160 200 300] + + # r_lu signifies leaf cuticles in healthy vegetation and otherwise the outer surfaces in the upper canopy. + const r_lu = SA_F32[ + inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000 + inf 9000 9000 9000 4000 8000 inf inf 9000 9000 9000 + inf inf 9000 9000 4000 8000 inf inf 9000 9000 9000 + inf inf inf inf 6000 9000 inf inf 9000 9000 9000 + inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000] + + # r_ac signifies transfer that depends only on canopy height and density. + const r_ac = SA_F32[ + 100.0 200 100 2000 2000 2000 0 0 300 150 200 + 100 150 100 1500 2000 1700 0 0 200 120 140 + 100 10 100 1000 2000 1500 0 0 100 50 120 + 100 10 10 1000 2000 1500 0 0 50 10 50 + 100 50 80 1200 2000 1500 0 0 200 60 120] + + # r_gs signifies uptake at the "ground" by soil, leaf litter, snow, water etc. 'S' and 'O' stand for SO2 and O3 respectively. + const r_gsS = SA_F32[ + 400.0 150 350 500 500 100 0 1000 0 220 40 + 400 200 350 500 500 100 0 1000 0 300 400 + 400 150 350 500 500 200 0 1000 0 200 400 + 100 100 100 100 100 100 0 1000 100 100 50 + 500 150 350 500 500 200 0 1000 0 250 40] + + const r_gsO = SA_F32[ + 300.0 150 200 200 200 300 2000 400 1000 180 200 + 300 150 200 200 200 300 2000 400 800 180 200 + 300 150 200 200 200 300 2000 400 1000 180 20 + 600 3500 3500 3500 3500 3500 2000 400 3500 3500 3500 + 300 150 200 200 200 300 2000 400 1000 180 200] + + # r_cl is meant to account for uptake pathways at the leaves, bark, etc. 'S' and 'O' stand for SO2 and O3 respectively. + const r_clS = SA_F32[ + inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000 + inf 9000 9000 9000 2000 4000 inf inf 9000 9000 9000 + inf inf 9000 9000 3000 6000 inf inf 9000 9000 9000 + inf inf inf 9000 200 400 inf inf 9000 inf 9000 + inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000] + + const r_clO = SA_F32[ + inf 1000 1000 1000 1000 1000 inf inf 1000 1000 1000 + inf 400 400 400 1000 600 inf inf 400 400 400 + inf 1000 400 400 1000 600 inf inf 800 600 600 + inf 1000 1000 400 1500 600 inf inf 800 1000 800 + inf 1000 500 500 1500 700 inf inf 600 800 800] + + # Holder for gas properties from Wesely (1989) Table 2.' + struct GasData + Dh2oPerDx::AbstractFloat + Hstar::AbstractFloat + Fo::AbstractFloat + end + + const So2Data = GasData(1.9, 1.e5, 0) + const O3Data = GasData(1.6, 0.01, 1) + const No2Data = GasData(1.6, 0.01, 0.1) # Wesely (1989) suggests that, + # in general, the sum of NO and NO2 should be considered rather + # than NO2 alone because rapid in-air chemical reactions can cause + # a significant change of NO and NO2 vertical fluxes between the + # surface and the point at which deposition velocities are applied, + # but the sum of NO and NO2 fluxes should be practically unchanged. + const NoData = GasData(1.3, 3.e-3, 0) # Changed according to Walmsley (1996) + const Hno3Data = GasData(1.9, 1.e14, 0) + const H2o2Data = GasData(1.4, 1.e5, 1) + const AldData = GasData(1.6, 15, 0) # Acetaldehyde (aldehyde class) + const HchoData = GasData(1.3, 6.e3, 0) # Formaldehyde + const OpData = GasData(1.6, 240, 0.1) + # Methyl hydroperoxide (organic peroxide class) + const PaaData = GasData(2.0, 540, 0.1) # Peroxyacetyl nitrate + const OraData = GasData(1.6, 4.e6, 0) # Formic acid (organic acid class) + const Nh3Data = GasData(0.97, 2.e4, 0) # Changed according to Walmsley (1996) + const PanData = GasData(2.6, 3.6, 0.1) # Peroxyacetyl nitrate + const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid + + # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. + function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) + rs = 0.0 + if Ts >= 39.9 || Ts <= 0.1 + rs = inf + else + rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * + (400.0 * 1.0 / (Ts * (40.0 - Ts))) + end + # Adjust for dew and rain (from "Effects of dew and rain" section). + if rainOrDew + rs *= 3 + end + return rs + end + + # Calculate the resistance from the effects of mixing forced by buoyant convection when sunlight heats the ground or lower canopy and by penetration of wind into canopies on the sides of hills [s m-1] when given the solar irradiation (G [W m-2]) and the slope of the local terrain (θ [radians]). From Wesely (1989) equation 5. + function r_dc(G, θ) + return 100.0 * (1.0 + 1000.0/(G+10.0)) / (1.0 + 1000.0*θ) + end + + #Calculate mesophyll resistance [s m-1] based on Wesely (1989) equation 6 when given the effective Henry's law coefficient (Hstar [M atm-1]) and the reactivity factor (fo [-]), both available in Wesely (1989) table 2. + function r_mx(Hstar::T, fo::T) where T<:AbstractFloat + return 1.0 / (Hstar/3000.0 + 100.0*fo) + end + + #Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]). + function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat + return r_s * Dh2oPerDx + r_mx + end + + # Calculate the resistance of the outer surfaces in the upper canopy (leaf cuticular resistance in healthy vegetation) based on Wesely (1989) equations 7 and 10-14 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 or O3. + function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where T<:AbstractFloat + rlux = 0.0 + if dew && (iSeason != 4) # Dew doesn't have any effect in the winter + if isSO2 + if iLandUse == 1 + rlux = 50.0 # equation 13 and a half + else + rlux = 100 # equation 10. + end + elseif isO3 + # equation 11 + rlux = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) + else + rluO = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) #equation 11 + rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.0e-5*Hstar+fo)) + + 1.0e-7*Hstar + fo/rluO) # equation 14, modified to match Walmsley eq. 5g + end + elseif rain && (iSeason != 4) + if isSO2 + if iLandUse == 1 + rlux = 50 #equation 13 and a half + else + # equation 12 + rlux = 1.0 / (1.0/5000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) + end + elseif isO3 + # equation 13 + rlux = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) + else + rluO = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) # equation 13 + rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.e-5*Hstar+fo)) + 1.0e-7*Hstar +fo/rluO) # equation 14, modified to match Walmsley eq. 5g + end + else + rlux = r_lu[iSeason, iLandUse] / (1.0e-5*Hstar + fo) + end + return rlux + end + + # Calculate the resistance of the exposed surfaces in the lower portions of structures (canopies, buildings) above the ground based on Wesely (1989) equation 8 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). + function r_clx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat + return 1.0 / (Hstar/(1.0e5*r_clS[iSeason, iLandUse]) + + fo/r_clO[iSeason, iLandUse]) + end + + # Calculate the resistance to uptake at the 'ground' surface based on Wesely (1989) equation 9 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). + function r_gsx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat + return 1.0 / (Hstar/(1.0e5*r_gsS[iSeason, iLandUse]) + + fo/r_gsO[iSeason, iLandUse]) + end + + # Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3). + # From Wesely (1989) regarding rain and dew inputs: "A direct computation of the surface wetness would be most desirable, e.g. by estimating the amount of free surface water accumulated and then evaporated. Alternatively, surface relative humidity might be a useful index. After dewfall and rainfall events are completed, surface wetness often disappears as a result of evaporation after approximately 2 hours of good atmospheric mixing, the period of time recommended earlier (Sheih et al., 1986)". + function SurfaceResistance(gasData::GasData, G, Ts, θ, + iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) + rs = r_s(G, Ts, iSeason, iLandUse, rain || dew) + rmx = r_mx(gasData.Hstar, gasData.Fo) + rsmx = r_smx(rs, gasData.Dh2oPerDx, rmx) + rdc = r_dc(G, θ) + rlux = r_lux(gasData.Hstar, gasData.Fo, iSeason, iLandUse, + rain, dew, isSO2, isO3) + rclx = 0.0 + rgsx = 0.0 + if isSO2 + rclx = r_clS[iSeason, iLandUse] + rgsx = r_gsS[iSeason, iLandUse] + elseif isO3 + rclx = r_clO[iSeason, iLandUse] + rgsx = r_gsO[iSeason, iLandUse] + else + rclx = r_clx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) + rgsx = r_gsx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) + end + + rac = r_ac[iSeason, iLandUse] + + # Correction for cold temperatures from page 4 column 1. + if Ts < 0.0 + correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark + rlux += correction + rclx += correction + rgsx += correction + end + + r_c = 1.0 / (1.0/(rsmx) + 1.0/rlux + 1.0/(rdc+rclx) + 1.0/(rac+rgsx)) + r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces. + r_c = min(r_c, 9999.0) + return r_c + end diff --git a/test/drydep_test.jl b/test/drydep_test.jl new file mode 100644 index 00000000..10b47d35 --- /dev/null +++ b/test/drydep_test.jl @@ -0,0 +1,79 @@ +using DepositionMTK +using Test + +@testset "mfp" begin + @test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" +end + +@testset "unit" begin + @test unit(dH2O(300u"K"))==u"m^2/s" + @test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" + @test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" +end + +@testset "viscosity" begin + T = [100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 1200, 1300, 1400, 1500, 1600] + μ_list = [0.6924, 1.0283, 1.3289, 1.488, 1.983, 2.075, 2.286, 2.484, 2.671, 2.848, 3.018, 3.177, 3.332, 3.481, 3.625, 3.765, 3.899, 4.023, 4.152, 4.44, 4.69, 4.93, 5.17, 5.4, 5.63] + μ_test = [] + for i in 1:25 + push!(μ_test, (mu((T[i])u"K"))u"m*s/kg") + end + @test μ_list ≈ μ_list rtol=1e-2 +end + +@testset "Cc" begin + Dp = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] + Cc_list = [216, 108, 43.6, 22.2, 11.4, 4.95, 2.85, 1.865, 1.326, 1.164, 1.082, 1.032, 1.016, 1.008, 1.003, 1.0016] + Cc_test = [] + for i in 1:16 + push!(Cc_test, cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K"))) + end + @test Cc_test ≈ Cc_list rtol=1e-2 +end + +@testset "Vs" begin + Dp = [0.01, 0.1, 1, 10] + Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 + Vs_test = [] + for i in 1:4 + push!(Vs_test, (vs((Dp[i]*1e-6)u"m", 1000u"kg*m^-3", cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K")), mu(298u"K")))u"s/m") + end + @test Vs_test ≈ Vs_list rtol=0.1 +end + +@testset "DryDepGas" begin + z = 50u"m" + z₀ = 0.04u"m" + u_star = 0.44u"m/s" + L = 0u"m" + T = 298u"K" + ρA = 1.2u"kg*m^-3" + G = 300u"W*m^-2" + θ = 0 + vd_true = 0.003u"m/s" + @test DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, 1, 10, false, false, false, false) ≈ vd_true rtol=0.33 +end + +@testset "DryDepParticle" begin + z = 20u"m" + z₀ = 0.02u"m" + u_star = 0.44u"m/s" + L = 0u"m" + T = 298u"K" + P = 101325u"Pa" + ρA = 1.2u"kg*m^-3" + ρParticle = 1000u"kg*m^-3" + Dp = [1.e-8, 1.e-7, 1.e-6, 1.e-5] + vd_true = [0.5, 0.012, 0.02, 0.6] # [cm/s] + vd_list = [] + for i in 1:4 + push!(vd_list, (DryDepParticle(z, z₀, u_star, L, (Dp[i])u"m", T, P, ρParticle, ρA, 1, 4))*100u"s/m") + end + @test vd_list ≈ vd_true rtol=0.8 # 80% difference, not ideal +end + + + + + + \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 69b5bbf4..3f0d4375 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,161 +1,12 @@ using DepositionMTK -using Test -using StaticArrays - -# include("../src/DepositionMTK.jl") - -# Results from Wesely (1989) table 3; updated to values from Walmsley (1996) table 1. - -const SO2 = SA_F32[ - 130.0 140 160 380 1000 100 1200 - 1400 1400 1400 1400 1500 100 1300 - 1100 1100 1100 1100 1200 90 1000 - 1000 1000 1000 1000 1100 1100 1100 - 270 290 330 620 1100 90 1000] - -const O3 = SA_F32[ - 100.0 110 130 320 960 960 580 - 430 470 520 710 1300 950 580 - 390 420 460 610 960 770 510 - 560 620 710 1100 3200 3200 3200 - 180 200 230 440 950 820 530] - -const NO2 = SA_F32[ - 120.0 130 160 480 2900 2700 2300 - 1900 1900 1900 2000 2700 2500 2200 - 1700 1700 1800 1900 2400 2300 2000 - 3900 4000 4100 4500 9999 9999 9999 - 270 290 350 850 2500 2300 2000] - -const H2O2 = SA_F32[ - 90.0 90 110 250 640 90 80 - 400 430 480 650 1100 90 90 - 370 390 430 550 840 90 80 - 400 430 470 620 1000 1000 1000 - 160 170 200 370 750 90 80] - -const ALD = SA_F32[ - 330.0 340 370 800 9999 9999 9999 - 9999 9999 9999 9999 9999 9999 9999 - 9999 9999 9999 9999 9999 9999 9999 - 9999 9999 9999 9999 9999 9999 9999 - 520 550 630 1700 9999 9999 9999] - -const HCHO = SA_F32[ - 100.0 110 140 450 6700 1400 1400 - 8700 8700 8700 8700 8700 1400 1400 - 8300 8300 8300 8300 8400 1400 1400 - 2900 2900 2900 2900 2900 2900 2900 - 250 270 340 1000 7500 1400 1400] - -const OP = SA_F32[ - 120.0 130 160 480 2800 2500 2200 - 1900 1900 1900 2000 2700 2400 2000 - 1700 1700 1800 1800 2400 2100 1900 - 3700 3700 3800 4200 8600 8600 8600 - 270 290 350 850 2500 2200 1900] - -const PAA = SA_F32[ - 150.0 160 200 580 2800 2400 2000 - 1900 1900 1900 2000 2700 2200 1900 - 1700 1700 1700 1800 2400 2000 1800 - 3400 3400 3500 3800 7200 7200 7200 - 330 350 420 960 2400 2100 1800] - -const ORA = SA_F32[ - 30.0 30 30 40 50 10 10 - 140 140 150 170 190 10 10 - 130 140 140 160 180 10 10 - 310 340 390 550 910 910 910 - 60 60 70 80 90 10 10] - -const NH3 = SA_F32[ - 80.0 80 100 320 2700 430 430 - 3400 3400 3400 3400 3400 440 440 - 3000 3000 3000 3000 3100 430 430 - 1500 1500 1500 1500 1500 1500 1500 - 180 200 240 680 2800 430 430] -const PAN = SA_F32[ - 190.0 210 250 700 2900 2700 2300 - 1900 1900 1900 2000 2700 2500 2200 - 1700 1700 1800 1900 2400 2300 2000 - 3900 4000 4100 4500 9999 9999 9999 - 410 430 510 1100 2500 2300 2000] -const HNO2 = SA_F32[ - 110.0 120 140 330 950 90 90 - 1000 1000 1000 1100 1400 90 90 - 860 860 870 910 1100 90 90 - 820 830 830 870 1000 1000 1000 - 220 240 280 530 1000 90 90] - -function different(a::Float64, b::Float32) where T<:AbstractFloat - c = abs(a - b) - return c/b > 0.1 && c >= 11.0 -end - -function TestWesely() - iLandUse = 4 # deciduous forest - Ts = [25.0, 10, 2, 0, 10] # Surface Temperature [C] - Garr = [800.0, 500, 300, 100, 0] # Solar radiation [W m-2] - θ = 0.0 # Slope [radians] - - polNames = [ - "SO2", "O3", "NO2", "H2O2", - "ALD", "HCHO", "OP", "PAA", - "ORA", "NH3", "PAN", "HNO2"] - - testData = [SO2, O3, NO2, H2O2, ALD, HCHO, OP, PAA, ORA, NH3, PAN, HNO2] - gasData = [ - So2Data, O3Data, No2Data, - H2o2Data, AldData, HchoData, - OpData, PaaData, OraData, - Nh3Data, PanData, Hno2Data] - - for i in 1:12 - pol = polNames[i] - polData = testData[i] - isSO2, isO3 = false, false - if pol == "SO2" - isSO2 = true - end - if pol == "O3" - isO3 = true - end - for iSeason in 1:5 - for ig in 1:5 - G = Garr[ig] - r_c = SurfaceResistance(gasData[i], G, Ts[iSeason], θ, - iSeason, iLandUse, false, false, isSO2, isO3) - if different(r_c, polData[iSeason, ig]) - println(pol, iSeason, G, r_c, polData[iSeason, ig]) - return false - end - end - r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, - iSeason, iLandUse, false, true, isSO2, isO3) # dew - if different(r_c, polData[iSeason, 6]) - println(pol, iSeason, "dew", r_c, polData[iSeason, 6]) - return false; - end - r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, - iSeason, iLandUse, true, false, isSO2, isO3) # rain - if different(r_c, polData[iSeason, 7]) - println(pol, iSeason, "rain", r_c, polData[iSeason, 7]) - return false; - end - end - end - return true -end +using Test, SafeTestsets @testset "DepositionMTK.jl" begin - @test TestWesely() == true - @test r_s(1.0, 1.0, 1, 1, true) ≈ 3.0772306344553016e30 - @test r_dc(1.0, 1.0) ≈ 9.181727363545544 - @test r_mx(1.0, 1.0) ≈ 0.009999966666777778 - @test r_smx(1.0, 1.0, 1.0) ≈ 2.0 - @test r_lux(1.0, 1.0, 1, 1, true, false, true, false) ≈ 50 - @test r_clx(1.0, 1.0, 1, 1) ≈ 9.999899563027895e24 - @test r_gsx(1.0, 1.0, 1, 1) ≈ 299.9977500168749 - @test SurfaceResistance(So2Data, 1.0, 1.0, 1.0, 1, 1, true, true, true, false) ≈ 45.45454545454546 + @safetestset "drydep_test.jl" begin + include("drydep_test.jl") + end + + @safetestset "wesely1989_test.jl" begin + include("wesely1989_test.jl") + end end diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl new file mode 100644 index 00000000..e0df1305 --- /dev/null +++ b/test/wesely1989_test.jl @@ -0,0 +1,163 @@ +#using DepositionMTK +#using Test +#using StaticArrays + +# include("../src/wesley1989.jl") + +# Results from Wesely (1989) table 3; updated to values from Walmsley (1996) table 1. +using DepositionMTK +using Test + +const SO2 = SA_F32[ + 130.0 140 160 380 1000 100 1200 + 1400 1400 1400 1400 1500 100 1300 + 1100 1100 1100 1100 1200 90 1000 + 1000 1000 1000 1000 1100 1100 1100 + 270 290 330 620 1100 90 1000] + +const O3 = SA_F32[ + 100.0 110 130 320 960 960 580 + 430 470 520 710 1300 950 580 + 390 420 460 610 960 770 510 + 560 620 710 1100 3200 3200 3200 + 180 200 230 440 950 820 530] + +const NO2 = SA_F32[ + 120.0 130 160 480 2900 2700 2300 + 1900 1900 1900 2000 2700 2500 2200 + 1700 1700 1800 1900 2400 2300 2000 + 3900 4000 4100 4500 9999 9999 9999 + 270 290 350 850 2500 2300 2000] + +const H2O2 = SA_F32[ + 90.0 90 110 250 640 90 80 + 400 430 480 650 1100 90 90 + 370 390 430 550 840 90 80 + 400 430 470 620 1000 1000 1000 + 160 170 200 370 750 90 80] + +const ALD = SA_F32[ + 330.0 340 370 800 9999 9999 9999 + 9999 9999 9999 9999 9999 9999 9999 + 9999 9999 9999 9999 9999 9999 9999 + 9999 9999 9999 9999 9999 9999 9999 + 520 550 630 1700 9999 9999 9999] + +const HCHO = SA_F32[ + 100.0 110 140 450 6700 1400 1400 + 8700 8700 8700 8700 8700 1400 1400 + 8300 8300 8300 8300 8400 1400 1400 + 2900 2900 2900 2900 2900 2900 2900 + 250 270 340 1000 7500 1400 1400] + +const OP = SA_F32[ + 120.0 130 160 480 2800 2500 2200 + 1900 1900 1900 2000 2700 2400 2000 + 1700 1700 1800 1800 2400 2100 1900 + 3700 3700 3800 4200 8600 8600 8600 + 270 290 350 850 2500 2200 1900] + +const PAA = SA_F32[ + 150.0 160 200 580 2800 2400 2000 + 1900 1900 1900 2000 2700 2200 1900 + 1700 1700 1700 1800 2400 2000 1800 + 3400 3400 3500 3800 7200 7200 7200 + 330 350 420 960 2400 2100 1800] + +const ORA = SA_F32[ + 30.0 30 30 40 50 10 10 + 140 140 150 170 190 10 10 + 130 140 140 160 180 10 10 + 310 340 390 550 910 910 910 + 60 60 70 80 90 10 10] + +const NH3 = SA_F32[ + 80.0 80 100 320 2700 430 430 + 3400 3400 3400 3400 3400 440 440 + 3000 3000 3000 3000 3100 430 430 + 1500 1500 1500 1500 1500 1500 1500 + 180 200 240 680 2800 430 430] +const PAN = SA_F32[ + 190.0 210 250 700 2900 2700 2300 + 1900 1900 1900 2000 2700 2500 2200 + 1700 1700 1800 1900 2400 2300 2000 + 3900 4000 4100 4500 9999 9999 9999 + 410 430 510 1100 2500 2300 2000] +const HNO2 = SA_F32[ + 110.0 120 140 330 950 90 90 + 1000 1000 1000 1100 1400 90 90 + 860 860 870 910 1100 90 90 + 820 830 830 870 1000 1000 1000 + 220 240 280 530 1000 90 90] + +function different(a::Float64, b::Float32) where T<:AbstractFloat + c = abs(a - b) + return c/b > 0.1 && c >= 11.0 +end + +function TestWesely() + iLandUse = 4 # deciduous forest + Ts = [25.0, 10, 2, 0, 10] # Surface Temperature [C] + Garr = [800.0, 500, 300, 100, 0] # Solar radiation [W m-2] + θ = 0.0 # Slope [radians] + + polNames = [ + "SO2", "O3", "NO2", "H2O2", + "ALD", "HCHO", "OP", "PAA", + "ORA", "NH3", "PAN", "HNO2"] + + testData = [SO2, O3, NO2, H2O2, ALD, HCHO, OP, PAA, ORA, NH3, PAN, HNO2] + gasData = [ + So2Data, O3Data, No2Data, + H2o2Data, AldData, HchoData, + OpData, PaaData, OraData, + Nh3Data, PanData, Hno2Data] + + for i in 1:12 + pol = polNames[i] + polData = testData[i] + isSO2, isO3 = false, false + if pol == "SO2" + isSO2 = true + end + if pol == "O3" + isO3 = true + end + for iSeason in 1:5 + for ig in 1:5 + G = Garr[ig] + r_c = SurfaceResistance(gasData[i], G, Ts[iSeason], θ, + iSeason, iLandUse, false, false, isSO2, isO3) + if different(r_c, polData[iSeason, ig]) + println(pol, iSeason, G, r_c, polData[iSeason, ig]) + return false + end + end + r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, + iSeason, iLandUse, false, true, isSO2, isO3) # dew + if different(r_c, polData[iSeason, 6]) + println(pol, iSeason, "dew", r_c, polData[iSeason, 6]) + return false; + end + r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, + iSeason, iLandUse, true, false, isSO2, isO3) # rain + if different(r_c, polData[iSeason, 7]) + println(pol, iSeason, "rain", r_c, polData[iSeason, 7]) + return false; + end + end + end + return true +end + +@testset "wesley1989.jl" begin + @test TestWesely() == true + @test r_s(1.0, 1.0, 1, 1, true) ≈ 3.0772306344553016e30 + @test r_dc(1.0, 1.0) ≈ 9.181727363545544 + @test r_mx(1.0, 1.0) ≈ 0.009999966666777778 + @test r_smx(1.0, 1.0, 1.0) ≈ 2.0 + @test r_lux(1.0, 1.0, 1, 1, true, false, true, false) ≈ 50 + @test r_clx(1.0, 1.0, 1, 1) ≈ 9.999899563027895e24 + @test r_gsx(1.0, 1.0, 1, 1) ≈ 299.9977500168749 + @test SurfaceResistance(So2Data, 1.0, 1.0, 1.0, 1, 1, true, true, true, false) ≈ 45.45454545454546 +end From 0ad2e5b1b2826d8acfb46755d4a2660bc132b1c6 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 18:47:27 -0500 Subject: [PATCH 03/32] Modify Project.toml --- Project.toml | 8 +++++++- src/dry_deposition.jl | 4 ---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index b8cb1c2b..324d19e5 100644 --- a/Project.toml +++ b/Project.toml @@ -8,12 +8,18 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] +EarthSciMLBase = "0.4" +SafeTestsets = "0.0.1" +StaticArrays = "1" +Unitful = "1" julia = "1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" [targets] -test = ["Test"] +test = ["Test","SafeTestsets"] diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 5b76abb5..f94b3059 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,7 +1,3 @@ -using Unitful -using ModelingToolkit -using StaticArrays -include("wesley1989.jl") export ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle g = 9.81u"m*s^-2" # gravitational acceleration [m/s2] From 07d427729c027fc7daad3a364b7ff96e058dcdb2 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 18:48:58 -0500 Subject: [PATCH 04/32] Modifying Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 324d19e5..e149fa6d 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] -EarthSciMLBase = "0.4" SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" From 1fb48dbde2433ab1d441f17627ecff4f5cc189b0 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 18:51:57 -0500 Subject: [PATCH 05/32] remove ModelingToolkit dependency --- src/DepositionMTK.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index a45e0f33..2bd5f98d 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -1,7 +1,6 @@ module DepositionMTK using Unitful -using ModelingToolkit using StaticArrays include("wesley1989.jl") From 376ab162ea81b56472ed10389d3fdd2b166461fb Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 18 May 2023 19:06:29 -0500 Subject: [PATCH 06/32] Add package for test files --- test/drydep_test.jl | 2 +- test/wesely1989_test.jl | 12 ++++-------- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/test/drydep_test.jl b/test/drydep_test.jl index 10b47d35..c048ddc5 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -1,5 +1,5 @@ using DepositionMTK -using Test +using Test,Unitful @testset "mfp" begin @test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl index e0df1305..b83f143d 100644 --- a/test/wesely1989_test.jl +++ b/test/wesely1989_test.jl @@ -1,13 +1,9 @@ -#using DepositionMTK -#using Test -#using StaticArrays - -# include("../src/wesley1989.jl") - -# Results from Wesely (1989) table 3; updated to values from Walmsley (1996) table 1. using DepositionMTK -using Test +using Test,StaticArrays +""" +Results from Wesely (1989) table 3; updated to values from Walmsley (1996) table 1. +""" const SO2 = SA_F32[ 130.0 140 160 380 1000 100 1200 1400 1400 1400 1400 1500 100 1300 From b2f036dcbc729b53420a48160c56e4a4dc9543a2 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 3 Jun 2023 11:12:06 +0800 Subject: [PATCH 07/32] Add wet deposition --- src/DepositionMTK.jl | 1 + src/wet_deposition.jl | 40 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 4 ++++ test/wetdep_test.jl | 8 ++++++++ 4 files changed, 53 insertions(+) create mode 100644 src/wet_deposition.jl create mode 100644 test/wetdep_test.jl diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index 2bd5f98d..0e3a43cc 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -5,5 +5,6 @@ using StaticArrays include("wesley1989.jl") include("dry_deposition.jl") +include("wet_deposition.jl") end \ No newline at end of file diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl new file mode 100644 index 00000000..94a5977f --- /dev/null +++ b/src/wet_deposition.jl @@ -0,0 +1,40 @@ +export WetDeposition + +""" +Calculate wet deposition based on formulas at +www.emep.int/UniDoc/node12.html. +Inputs are fraction of grid cell covered by clouds (cloudFrac), +rain mixing ratio (qrain), air density (ρ_air [kg/m3]), +and fall distance (Δz [m]). +Outputs are wet deposition rates for PM2.5, SO2, and other gases +(wdParticle, wdSO2, and wdOtherGas [1/s]). +""" +function WetDeposition(cloudFrac, qrain, ρ_air, Δz) + A = 5.2u"m^3/kg/s" # m3 kg-1 s-1; Empirical coefficient + E = 0.1 # size-dependent collection efficiency of aerosols by the raindrops + wSubSO2 = 0.15 # sub-cloud scavanging ratio + wSubOther = 0.5 # sub-cloud scavanging ratio + wInSO2 = 0.3 # in-cloud scavanging ratio + wInParticle = 1.0 # in-cloud scavanging ratio + wInOther = 1.4 # in-cloud scavanging ratio + ρwater = 1000.0u"kg*m^-3" # kg/m3 + Vdr = 5.0u"m/s" # raindrop fall speed, m/s + + # precalculated constant combinations + AE = A * E + wSubSO2VdrPerρwater = wSubSO2 * Vdr / ρwater + wSubOtherVdrPerρwater = wSubOther * Vdr / ρwater + wInSO2VdrPerρwater = wInSO2 * Vdr / ρwater + wInParticleVdrPerρwater = wInParticle * Vdr / ρwater + wInOtherVdrPerρwater = wInOther * Vdr / ρwater + + # wdParticle (subcloud) = A * P / Vdr * E; P = QRAIN * Vdr * ρgas => wdParticle = A * QRAIN * ρgas * E + # wdGas (subcloud) = wSub * P / Δz / ρwater = wSub * QRAIN * Vdr * ρgas / Δz / ρwater + # wd (in-cloud) = wIn * P / Δz / ρwater = wIn * QRAIN * Vdr * ρgas / Δz / ρwater + + wdParticle = qrain * ρ_air * (AE + cloudFrac*(wInParticleVdrPerρwater/Δz)) + wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz + wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz + + return wdParticle, wdSO2, wdOtherGas +end diff --git a/test/runtests.jl b/test/runtests.jl index 3f0d4375..08c3e71d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,4 +9,8 @@ using Test, SafeTestsets @safetestset "wesely1989_test.jl" begin include("wesely1989_test.jl") end + + @safetestset "wetdep_test.jl" begin + includ("wetdep_test.jl") + end end diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl new file mode 100644 index 00000000..07f355b8 --- /dev/null +++ b/test/wetdep_test.jl @@ -0,0 +1,8 @@ +using DepositionMTK +using Test,Unitful + +@testset "unit" begin + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[1])==u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[2])==u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[3])==u"s^-1" +end \ No newline at end of file From 430e9c8cf5c1f727d6b9869c2bf2eaf78b637ceb Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 3 Jun 2023 11:44:22 +0800 Subject: [PATCH 08/32] Change Unit Test --- test/runtests.jl | 4 ++-- test/wetdep_test.jl | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 08c3e71d..86a240d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,6 @@ using Test, SafeTestsets end @safetestset "wetdep_test.jl" begin - includ("wetdep_test.jl") + include("wetdep_test.jl") end -end +end \ No newline at end of file diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index 07f355b8..a2bc1c92 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -2,7 +2,7 @@ using DepositionMTK using Test,Unitful @testset "unit" begin - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[1])==u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[2])==u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[3])==u"s^-1" -end \ No newline at end of file + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[1]) == u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[2]) == u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[3]) == u"s^-1" +end From 855146e202342b9e6fcaf6bbefe5ed20307774dd Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Fri, 16 Jun 2023 22:40:35 +0800 Subject: [PATCH 09/32] Add ODE module --- Project.toml | 8 +++++-- src/wet_deposition.jl | 51 +++++++++++++++++++++++++++++++++++++++---- test/wetdep_test.jl | 6 ++--- 3 files changed, 56 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index e149fa6d..82cae4f2 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,12 @@ authors = ["EarthSciML authors and contributors"] version = "0.1.0" [deps] +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" @@ -17,8 +21,8 @@ Unitful = "1" julia = "1" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test","SafeTestsets"] +test = ["Test", "SafeTestsets"] diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 94a5977f..25583067 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -1,4 +1,4 @@ -export WetDeposition +export WetDeposition, Wetdeposition """ Calculate wet deposition based on formulas at @@ -9,7 +9,7 @@ and fall distance (Δz [m]). Outputs are wet deposition rates for PM2.5, SO2, and other gases (wdParticle, wdSO2, and wdOtherGas [1/s]). """ -function WetDeposition(cloudFrac, qrain, ρ_air, Δz) +function WetDeposition(cloudFrac, qrain, ρ_air, Δz, output) A = 5.2u"m^3/kg/s" # m3 kg-1 s-1; Empirical coefficient E = 0.1 # size-dependent collection efficiency of aerosols by the raindrops wSubSO2 = 0.15 # sub-cloud scavanging ratio @@ -32,9 +32,52 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) # wdGas (subcloud) = wSub * P / Δz / ρwater = wSub * QRAIN * Vdr * ρgas / Δz / ρwater # wd (in-cloud) = wIn * P / Δz / ρwater = wIn * QRAIN * Vdr * ρgas / Δz / ρwater - wdParticle = qrain * ρ_air * (AE + cloudFrac*(wInParticleVdrPerρwater/Δz)) + wdParticle = qrain * ρ_air * (AE + cloudFrac * (wInParticleVdrPerρwater / Δz)) wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz - return wdParticle, wdSO2, wdOtherGas + if output == 1 + return wdParticle + elseif output == 2 + return wdSO2 + else + return wdOtherGas + end end + +using EarthSciMLBase +using ModelingToolkit +using Unitful +# Add unit "ppb" to Unitful +module MyUnits +using Unitful +@unit ppb "ppb" Number 1 / 1000000000 false +end +Unitful.register(MyUnits) + +struct Wetdeposition <: EarthSciMLODESystem + sys::ODESystem + function Wetdeposition(t) + cloudFrac = 0.5 + qrain = 0.5 + ρ_air = 1.204u"kg*m^-3" + Δz = 200u"m" + @parameters k1 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 2) * 1u"s" [unit = u"s^-1"] + @parameters k2 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 3) * 1u"s" [unit = u"s^-1"] + @parameters t [unit = u"s"] + + D = Differential(t) + + @variables wet_SO2(t) = 2.0 [unit = u"ppb"] + @variables wet_other(t) = 10.0 [unit = u"ppb"] + + eqs = [ + D(wet_SO2) ~ -k1 * wet_SO2 + D(wet_other) ~ -k2 * wet_other + ] + + new(ODESystem(eqs, t, [wet_SO2, wet_other], [k1,k2]; name=:Wetdeposition)) + #ODESystem(eqs, t, [wet_SO2, wet_other], [k1, k2]; name=:WetDeposition) + end +end + diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index a2bc1c92..0b94f91c 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -2,7 +2,7 @@ using DepositionMTK using Test,Unitful @testset "unit" begin - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[1]) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[2]) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[3]) == u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",1)) == u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",2)) == u"s^-1" + @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",3)) == u"s^-1" end From 99ea0f3eb20c0bffcede9e489e41ee91bb7a0eb5 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Fri, 16 Jun 2023 23:22:54 +0800 Subject: [PATCH 10/32] Update of Wetdeposition --- src/wet_deposition.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 25583067..8014a156 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -57,27 +57,22 @@ Unitful.register(MyUnits) struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem - function Wetdeposition(t) - cloudFrac = 0.5 - qrain = 0.5 - ρ_air = 1.204u"kg*m^-3" - Δz = 200u"m" + function Wetdeposition(t, cloudFrac, qrain, ρ_air, Δz) @parameters k1 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 2) * 1u"s" [unit = u"s^-1"] @parameters k2 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 3) * 1u"s" [unit = u"s^-1"] @parameters t [unit = u"s"] D = Differential(t) - @variables wet_SO2(t) = 2.0 [unit = u"ppb"] - @variables wet_other(t) = 10.0 [unit = u"ppb"] + @variables SO2(t) [unit = u"ppb"] + @variables O3(t) [unit = u"ppb"] eqs = [ - D(wet_SO2) ~ -k1 * wet_SO2 - D(wet_other) ~ -k2 * wet_other + D(SO2) ~ -k1 * SO2 + D(O3) ~ -k2 * O3 ] - new(ODESystem(eqs, t, [wet_SO2, wet_other], [k1,k2]; name=:Wetdeposition)) - #ODESystem(eqs, t, [wet_SO2, wet_other], [k1, k2]; name=:WetDeposition) + new(ODESystem(eqs, t, [SO2, O3], [k1,k2]; name=:Wetdeposition)) end end From 02621778bab5a3f7da22246ba1769215ab61ed0c Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 17 Jun 2023 00:07:15 +0800 Subject: [PATCH 11/32] Delete @parameters t --- src/wet_deposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 8014a156..293c664c 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -60,7 +60,7 @@ struct Wetdeposition <: EarthSciMLODESystem function Wetdeposition(t, cloudFrac, qrain, ρ_air, Δz) @parameters k1 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 2) * 1u"s" [unit = u"s^-1"] @parameters k2 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 3) * 1u"s" [unit = u"s^-1"] - @parameters t [unit = u"s"] + #@parameters t [unit = u"s"] D = Differential(t) From 3a7721d2c8a0fcf781f6a3e642ee21f212f1517b Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 17 Jun 2023 00:15:00 +0800 Subject: [PATCH 12/32] update --- src/wet_deposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 293c664c..8014a156 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -60,7 +60,7 @@ struct Wetdeposition <: EarthSciMLODESystem function Wetdeposition(t, cloudFrac, qrain, ρ_air, Δz) @parameters k1 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 2) * 1u"s" [unit = u"s^-1"] @parameters k2 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 3) * 1u"s" [unit = u"s^-1"] - #@parameters t [unit = u"s"] + @parameters t [unit = u"s"] D = Differential(t) From 0dba85f2af0233f9f4436e2106f5c830703f5cea Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Wed, 21 Jun 2023 17:17:55 +0800 Subject: [PATCH 13/32] Change of @onstants --- src/DepositionMTK.jl | 2 ++ src/wet_deposition.jl | 46 +++++++++++++++++++++---------------------- test/wetdep_test.jl | 13 +++++++++--- 3 files changed, 34 insertions(+), 27 deletions(-) diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index 0e3a43cc..4fb00139 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -2,6 +2,8 @@ module DepositionMTK using Unitful using StaticArrays +using ModelingToolkit +using EarthSciMLBase include("wesley1989.jl") include("dry_deposition.jl") diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 8014a156..f92e22cc 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -1,4 +1,4 @@ -export WetDeposition, Wetdeposition +export WetDeposition, Wetdeposition, wd_defaults """ Calculate wet deposition based on formulas at @@ -9,19 +9,24 @@ and fall distance (Δz [m]). Outputs are wet deposition rates for PM2.5, SO2, and other gases (wdParticle, wdSO2, and wdOtherGas [1/s]). """ -function WetDeposition(cloudFrac, qrain, ρ_air, Δz, output) - A = 5.2u"m^3/kg/s" # m3 kg-1 s-1; Empirical coefficient + +@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient +@constants ρwater = 1000.0 [unit = u"kg*m^-3"] # kg/m3 +@constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s + +function WetDeposition(cloudFrac, qrain, ρ_air, Δz) + #@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient E = 0.1 # size-dependent collection efficiency of aerosols by the raindrops wSubSO2 = 0.15 # sub-cloud scavanging ratio wSubOther = 0.5 # sub-cloud scavanging ratio wInSO2 = 0.3 # in-cloud scavanging ratio wInParticle = 1.0 # in-cloud scavanging ratio wInOther = 1.4 # in-cloud scavanging ratio - ρwater = 1000.0u"kg*m^-3" # kg/m3 - Vdr = 5.0u"m/s" # raindrop fall speed, m/s + # @constants ρwater = 1000.0 [unit = u"kg*m^-3"] # kg/m3 + # @constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s # precalculated constant combinations - AE = A * E + AE = A_wd * E wSubSO2VdrPerρwater = wSubSO2 * Vdr / ρwater wSubOtherVdrPerρwater = wSubOther * Vdr / ρwater wInSO2VdrPerρwater = wInSO2 * Vdr / ρwater @@ -36,18 +41,10 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz, output) wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz - if output == 1 - return wdParticle - elseif output == 2 - return wdSO2 - else - return wdOtherGas - end + return wdParticle, wdSO2, wdOtherGas end +wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0] -using EarthSciMLBase -using ModelingToolkit -using Unitful # Add unit "ppb" to Unitful module MyUnits using Unitful @@ -57,9 +54,11 @@ Unitful.register(MyUnits) struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem - function Wetdeposition(t, cloudFrac, qrain, ρ_air, Δz) - @parameters k1 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 2) * 1u"s" [unit = u"s^-1"] - @parameters k2 = WetDeposition(cloudFrac, qrain, ρ_air, Δz, 3) * 1u"s" [unit = u"s^-1"] + function Wetdeposition(t) + @parameters cloudFrac = 0.5 + @parameters qrain = 0.5 + @parameters ρ_air = 1.204 [unit = u"kg*m^-3"] + @parameters Δz = 200 [unit = u"m"] @parameters t [unit = u"s"] D = Differential(t) @@ -68,11 +67,10 @@ struct Wetdeposition <: EarthSciMLODESystem @variables O3(t) [unit = u"ppb"] eqs = [ - D(SO2) ~ -k1 * SO2 - D(O3) ~ -k2 * O3 + D(SO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2 + D(O3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * O3 ] - new(ODESystem(eqs, t, [SO2, O3], [k1,k2]; name=:Wetdeposition)) + new(ODESystem(eqs, t, [SO2, O3], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) end -end - +end \ No newline at end of file diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index 0b94f91c..661ce352 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -1,8 +1,15 @@ using DepositionMTK using Test,Unitful +@parameters cloudFrac = 0.5 +@parameters qrain = 0.5 +@parameters ρ_air = 1.204 [unit = u"kg*m^-3"] +@parameters Δz = 200 [unit = u"m"] + @testset "unit" begin - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",1)) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",2)) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m",3)) == u"s^-1" + @test substitute(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[1], Dict(cloudFrac => 0.5, qrain => 0.5, ρ_air => 1.204, Δz => 200, wd_defaults...)) ≈ 0.313047525 + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[1]) == u"s^-1" + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2]) == u"s^-1" + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3]) == u"s^-1" end + From bc99b8cbe1ef46a32e7ae8a3dbc847414a607e3c Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Wed, 21 Jun 2023 19:28:22 +0800 Subject: [PATCH 14/32] delete in function parameter t --- src/wet_deposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index f92e22cc..d70a5a2e 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -59,7 +59,7 @@ struct Wetdeposition <: EarthSciMLODESystem @parameters qrain = 0.5 @parameters ρ_air = 1.204 [unit = u"kg*m^-3"] @parameters Δz = 200 [unit = u"m"] - @parameters t [unit = u"s"] + #@parameters t [unit = u"s"] D = Differential(t) From d02045609a66417c3fa33c204502940b3910d098 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Wed, 21 Jun 2023 19:37:02 +0800 Subject: [PATCH 15/32] Change back to add parameter t in function --- src/wet_deposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index d70a5a2e..f92e22cc 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -59,7 +59,7 @@ struct Wetdeposition <: EarthSciMLODESystem @parameters qrain = 0.5 @parameters ρ_air = 1.204 [unit = u"kg*m^-3"] @parameters Δz = 200 [unit = u"m"] - #@parameters t [unit = u"s"] + @parameters t [unit = u"s"] D = Differential(t) From e4b351bd51ed94d438e1de7e857a361cb0dc2a8e Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 18:58:54 +0800 Subject: [PATCH 16/32] Change functions and unit test to MTK --- Project.toml | 1 + src/DepositionMTK.jl | 1 + src/dry_deposition.jl | 76 ++++++++++++++++---------------- src/wesley1989.jl | 32 ++++++++------ test/drydep_test.jl | 97 ++++++++++++++++++++++++----------------- test/wesely1989_test.jl | 2 +- test/wetdep_test.jl | 2 +- 7 files changed, 118 insertions(+), 93 deletions(-) diff --git a/Project.toml b/Project.toml index 82cae4f2..8e7e5b7b 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" +IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index 4fb00139..6e461ca8 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -4,6 +4,7 @@ using Unitful using StaticArrays using ModelingToolkit using EarthSciMLBase +using IfElse include("wesley1989.jl") include("dry_deposition.jl") diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index f94b3059..6278ece1 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,10 +1,10 @@ -export ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle +export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle -g = 9.81u"m*s^-2" # gravitational acceleration [m/s2] -κ = 0.4 # von Karman constant -k = 1.3806488e-23u"m^2*kg*s^-2/K" # Boltzmann constant -M_air = 28.97e-3u"kg/mol" # molecular weight of air -R = 8.3144621u"J/K/mol" # Gas constant +@constants g = 9.81 [unit = u"m*s^-2"] # gravitational acceleration [m/s2] +@constants κ = 0.4 # von Karman constant +@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K"] # Boltzmann constant +@constants M_air = 28.97e-3 [unit = u"kg/mol"] # molecular weight of air +@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1"] # Gas constant """ Function Ra calculates aerodynamic resistance to dry deposition @@ -12,34 +12,26 @@ where z is the top of the surface layer [m], z₀ is the roughness length [m], u Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 & 19.14. """ +@constants unit_m = 1 [unit = u"m"] function ra(z, z₀, u_star, L) - if L == 0u"m" - ζ = 0 - ζ₀= 0 - else - ζ = z/L - ζ₀= z₀/L - end - print(ζ) - if 0 < ζ < 1 - rₐ = 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)) - elseif ζ == 0 - rₐ = 1/(κ*u_star)*log(z/z₀) - elseif -1 < ζ < 0 - η₀=(1-15*ζ₀)^1/4 - η =(1-15*ζ)^1/4 - rₐ = 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))] - else - print("wrong rₐ") - end + ζ = IfElse.ifelse((L/unit_m == 0), 0, z/L) + ζ₀ = IfElse.ifelse((L/unit_m == 0), 0, z₀/L) + rₐ_1 = IfElse.ifelse((0 < ζ) & (ζ < 1), 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)), 1/(κ*u_star)*log(z/z₀)) + η₀=(1-15*ζ₀)^1/4 + η =(1-15*ζ)^1/4 + rₐ = IfElse.ifelse((-1 < ζ) & (ζ < 0), 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))][1], rₐ_1) return rₐ end """ Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. """ + +@constants unit_T = 1 [unit = u"K"] +@constants unit_convert_mu = 1 [unit = u"kg/m/s"] function mu(T) - return (1.8e-5*(T/298u"K")^0.85)u"kg/m/s" + return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu + #return (1.8e-5*(T/T₀)^0.85)*unit_convert_mu end """ @@ -66,12 +58,15 @@ Function vs calculates the terminal setting velocity of a particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. From equation 9.42 in Seinfeld and Pandis (2006). """ +# Particle diameter Dₚ greater than 20um; Stokes settling no longer applies. +@constants unit_v = 1 [unit = u"m/s"] function vs(Dₚ, ρₚ, Cc, μ) - if Dₚ > 20.e-6u"m" - print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.") - else - return Dₚ^2*ρₚ*g*Cc/(18*μ) - end + IfElse.ifelse((Dₚ > 20.e-6*unit_m), 99999999*unit_v, Dₚ^2*ρₚ*g*Cc/(18*μ)) + # if Dₚ > 20.e-6u"m" + # print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.") + # else + # return Dₚ^2*ρₚ*g*Cc/(18*μ) + # end end """ @@ -87,9 +82,12 @@ end Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K] using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf """ + +@constants T_unitless = 1 [unit = u"K^-1"] +@constants unit_dH2O = 1 [unit = u"m^2/s"] function dH2O(T) - T_unitless = T*1u"K^-1" - return (-2.775e-6 + 4.479e-8*T_unitless + 1.656e-10*T_unitless^2)u"m^2/s" + #T_unitless = T*1u"K^-1" + return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O end """ @@ -178,7 +176,7 @@ From Seinfeld and Pandis (2006) equation 19.27. function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) α = α_table[iLandUse] γ = γ_table[iLandUse] - A = (A_table[iSeason,iLandUse]*10^(-3))u"m" + A = (A_table[iSeason,iLandUse]*10^(-3))*unit_m R1 = exp(-St^0.5) term_1 = Sc^(-γ) term_2 = (St/(α+St))^2 @@ -195,13 +193,16 @@ irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. Based on Seinfeld and Pandis (2006) equation 19.2. """ + +@constants G_unitless = 1 [unit = u"m^2/W"] +@constants Rc_unit = 1 [unit = u"s/m"] function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s] Sc = sc(μ,ρA, Dg) Rb = RbGas(Sc, u_star) - Rc = SurfaceResistance(gasData, G/u"W*m^-2", (Ts/u"K"-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m" + Rc = SurfaceResistance(gasData, G*G_unitless, (Ts*T_unitless-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)*Rc_unit return 1/(Ra+Rb+Rc) end @@ -218,9 +219,9 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: Cc = cc(Dp, Ts, P, μ) Vs = vs(Dp, ρParticle, Cc, μ) if iLandUse == 4 # dessert - St = stSmooth(Vs, u_star, μ, ρA) + St = stSmooth(Vs, u_star, μ, ρA) else - St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))u"m") + St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))*unit_m) end D = dParticle(Ts, P, Dp, Cc, μ) Sc = sc(μ, ρA, D) @@ -228,3 +229,4 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: return 1/(Ra+Rb+Ra*Rb*Vs)+Vs end +defaults = [g => 9.81, κ => 0.4, k => 1.3806488e-23, M_air => 28.97e-3, R => 8.3144621, unit_T => 1, unit_convert_mu => 1, T_unitless => 1, unit_dH2O => 1, unit_m => 1, G_unitless => 1, Rc_unit => 1, unit_v => 1] \ No newline at end of file diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 8835ea43..432ee806 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -88,12 +88,14 @@ # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) rs = 0.0 - if Ts >= 39.9 || Ts <= 0.1 - rs = inf - else - rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * - (400.0 * 1.0 / (Ts * (40.0 - Ts))) - end + rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * + (400.0 * 1.0 / (Ts * (40.0 - Ts)))) + # if Ts >= 39.9 || Ts <= 0.1 + # rs = inf + # else + # rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * + # (400.0 * 1.0 / (Ts * (40.0 - Ts))) + # end # Adjust for dew and rain (from "Effects of dew and rain" section). if rainOrDew rs *= 3 @@ -112,7 +114,7 @@ end #Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]). - function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat + function r_smx(r_s, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat return r_s * Dh2oPerDx + r_mx end @@ -193,12 +195,16 @@ rac = r_ac[iSeason, iLandUse] # Correction for cold temperatures from page 4 column 1. - if Ts < 0.0 - correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark - rlux += correction - rclx += correction - rgsx += correction - end + correction = IfElse.ifelse((Ts < 0.0), 1000.0 * exp(-Ts-4), 0) # [s m-1] #mark + rlux += correction + rclx += correction + rgsx += correction + # if Ts < 0.0 + # correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark + # rlux += correction + # rclx += correction + # rgsx += correction + # end r_c = 1.0 / (1.0/(rsmx) + 1.0/rlux + 1.0/(rdc+rclx) + 1.0/(rac+rgsx)) r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces. diff --git a/test/drydep_test.jl b/test/drydep_test.jl index c048ddc5..0d3ae8af 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -1,77 +1,92 @@ using DepositionMTK -using Test,Unitful +using Test,Unitful, ModelingToolkit + +begin + @parameters T [unit = u"K"] + @parameters P [unit = u"kg*m^-1*s^-2"] # 1Pa = 1kg/m/s^-2 + @parameters μ [unit = u"kg/m/s"] + @parameters z [unit = u"m"] + @parameters z₀ [unit = u"m"] + @parameters u_star [unit = u"m/s"] + @parameters L [unit = u"m"] + @parameters Dp [unit = u"m"] + @parameters ρParticle [unit = u"kg*m^-3"] + @parameters ρA [unit = u"kg*m^-3"] + @parameters G [unit = u"W*m^-2"] + @parameters θ +end @testset "mfp" begin - @test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" + @test substitute(mfp(T,P,μ), Dict(T => 298, P => 101300, μ => 1.8e-5, defaults...)) ≈ 6.512893276888993e-8 + @test ModelingToolkit.get_unit(mfp(T,P,μ)) == u"m" + #@test unit(ModelingToolkit.get_unit(mfp(T,P,μ))*1 - 1u"m") == u"m" + #@test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" end @testset "unit" begin - @test unit(dH2O(300u"K"))==u"m^2/s" - @test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" - @test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" + @test ModelingToolkit.get_unit(dH2O(T)) == u"m^2/s" + #@test unit(dH2O(300u"K"))==u"m^2/s" + @test ModelingToolkit.get_unit(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 1)) == u"m/s" + #@test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" + @test ModelingToolkit.get_unit(DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, 1, 1, false, false, true, false)) == u"m/s" + #@test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" end @testset "viscosity" begin - T = [100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 1200, 1300, 1400, 1500, 1600] - μ_list = [0.6924, 1.0283, 1.3289, 1.488, 1.983, 2.075, 2.286, 2.484, 2.671, 2.848, 3.018, 3.177, 3.332, 3.481, 3.625, 3.765, 3.899, 4.023, 4.152, 4.44, 4.69, 4.93, 5.17, 5.4, 5.63] + T_ = [275, 300, 325, 350, 375, 400] + μ_list = [1.725, 1.846, 1.962, 2.075, 2.181, 2.286].*10^-5 μ_test = [] - for i in 1:25 - push!(μ_test, (mu((T[i])u"K"))u"m*s/kg") + for i in 1:6 + push!(μ_test, substitute(mu(T), Dict(T => T_[i], defaults...))) + end + for i in 1:6 + @test (μ_test[i] - μ_list[i])/ μ_list[i] < 0.01 end - @test μ_list ≈ μ_list rtol=1e-2 end @testset "Cc" begin - Dp = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] + Dp_ = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] Cc_list = [216, 108, 43.6, 22.2, 11.4, 4.95, 2.85, 1.865, 1.326, 1.164, 1.082, 1.032, 1.016, 1.008, 1.003, 1.0016] Cc_test = [] for i in 1:16 - push!(Cc_test, cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K"))) + push!(Cc_test, substitute(cc(Dp,T,P,μ), Dict(Dp => Dp_[i]*1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, defaults...)), defaults...))) + end + for i in 1: 16 + @test (Cc_test[i] - Cc_list[i])/ Cc_list[i] < 0.03 end - @test Cc_test ≈ Cc_list rtol=1e-2 end +@parameters Cc @testset "Vs" begin - Dp = [0.01, 0.1, 1, 10] - Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 + Dp_ = [0.01, 0.1, 1, 10] + Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 # convert to m/s Vs_test = [] + Cc_list = [] for i in 1:4 - push!(Vs_test, (vs((Dp[i]*1e-6)u"m", 1000u"kg*m^-3", cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K")), mu(298u"K")))u"s/m") - end - @test Vs_test ≈ Vs_list rtol=0.1 + push!(Cc_list, substitute(cc(Dp,T,P,μ), Dict(Dp => Dp_[i]*1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, defaults...)), defaults...))) + push!(Vs_test, substitute(vs(Dp, ρParticle, Cc, μ), Dict(Dp => Dp_[i]*1e-6, ρParticle => 1000, Cc => Cc_list[i], μ => 1.836522217711828e-5, defaults...))) + end + for i in 1:4 + @test (Vs_test[i]-Vs_list[i])/Vs_list[i] < 1 + end end @testset "DryDepGas" begin - z = 50u"m" - z₀ = 0.04u"m" - u_star = 0.44u"m/s" - L = 0u"m" - T = 298u"K" - ρA = 1.2u"kg*m^-3" - G = 300u"W*m^-2" - θ = 0 - vd_true = 0.003u"m/s" - @test DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, 1, 10, false, false, false, false) ≈ vd_true rtol=0.33 + vd_true = 0.03 # m/s + @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, 0, 1, 10, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, defaults...)) - vd_true)/ vd_true < 0.33 end @testset "DryDepParticle" begin - z = 20u"m" - z₀ = 0.02u"m" - u_star = 0.44u"m/s" - L = 0u"m" - T = 298u"K" - P = 101325u"Pa" - ρA = 1.2u"kg*m^-3" - ρParticle = 1000u"kg*m^-3" - Dp = [1.e-8, 1.e-7, 1.e-6, 1.e-5] - vd_true = [0.5, 0.012, 0.02, 0.6] # [cm/s] + Dp_ = [1.e-8, 1.e-7, 1.e-6, 1.e-5] + vd_true = [0.5, 0.012, 0.02, 0.6] ./100 # [m/s] vd_list = [] for i in 1:4 - push!(vd_list, (DryDepParticle(z, z₀, u_star, L, (Dp[i])u"m", T, P, ρParticle, ρA, 1, 4))*100u"s/m") + push!(vd_list, substitute(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 4), Dict(z => 20, z₀ => 0.02, u_star => 0.44, L => 0, T => 298, P=> 101325, ρA => 1.2, ρParticle => 1000, Dp => Dp_[i], defaults...))) + end + for i in 1:4 + @test vd_list[i]-vd_true[i] < 0.015 end - @test vd_list ≈ vd_true rtol=0.8 # 80% difference, not ideal end - diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl index b83f143d..d0639507 100644 --- a/test/wesely1989_test.jl +++ b/test/wesely1989_test.jl @@ -86,7 +86,7 @@ const HNO2 = SA_F32[ 820 830 830 870 1000 1000 1000 220 240 280 530 1000 90 90] -function different(a::Float64, b::Float32) where T<:AbstractFloat +function different(a::Float64, b::Float32) #where T<:AbstractFloat c = abs(a - b) return c/b > 0.1 && c >= 11.0 end diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index 661ce352..a8ed2c8f 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -1,5 +1,5 @@ using DepositionMTK -using Test,Unitful +using Test,Unitful, ModelingToolkit @parameters cloudFrac = 0.5 @parameters qrain = 0.5 From 93433590231990a9e799e26bd7fc7abf8f9a9822 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 19:45:37 +0800 Subject: [PATCH 17/32] Change of Project.toml --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index 8e7e5b7b..e1b9f701 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,9 @@ SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" julia = "1" +EarthSciMLBase = "0.4" +ModelingToolkit = "8" + [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" From 0bedc50ee5ef0b362d8dc8115705b0d6b24030e4 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 20:19:01 +0800 Subject: [PATCH 18/32] Try --- Project.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e1b9f701..41a26e8a 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,8 @@ ModelingToolkit = "8" [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" [targets] -test = ["Test", "SafeTestsets"] +test = ["Test", "SafeTestsets", "Unitful", "ModelingToolkit"] From 62e2e9ec2c603cb48b4fdb9d5e168ccd3df69536 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 20:37:57 +0800 Subject: [PATCH 19/32] Why jobs fail --- Project.toml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 41a26e8a..e1b9f701 100644 --- a/Project.toml +++ b/Project.toml @@ -27,8 +27,6 @@ ModelingToolkit = "8" [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" [targets] -test = ["Test", "SafeTestsets", "Unitful", "ModelingToolkit"] +test = ["Test", "SafeTestsets"] From c3907001a085f1d5c163b0d83c1ef178099711f2 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 20:48:03 +0800 Subject: [PATCH 20/32] Change Modelingtoolkit version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e1b9f701..2b0d7444 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ StaticArrays = "1" Unitful = "1" julia = "1" EarthSciMLBase = "0.4" -ModelingToolkit = "8" +ModelingToolkit = "8.0.0" [extras] From 353df372fadd8b00a18e9fd78994e78e4ee7ab14 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 22:53:26 +0800 Subject: [PATCH 21/32] Changde compat --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2b0d7444..624e12fe 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,6 @@ StaticArrays = "1" Unitful = "1" julia = "1" EarthSciMLBase = "0.4" -ModelingToolkit = "8.0.0" [extras] From 79fd6fc79b4f87a263ad1437b87dabe6334985bd Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Thu, 22 Jun 2023 22:56:23 +0800 Subject: [PATCH 22/32] Change of compat EarthSciMLBase --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 624e12fe..5b2d281e 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,6 @@ SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" julia = "1" -EarthSciMLBase = "0.4" [extras] From 70c5a3d5aa84e2bc6cd3bf5d1f94cc52c3f6aa10 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 1 Jul 2023 19:41:58 +0800 Subject: [PATCH 23/32] Add deposition EarthsciMl function and downgrade Catalyst to version 12 --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 5b2d281e..e1b9f701 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,8 @@ SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" julia = "1" +EarthSciMLBase = "0.4" +ModelingToolkit = "8" [extras] From 174ff8c91cac8e40c51db187059280f221c3ebe4 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 1 Jul 2023 19:59:34 +0800 Subject: [PATCH 24/32] add DrydepositionG function --- src/dry_deposition.jl | 51 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 6278ece1..08223526 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,4 +1,4 @@ -export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle +export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle, DrydepositionG @constants g = 9.81 [unit = u"m*s^-2"] # gravitational acceleration [m/s2] @constants κ = 0.4 # von Karman constant @@ -229,4 +229,51 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: return 1/(Ra+Rb+Ra*Rb*Vs)+Vs end -defaults = [g => 9.81, κ => 0.4, k => 1.3806488e-23, M_air => 28.97e-3, R => 8.3144621, unit_T => 1, unit_convert_mu => 1, T_unitless => 1, unit_dH2O => 1, unit_m => 1, G_unitless => 1, Rc_unit => 1, unit_v => 1] \ No newline at end of file +defaults = [g => 9.81, κ => 0.4, k => 1.3806488e-23, M_air => 28.97e-3, R => 8.3144621, unit_T => 1, unit_convert_mu => 1, T_unitless => 1, unit_dH2O => 1, unit_m => 1, G_unitless => 1, Rc_unit => 1, unit_v => 1] + +# Add unit "ppb" to Unitful +module MyUnits +using Unitful +@unit ppb "ppb" Number 1 / 1000000000 false +end +Unitful.register(MyUnits) + +struct DrydepositionG <: EarthSciMLODESystem + sys::ODESystem + function DrydepositionG(t) + iSeason = 1 + iLandUse = 10 + rain = false + dew = false + @parameters z = 50 [unit = u"m"] + @parameters z₀ = 0.04 [unit = u"m"] + @parameters u_star = 0.44 [unit = u"m/s"] + @parameters L = 0 [unit = u"m"] + @parameters ρA = 1.2 [unit = u"kg*m^-3"] + @parameters G = 300 [unit = u"W*m^-2"] + @parameters T = 298 [unit = u"K"] + @parameters θ = 0 + @parameters t [unit = u"s"] + + D = Differential(t) + + @variables SO2(t) [unit = u"ppb"] + @variables O3(t) [unit = u"ppb"] + @variables NO2(t) [unit = u"ppb"] + @variables NO(t) [unit = u"ppb"] + @variables H2O2(t) [unit = u"ppb"] + @variables CH2O(t) [unit = u"ppb"] + + eqs = [ + D(SO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, iSeason, iLandUse, rain, dew, true, false) / z * SO2 + D(O3) ~ -DryDepGas(z, z₀, u_star, L, ρA, O3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, true) / z * O3 + D(NO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO2 + D(NO) ~ -DryDepGas(z, z₀, u_star, L, ρA, NoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO + D(H2O2) ~ -DryDepGas(z, z₀, u_star, L, ρA, H2o2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * H2O2 + D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O + ] + + new(ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG)) + end +end + From 36efe8aa938cb8b4d76784458ee78eb8820e79c3 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sat, 1 Jul 2023 20:12:50 +0800 Subject: [PATCH 25/32] change from Unitful.MyUnit to myUnit in wet_deposition --- src/wet_deposition.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index f92e22cc..a9031b50 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -46,11 +46,11 @@ end wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0] # Add unit "ppb" to Unitful -module MyUnits +module myUnits using Unitful @unit ppb "ppb" Number 1 / 1000000000 false end -Unitful.register(MyUnits) +Unitful.register(myUnits) struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem From c762435cc234a23ff073177968cc1558a49cd25c Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sun, 2 Jul 2023 17:27:23 +0800 Subject: [PATCH 26/32] git add src/wet_deposition.jl --- .github/workflows/CI.yml | 4 ++-- src/dry_deposition.jl | 10 ++++++++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f3733deb..e498016c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,8 +10,8 @@ jobs: fail-fast: false matrix: version: - - '1.0' - - '1.6' + - '1.8' + - '1.9' - 'nightly' os: - ubuntu-latest diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 08223526..1bf175ea 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -238,6 +238,16 @@ using Unitful end Unitful.register(MyUnits) +""" +Description: This is a box model used to calculate the gas species concentration rate changed by dry deposition. +Build Drydeposition model (gas) +# Example +``` julia + @parameters t + d = DrydepositionG(t) +``` +""" + struct DrydepositionG <: EarthSciMLODESystem sys::ODESystem function DrydepositionG(t) From dd791fde8404d31f54d765bd1236bd914deb18fc Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sun, 2 Jul 2023 17:28:05 +0800 Subject: [PATCH 27/32] Add Example description for dry and wet deposition function and change Cl.yml julia version --- src/wet_deposition.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index a9031b50..c4926c34 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -52,6 +52,16 @@ using Unitful end Unitful.register(myUnits) +""" +Description: This is a box model used to calculate the concentration rate changed by wet deposition. +Build Wetdeposition model +# Example +``` julia + @parameters t + wd = Wetdeposition(t) +``` +""" + struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem function Wetdeposition(t) From 50596649fc4920a809772fe624c7982d74250181 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Wed, 19 Jul 2023 23:31:38 +0800 Subject: [PATCH 28/32] Update EarthSciMlBase versiont to 0.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e1b9f701..4ed6fd58 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" julia = "1" -EarthSciMLBase = "0.4" +EarthSciMLBase = "0.6" ModelingToolkit = "8" From 03be7e6da38ce64107abfbea8eefc6f735b67346 Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Mon, 24 Jul 2023 23:55:57 +0800 Subject: [PATCH 29/32] Small changes following suggestions in pull request --- src/dry_deposition.jl | 28 ++++++++-------------------- src/wet_deposition.jl | 13 ++++--------- 2 files changed, 12 insertions(+), 29 deletions(-) diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 1bf175ea..a225ce7f 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,17 +1,17 @@ export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle, DrydepositionG -@constants g = 9.81 [unit = u"m*s^-2"] # gravitational acceleration [m/s2] +@constants g = 9.81 [unit = u"m*s^-2", description = "gravitational acceleration"] @constants κ = 0.4 # von Karman constant -@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K"] # Boltzmann constant -@constants M_air = 28.97e-3 [unit = u"kg/mol"] # molecular weight of air -@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1"] # Gas constant +@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K", description = "Boltzmann constant"] +@constants M_air = 28.97e-3 [unit = u"kg/mol", description = "molecular weight of air"] +@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1", description = "Gas constant"] """ Function Ra calculates aerodynamic resistance to dry deposition where z is the top of the surface layer [m], z₀ is the roughness length [m], u_star is friction velocity [m/s], and L is Monin-Obukhov length [m] -Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 & 19.14. +Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd Edition, John Wiley & Sons, New York.] +equation 19.13 & 19.14. """ - @constants unit_m = 1 [unit = u"m"] function ra(z, z₀, u_star, L) ζ = IfElse.ifelse((L/unit_m == 0), 0, z/L) @@ -26,12 +26,10 @@ end """ Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. """ - @constants unit_T = 1 [unit = u"K"] @constants unit_convert_mu = 1 [unit = u"kg/m/s"] function mu(T) return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu - #return (1.8e-5*(T/T₀)^0.85)*unit_convert_mu end """ @@ -58,15 +56,10 @@ Function vs calculates the terminal setting velocity of a particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. From equation 9.42 in Seinfeld and Pandis (2006). """ -# Particle diameter Dₚ greater than 20um; Stokes settling no longer applies. @constants unit_v = 1 [unit = u"m/s"] function vs(Dₚ, ρₚ, Cc, μ) IfElse.ifelse((Dₚ > 20.e-6*unit_m), 99999999*unit_v, Dₚ^2*ρₚ*g*Cc/(18*μ)) - # if Dₚ > 20.e-6u"m" - # print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.") - # else - # return Dₚ^2*ρₚ*g*Cc/(18*μ) - # end + # Particle diameter Dₚ greater than 20um; Stokes settling no longer applies. end """ @@ -82,11 +75,9 @@ end Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K] using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf """ - @constants T_unitless = 1 [unit = u"K^-1"] @constants unit_dH2O = 1 [unit = u"m^2/s"] function dH2O(T) - #T_unitless = T*1u"K^-1" return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O end @@ -172,7 +163,6 @@ where Sc is the dimensionless Schmidt number, u_star is the friction velocity [m Dp is particle diameter [m], and iSeason and iLandUse are season and land use indexes, respectively. From Seinfeld and Pandis (2006) equation 19.27. """ - function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) α = α_table[iLandUse] γ = γ_table[iLandUse] @@ -193,13 +183,12 @@ irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. Based on Seinfeld and Pandis (2006) equation 19.2. """ - @constants G_unitless = 1 [unit = u"m^2/W"] @constants Rc_unit = 1 [unit = u"s/m"] function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) - Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s] + Dg = dH2O(Ts)/gasData.Dh2oPerDx # Diffusivity of gas of interest [m2/s] Sc = sc(μ,ρA, Dg) Rb = RbGas(Sc, u_star) Rc = SurfaceResistance(gasData, G*G_unitless, (Ts*T_unitless-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)*Rc_unit @@ -247,7 +236,6 @@ Build Drydeposition model (gas) d = DrydepositionG(t) ``` """ - struct DrydepositionG <: EarthSciMLODESystem sys::ODESystem function DrydepositionG(t) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index c4926c34..a23d6023 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -1,5 +1,9 @@ export WetDeposition, Wetdeposition, wd_defaults +@constants A_wd = 5.2 [unit = u"m^3/kg/s", deposition = "Empirical coefficient"] +@constants ρwater = 1000.0 [unit = u"kg*m^-3"] +@constants Vdr = 5.0 [unit = u"m/s", deposition = "raindrop fall speed"] + """ Calculate wet deposition based on formulas at www.emep.int/UniDoc/node12.html. @@ -9,21 +13,13 @@ and fall distance (Δz [m]). Outputs are wet deposition rates for PM2.5, SO2, and other gases (wdParticle, wdSO2, and wdOtherGas [1/s]). """ - -@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient -@constants ρwater = 1000.0 [unit = u"kg*m^-3"] # kg/m3 -@constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s - function WetDeposition(cloudFrac, qrain, ρ_air, Δz) - #@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient E = 0.1 # size-dependent collection efficiency of aerosols by the raindrops wSubSO2 = 0.15 # sub-cloud scavanging ratio wSubOther = 0.5 # sub-cloud scavanging ratio wInSO2 = 0.3 # in-cloud scavanging ratio wInParticle = 1.0 # in-cloud scavanging ratio wInOther = 1.4 # in-cloud scavanging ratio - # @constants ρwater = 1000.0 [unit = u"kg*m^-3"] # kg/m3 - # @constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s # precalculated constant combinations AE = A_wd * E @@ -61,7 +57,6 @@ Build Wetdeposition model wd = Wetdeposition(t) ``` """ - struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem function Wetdeposition(t) From f7827a7520492d71448401c4c4a8c29b3743345d Mon Sep 17 00:00:00 2001 From: Jialin Liu Date: Sun, 6 Aug 2023 00:07:55 +0800 Subject: [PATCH 30/32] Format change debug --- src/dry_deposition.jl | 18 +++++++++--------- src/wet_deposition.jl | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index a225ce7f..a53c29c8 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -6,13 +6,13 @@ export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGa @constants M_air = 28.97e-3 [unit = u"kg/mol", description = "molecular weight of air"] @constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1", description = "Gas constant"] +@constants unit_m = 1 [unit = u"m"] """ Function Ra calculates aerodynamic resistance to dry deposition where z is the top of the surface layer [m], z₀ is the roughness length [m], u_star is friction velocity [m/s], and L is Monin-Obukhov length [m] Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd Edition, John Wiley & Sons, New York.] equation 19.13 & 19.14. """ -@constants unit_m = 1 [unit = u"m"] function ra(z, z₀, u_star, L) ζ = IfElse.ifelse((L/unit_m == 0), 0, z/L) ζ₀ = IfElse.ifelse((L/unit_m == 0), 0, z₀/L) @@ -23,11 +23,11 @@ function ra(z, z₀, u_star, L) return rₐ end +@constants unit_T = 1 [unit = u"K"] +@constants unit_convert_mu = 1 [unit = u"kg/m/s"] """ Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. """ -@constants unit_T = 1 [unit = u"K"] -@constants unit_convert_mu = 1 [unit = u"kg/m/s"] function mu(T) return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu end @@ -51,12 +51,12 @@ function cc(Dₚ,T,P,μ) return 1+2*λ/Dₚ*(1.257+0.4*exp(-1.1*Dₚ/(2*λ))) end +@constants unit_v = 1 [unit = u"m/s"] """ Function vs calculates the terminal setting velocity of a particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. From equation 9.42 in Seinfeld and Pandis (2006). """ -@constants unit_v = 1 [unit = u"m/s"] function vs(Dₚ, ρₚ, Cc, μ) IfElse.ifelse((Dₚ > 20.e-6*unit_m), 99999999*unit_v, Dₚ^2*ρₚ*g*Cc/(18*μ)) # Particle diameter Dₚ greater than 20um; Stokes settling no longer applies. @@ -71,12 +71,12 @@ function dParticle(T,P,Dₚ,Cc,μ) return k*T*Cc/(3*pi*μ*Dₚ) end +@constants T_unitless = 1 [unit = u"K^-1"] +@constants unit_dH2O = 1 [unit = u"m^2/s"] """ Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K] using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf -""" -@constants T_unitless = 1 [unit = u"K^-1"] -@constants unit_dH2O = 1 [unit = u"m^2/s"] +""" function dH2O(T) return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O end @@ -174,6 +174,8 @@ function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) return 1/(3*u_star*(term_1+term_2+term_3)*R1) end +@constants G_unitless = 1 [unit = u"m^2/W"] +@constants Rc_unit = 1 [unit = u"s/m"] """ Function DryDepGas calculates dry deposition velocity [m/s] for a gas species, where z is the height of the surface layer [m], zo is roughness length [m], u_star is friction velocity [m/s], @@ -183,8 +185,6 @@ irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. Based on Seinfeld and Pandis (2006) equation 19.2. """ -@constants G_unitless = 1 [unit = u"m^2/W"] -@constants Rc_unit = 1 [unit = u"s/m"] function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index a23d6023..4cbb4a01 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -1,8 +1,8 @@ export WetDeposition, Wetdeposition, wd_defaults -@constants A_wd = 5.2 [unit = u"m^3/kg/s", deposition = "Empirical coefficient"] +@constants A_wd = 5.2 [unit = u"m^3/kg/s", description = "Empirical coefficient"] @constants ρwater = 1000.0 [unit = u"kg*m^-3"] -@constants Vdr = 5.0 [unit = u"m/s", deposition = "raindrop fall speed"] +@constants Vdr = 5.0 [unit = u"m/s", description = "raindrop fall speed"] """ Calculate wet deposition based on formulas at From fb89fa536c83b52aa4af8efcfa4cff739162de2b Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Thu, 15 Feb 2024 14:20:57 -0600 Subject: [PATCH 31/32] Add more stable gases in the wet deposition functions --- src/wet_deposition.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 408932d7..102e6a49 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -21,7 +21,7 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) wInParticle = 1.0 # in-cloud scavanging ratio wInOther = 1.4 # in-cloud scavanging ratio - # precalculated constant combinations + # precalculated constant combinations AE = A_wd * E wSubSO2VdrPerρwater = wSubSO2 * Vdr / ρwater wSubOtherVdrPerρwater = wSubOther * Vdr / ρwater @@ -29,11 +29,11 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) wInParticleVdrPerρwater = wInParticle * Vdr / ρwater wInOtherVdrPerρwater = wInOther * Vdr / ρwater - # wdParticle (subcloud) = A * P / Vdr * E; P = QRAIN * Vdr * ρgas => wdParticle = A * QRAIN * ρgas * E + # wdParticle (subcloud) = A * P / Vdr * E; P = QRAIN * Vdr * ρgas => wdParticle = A * QRAIN * ρgas * E # wdGas (subcloud) = wSub * P / Δz / ρwater = wSub * QRAIN * Vdr * ρgas / Δz / ρwater # wd (in-cloud) = wIn * P / Δz / ρwater = wIn * QRAIN * Vdr * ρgas / Δz / ρwater - wdParticle = qrain * ρ_air * (AE + cloudFrac * (wInParticleVdrPerρwater / Δz)) + wdParticle = qrain * ρ_air * (AE + cloudFrac * (wInParticleVdrPerρwater / Δz)) wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz @@ -70,12 +70,22 @@ struct Wetdeposition <: EarthSciMLODESystem @variables SO2(t) [unit = u"ppb"] @variables O3(t) [unit = u"ppb"] + @variables NO2(t) [unit = u"ppb"] + @variables CH4(t) [unit = u"ppb"] + @variables CO(t) [unit = u"ppb"] + @variables DMS(t) [unit = u"ppb"] + @variables ISOP(t) [unit = u"ppb"] eqs = [ D(SO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2 D(O3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * O3 + D(NO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * NO2 + D(CH4) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH4 + D(CO) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CO + D(DMS) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * DMS + D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP ] new(ODESystem(eqs, t, [SO2, O3], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) end -end +end \ No newline at end of file From 83a45fc018c509ce0090d4661a68a531b3ed6d53 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Thu, 15 Feb 2024 15:47:45 -0600 Subject: [PATCH 32/32] Add more variables --- src/wet_deposition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 102e6a49..1f4d2bf8 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -86,6 +86,6 @@ struct Wetdeposition <: EarthSciMLODESystem D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP ] - new(ODESystem(eqs, t, [SO2, O3], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) + new(ODESystem(eqs, t, [SO2, O3, NO2, CH4, CO, DMS, ISOP], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) end end \ No newline at end of file