From 68db407da3cc750e163bc4f35bf87849a544a4b3 Mon Sep 17 00:00:00 2001 From: Christopher Tessum Date: Sun, 18 Feb 2024 13:33:07 -0600 Subject: [PATCH] Update name and format --- .github/workflows/CI.yml | 6 +- Project.toml | 2 +- README.md | 10 +- docs/Project.toml | 2 +- docs/make.jl | 16 +- docs/src/index.md | 8 +- src/dry_deposition.jl | 84 ++++---- src/wesley1989.jl | 426 +++++++++++++++++++-------------------- src/wet_deposition.jl | 70 +++---- test/drydep_test.jl | 51 +++-- test/runtests.jl | 6 +- test/wesely1989_test.jl | 232 ++++++++++----------- test/wetdep_test.jl | 4 +- 13 files changed, 458 insertions(+), 459 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index a553f92e..ee0b1de9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -55,9 +55,9 @@ jobs: - run: | julia --project=docs -e ' using Documenter: DocMeta, doctest - using DepositionMTK - DocMeta.setdocmeta!(DepositionMTK, :DocTestSetup, :(using DepositionMTK); recursive=true) - doctest(DepositionMTK)' + using AtmosphericDeposition + DocMeta.setdocmeta!(AtmosphericDeposition, :DocTestSetup, :(using AtmosphericDeposition); recursive=true) + doctest(AtmosphericDeposition)' - run: julia --project=docs docs/make.jl env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/Project.toml b/Project.toml index b9556369..03c83841 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,4 @@ -name = "DepositionMTK" +name = "AtmosphericDeposition" uuid = "1a52f20c-0d16-41d8-a00a-b2996d86a462" authors = ["EarthSciML authors and contributors"] version = "0.1.0" diff --git a/README.md b/README.md index 1ef9c178..e7afad1a 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ -# DepositionMTK +# AtmosphericDeposition -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://earthsciml.github.io/DepositionMTK.jl/stable) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://earthsciml.github.io/DepositionMTK.jl/dev) -[![Build Status](https://github.com/EarthSciML/DepositionMTK.jl/workflows/CI/badge.svg)](https://github.com/earthsciml/DepositionMTK.jl/actions) -[![Coverage](https://codecov.io/gh/EarthSciML/DepositionMTK.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/earthsciml/DepositionMTK.jl) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://earthsciml.github.io/AtmosphericDeposition.jl/stable) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://earthsciml.github.io/AtmosphericDeposition.jl/dev) +[![Build Status](https://github.com/EarthSciML/AtmosphericDeposition.jl/workflows/CI/badge.svg)](https://github.com/earthsciml/AtmosphericDeposition.jl/actions) +[![Coverage](https://codecov.io/gh/EarthSciML/AtmosphericDeposition.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/earthsciml/AtmosphericDeposition.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) diff --git a/docs/Project.toml b/docs/Project.toml index 7f5bbc39..479f84cb 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,3 @@ [deps] -DepositionMTK = "1a52f20c-0d16-41d8-a00a-b2996d86a462" +AtmosphericDeposition = "1a52f20c-0d16-41d8-a00a-b2996d86a462" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/docs/make.jl b/docs/make.jl index 6bfa55d8..e4ba3085 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,16 +1,16 @@ -using DepositionMTK +using AtmosphericDeposition using Documenter -DocMeta.setdocmeta!(DepositionMTK, :DocTestSetup, :(using DepositionMTK); recursive=true) +DocMeta.setdocmeta!(AtmosphericDeposition, :DocTestSetup, :(using AtmosphericDeposition); recursive=true) makedocs(; - modules=[DepositionMTK], - authors="Chris Tessum and contributors", - repo="https://github.com/ctessum/DepositionMTK.jl/blob/{commit}{path}#{line}", - sitename="DepositionMTK.jl", + modules=[AtmosphericDeposition], + authors="EarthSciML authors and contributors", + repo="https://github.com/earthsciml/AtmosphericDeposition.jl/blob/{commit}{path}#{line}", + sitename="AtmosphericDeposition.jl", format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://ctessum.github.io/DepositionMTK.jl", + canonical="https://earthciml.github.io/AtmosphericDeposition.jl", assets=String[], ), pages=[ @@ -19,5 +19,5 @@ makedocs(; ) deploydocs(; - repo="github.com/ctessum/DepositionMTK.jl", + repo="github.com/earthsciml/AtmosphericDeposition.jl", ) diff --git a/docs/src/index.md b/docs/src/index.md index 23fd3d09..81bc27a4 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,10 +1,10 @@ ```@meta -CurrentModule = DepositionMTK +CurrentModule = AtmosphericDeposition ``` -# DepositionMTK +# AtmosphericDeposition -Documentation for [DepositionMTK](https://github.com/EarthSciML/DepositionMTK.jl). +Documentation for [AtmosphericDeposition](https://github.com/EarthSciML/AtmosphericDeposition.jl). ```@index ``` @@ -40,5 +40,5 @@ and variations of uptake characteristics by individual plant species within the ```@autodocs -Modules = [DepositionMTK] +Modules = [AtmosphericDeposition] ``` diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 0a750154..9b1a0945 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -3,8 +3,8 @@ export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGa @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", 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"] +@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"] """ @@ -14,12 +14,12 @@ Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmo equation 19.13 & 19.14. """ 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) - 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) + ζ = 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 @@ -29,7 +29,7 @@ end Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. """ function mu(T) - return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu + return (1.458 * 10^-6 * (T / unit_T)^(3 / 2) / ((T / unit_T) + 110.4)) * unit_convert_mu end """ @@ -37,8 +37,8 @@ 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) +function mfp(T, P, μ) + return 2 * μ / (P * (8 * M_air / (pi * R * T))^0.5) end """ @@ -46,9 +46,9 @@ 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*λ))) +function cc(Dₚ, T, P, μ) + λ = mfp(T, P, μ) + return 1 + 2 * λ / Dₚ * (1.257 + 0.4 * exp(-1.1 * Dₚ / (2 * λ))) end @constants unit_v = 1 [unit = u"m/s"] @@ -58,7 +58,7 @@ particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], C From equation 9.42 in Seinfeld and Pandis (2006). """ function vs(Dₚ, ρₚ, Cc, μ) - IfElse.ifelse((Dₚ > 20.e-6*unit_m), 99999999*unit_v, Dₚ^2*ρₚ*g*Cc/(18*μ)) + 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. end @@ -67,8 +67,8 @@ Function dParticle calculates the brownian diffusivity of a particle [m2/s] usin (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ₚ) +function dParticle(T, P, Dₚ, Cc, μ) + return k * T * Cc / (3 * pi * μ * Dₚ) end @constants T_unitless = 1 [unit = u"K^-1"] @@ -76,9 +76,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 -""" +""" function dH2O(T) - return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O + return (-2.775e-6 + 4.479e-8 * T * T_unitless + 1.656e-10 * (T * T_unitless)^2) * unit_dH2O end """ @@ -86,7 +86,7 @@ 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) + return μ / (ρ * D) end """ @@ -95,7 +95,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 """ @@ -104,7 +104,7 @@ where vs is settling velocity [m/s], u_star is friction velocity [m/s], and A is based on Seinfeld and Pandis (2006) equation 19.24. """ function stVeg(vₛ, u_star, A) - return vₛ*u_star/(g*A) + return vₛ * u_star / (g * A) end """ @@ -113,7 +113,7 @@ where Sc is the dimensionless Schmidt number and u_star is the friction velocity From Seinfeld and Pandis (2006) equation 19.17. """ function RbGas(Sc, u_star) - return 5*Sc^(2/3)/u_star + return 5 * Sc^(2 / 3) / u_star end """ @@ -163,15 +163,15 @@ 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) +function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) α = α_table[iLandUse] γ = γ_table[iLandUse] - A = (A_table[iSeason,iLandUse]*10^(-3))*unit_m + A = (A_table[iSeason, iLandUse] * 10^(-3)) * unit_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) + 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 @constants G_unitless = 1 [unit = u"m^2/W"] @@ -185,14 +185,14 @@ 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. """ -function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) +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) + 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 - return 1/(Ra+Rb+Rc) + 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 """ @@ -208,14 +208,14 @@ 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))*unit_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) Rb = RbParticle(Sc, u_star, St, Dp, iSeason, iLandUse) - return 1/(Ra+Rb+Ra*Rb*Vs)+Vs + 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] @@ -243,10 +243,10 @@ struct DrydepositionG <: EarthSciMLODESystem 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 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"] @@ -262,7 +262,7 @@ struct DrydepositionG <: EarthSciMLODESystem @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(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 @@ -272,5 +272,5 @@ struct DrydepositionG <: EarthSciMLODESystem new(ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG)) end -end +end diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 432ee806..191cfc57 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -1,213 +1,213 @@ - 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 - 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 - 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, 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. - 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. - r_c = min(r_c, 9999.0) - return r_c - end +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 + 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 + 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, 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. + 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. + r_c = min(r_c, 9999.0) + return r_c +end diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 9d568556..ba62a4ac 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", description = "Empirical coefficient"] -@constants ρwater = 1000.0 [unit = u"kg*m^-3"] -@constants Vdr = 5.0 [unit = u"m/s", description = "raindrop fall speed"] +@constants ρwater = 1000.0 [unit = u"kg*m^-3"] +@constants Vdr = 5.0 [unit = u"m/s", description = "raindrop fall speed"] """ Calculate wet deposition based on formulas at @@ -14,30 +14,30 @@ Outputs are wet deposition rates for PM2.5, SO2, and other gases (wdParticle, wdSO2, and wdOtherGas [1/s]). """ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) - 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 + 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 # precalculated constant combinations - AE = A_wd * 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 + AE = A_wd * 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 + # 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 + wdSO2 = (wSubSO2VdrPerρwater + cloudFrac * wSubSO2VdrPerρwater) * qrain * ρ_air / Δz + wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac * wSubOtherVdrPerρwater) * qrain * ρ_air / Δz - return wdParticle, wdSO2, wdOtherGas + return wdParticle, wdSO2, wdOtherGas end wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0] @@ -60,31 +60,31 @@ Build Wetdeposition model struct Wetdeposition <: EarthSciMLODESystem sys::ODESystem 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 cloudFrac = 0.5 + @parameters qrain = 0.5 + @parameters ρ_air = 1.204 [unit = u"kg*m^-3"] + @parameters Δz = 200 [unit = u"m"] D = Differential(t) @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"] + @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(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 + 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, NO2, CH4, CO, DMS, ISOP], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) end -end \ No newline at end of file +end \ No newline at end of file diff --git a/test/drydep_test.jl b/test/drydep_test.jl index 0d3ae8af..3bb53a50 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -1,5 +1,5 @@ -using DepositionMTK -using Test,Unitful, ModelingToolkit +using AtmosphericDeposition +using Test, Unitful, ModelingToolkit begin @parameters T [unit = u"K"] @@ -17,8 +17,8 @@ begin end @testset "mfp" begin - @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 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 @@ -34,13 +34,13 @@ end @testset "viscosity" begin T_ = [275, 300, 325, 350, 375, 400] - μ_list = [1.725, 1.846, 1.962, 2.075, 2.181, 2.286].*10^-5 + μ_list = [1.725, 1.846, 1.962, 2.075, 2.181, 2.286] .* 10^-5 μ_test = [] 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 + @test (μ_test[i] - μ_list[i]) / μ_list[i] < 0.01 end end @@ -49,46 +49,45 @@ end 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, substitute(cc(Dp,T,P,μ), Dict(Dp => Dp_[i]*1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, defaults...)), defaults...))) + 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 + for i in 1:16 + @test (Cc_test[i] - Cc_list[i]) / Cc_list[i] < 0.03 end end -@parameters Cc +@parameters Cc @testset "Vs" begin Dp_ = [0.01, 0.1, 1, 10] - Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 # convert to m/s + 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!(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 + 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 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 + @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 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_true = [0.5, 0.012, 0.02, 0.6] ./ 100 # [m/s] vd_list = [] for i in 1:4 - 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...))) + 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 + @test vd_list[i] - vd_true[i] < 0.015 end end - - - - \ No newline at end of file + + + diff --git a/test/runtests.jl b/test/runtests.jl index 86a240d4..c36e7f4f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,11 +1,11 @@ -using DepositionMTK +using AtmosphericDeposition using Test, SafeTestsets -@testset "DepositionMTK.jl" begin +@testset "AtmosphericDeposition.jl" begin @safetestset "drydep_test.jl" begin include("drydep_test.jl") end - + @safetestset "wesely1989_test.jl" begin include("wesely1989_test.jl") end diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl index d0639507..c3995a00 100644 --- a/test/wesely1989_test.jl +++ b/test/wesely1989_test.jl @@ -1,149 +1,149 @@ -using DepositionMTK -using Test,StaticArrays +using AtmosphericDeposition +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 - 1100 1100 1100 1100 1200 90 1000 - 1000 1000 1000 1000 1100 1100 1100 - 270 290 330 620 1100 90 1000] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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] + 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 + 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 + 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 diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index f60b465c..8157aeb6 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -1,5 +1,5 @@ -using DepositionMTK -using Test,Unitful, ModelingToolkit +using AtmosphericDeposition +using Test, Unitful, ModelingToolkit @parameters cloudFrac = 0.5 @parameters qrain = 0.5