Update documentation andinclude species HNO3, H2O2, CH2O

docs/Project.toml

@@ -4,5 +4,5 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/src/dry_deposition.md

@@ -3,24 +3,24 @@
 This is an implementation of a box model used to calculate changes in gas species concentration due to dry deposition.
 
 ## Running the model
-Here's an example of how concentration of different species, such as SO₂, O₃, NO₂, NO, H₂O₂ and CH₂O change due to dry deposition. 
+Here's an example of how concentration of different species, such as SO₂, O₃, NO₂, NO, H₂O₂, CH₂O and HNO₃change due to dry deposition. 
 
 We can create an instance of the model in the following manner:
 ```julia
 using AtmosphericDeposition
 using ModelingToolkit
 using DifferentialEquations
 using EarthSciMLBase
-using Unitful
+using DynamicQuantities
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description="Time"]
-model = DrydepositionG(t)
+model = DrydepositionG()
 ```
 Before running any simulations with the model we need to convert it into a system of differential equations.
 ```julia
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24) 
-u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O
+u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃
 sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```
 
@@ -29,14 +29,14 @@ using AtmosphericDeposition
 using ModelingToolkit
 using DifferentialEquations
 using EarthSciMLBase
-using Unitful
+using DynamicQuantities
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description="Time"]
-model = DrydepositionG(t)
+model = DrydepositionG()
 
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24) 
-u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O
+u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃
 sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```
 
@@ -49,16 +49,16 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig
 ## Parameters
 The parameters in the model are:
 ```julia
-parameters(sys) # [z, z₀, u_star, L, ρA, G, T, θ]
+parameters(sys) # [iLandUse, z, z₀, u_star, L, ρA, G, iSeason, T, θ]
 ```
-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], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians].
+where ```iSeason``` and ```iLandUse``` are the index for the season and landuse,```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], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians].
 
 Let's run some simulation with different value for parameter ```z```. 
 ```@example 1
 @unpack O3 = sys
 
-p1 = [50,0.04,0.44,0,1.2,300,298,0]
-p2 = [10,0.04,0.44,0,1.2,300,298,0]
+p1 = [10, 50,0.04,0.44,0,1.2,300,1,298,0]
+p2 = [10, 10,0.04,0.44,0,1.2,300,1,298,0]
 sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
 sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 

docs/src/emep.md

@@ -8,17 +8,17 @@ using AtmosphericDeposition
 using ModelingToolkit
 using DifferentialEquations
 using EarthSciMLBase
-using Unitful
+using DynamicQuantities
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description="Time"]
-model = Wetdeposition(t)
+model = Wetdeposition()
 ```
 
 Before running any simulations with the model we need to convert it into a system of differential equations.
 ```julia 
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24)
-u0 = [2.0,10.0,5,1400,275,50,0.15]  # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP
+u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15]  # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O
 prob = ODEProblem(sys, u0, tspan, [])
 sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```
@@ -28,14 +28,14 @@ using AtmosphericDeposition
 using ModelingToolkit
 using DifferentialEquations
 using EarthSciMLBase
-using Unitful
+using DynamicQuantities
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description="Time"]
-model = Wetdeposition(t)
+model = Wetdeposition()
 
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24)
-u0 = [2.0,10.0,5,1400,275,50,0.15]  # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP
+u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15]  # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O
 prob = ODEProblem(sys, u0, tspan, [])
 sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```

ext/EarthSciDataExt.jl

@@ -13,7 +13,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler,
         gg = 9.81, [unit = u"m*s^-2", description="gravitational acceleration"],
     )
 
-    # ρA in DrydepositionG(t) are in units of "kg/m3".
+    # ρA in DrydepositionG() are in units of "kg/m3".
     # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3
     # Overall, ρA = P*M/(R*T)*kgperg
 
@@ -40,7 +40,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g
         Vdr = 5.0, [unit = u"m/s", description="droplet velocity"],
     )
 
-    # ρ_air in Wetdeposition(t) are in units of "kg/m3".
+    # ρ_air in Wetdeposition() are in units of "kg/m3".
     # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3
     # Overall, ρ_air = P*M/(R*T)*kgperg
 

ext/GasChemExt.jl

@@ -26,6 +26,8 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos
         c.CO => d.CO,
         c.DMS => d.DMS,
         c.ISOP => d.ISOP,
+        c.H2O2 => d.H2O2,
+        c.CH2O => d.CH2O,
     ))
 end
 

src/dry_deposition.jl

@@ -260,6 +260,7 @@ function DrydepositionG(; name=:DrydepositionG)
         NO(t), [unit = u"ppb"],
         H2O2(t), [unit = u"ppb"],
         CH2O(t), [unit = u"ppb"],
+        HNO3(t), [unit = u"ppb"],
     )
 
     eqs = [
@@ -269,6 +270,7 @@ function DrydepositionG(; name=:DrydepositionG)
         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
+        D(HNO3) ~ -DryDepGas(z, z₀, u_star, L, ρA, Hno3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * HNO3
     ]
 
     ODESystem(eqs, t, vars, params; name=name,

src/wet_deposition.jl

@@ -72,6 +72,9 @@ function Wetdeposition(; name=:Wetdeposition)
         CO(t), [unit = u"ppb"],
         DMS(t), [unit = u"ppb"],
         ISOP(t), [unit = u"ppb"],
+        H2O2(t), [unit = u"ppb"],
+        HNO3(t), [unit = u"ppb"],
+        CH2O(t), [unit = u"ppb"],
     )
 
     eqs = [
@@ -82,6 +85,9 @@ function Wetdeposition(; name=:Wetdeposition)
         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(H2O2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * H2O2
+        D(HNO3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * HNO3
+        D(CH2O) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH2O
     ]
 
     ODESystem(eqs, t, vars, params; name=name,

test/connector_test.jl

@@ -7,7 +7,9 @@ using ModelingToolkit:t
     composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition())
     combined_mtk = convert(ODESystem, composed_ode)
     sys = structural_simplify(combined_mtk)
-    @test length(unknowns(sys)) ≈ 18
+    print(unknowns(sys))
+    @test length(unknowns(sys)) ≈ 20
+    #TODO Change 20 to 18 after the latest GasChem package include species HNO3
 
     eqs = string(equations(sys))
     wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"]
@@ -23,7 +25,8 @@ end
     model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG())
 
     sys = structural_simplify(convert(ODESystem, model))
-    @test length(unknowns(sys)) ≈ 18
+    @test length(unknowns(sys)) ≈ 20
+    #TODO Change 20 to 18 after the latest GasChem package include species HNO3
 
     eqs = string(equations(convert(ODESystem, model)))
     wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)"