From 3a25667e0427d9799a3629aab2d8d47dd404440e Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Tue, 20 Aug 2024 22:24:31 -0700
Subject: [PATCH 01/11] Update with EarthSciMLBasev0.15

---
 Project.toml           |  8 ++++--
 ext/EarthSciDataExt.jl | 51 ++++++++++++++++++++++++++++++++++
 ext/GasChemExt.jl      | 63 +++++++++++++++---------------------------
 src/dry_deposition.jl  | 44 +++++++++++++++--------------
 src/wet_deposition.jl  | 38 +++++++++++++------------
 test/connector_test.jl | 31 +++++++++++++++++----
 6 files changed, 149 insertions(+), 86 deletions(-)
 create mode 100644 ext/EarthSciDataExt.jl

diff --git a/Project.toml b/Project.toml
index 2e2d7eec..f70640de 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,18 +11,21 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
+EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
 GasChem = "58070593-4751-4c87-a5d1-63807d11d76c"
 
 [extensions]
 GasChemExt = "GasChem"
+EarthSciDataExt = "EarthSciData"
 
 [compat]
 DifferentialEquations = "7"
-EarthSciMLBase = "0.10"
+EarthSciMLBase = "0.15"
 GasChem = "0.6"
 IfElse = "0.1"
 ModelingToolkit = "8"
@@ -33,9 +36,10 @@ julia = "1"
 
 [extras]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
 GasChem = "58070593-4751-4c87-a5d1-63807d11d76c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "GasChem", "Dates"]
+test = ["Test", "SafeTestsets", "GasChem", "Dates", "EarthSciData"]
diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl
new file mode 100644
index 00000000..08babae6
--- /dev/null
+++ b/ext/EarthSciDataExt.jl
@@ -0,0 +1,51 @@
+module EarthSciDataExt
+
+using AtmosphericDeposition, EarthSciData, EarthSciMLBase, Unitful, ModelingToolkit
+
+function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, g::EarthSciData.GEOSFPCoupler)
+    d, g = d.sys, g.sys
+
+    @constants(
+        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
+        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
+        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+    )
+
+    # ρA in DrydepositionG(t) 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
+
+    d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA)
+    ConnectorSystem([
+            d.T ~ g.I3₊T,
+            d.z ~ g.A1₊PBLH,
+            d.z₀ ~ g.A1₊Z0M,
+            d.u_star ~ g.A1₊USTAR,
+            d.G ~ g. A1₊SWGDN,
+            d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+        ], d, g)
+end
+
+function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g::EarthSciData.GEOSFPCoupler)
+    d, g = d.sys, g.sys
+
+    @constants(
+        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
+        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
+        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+    )
+
+    # ρ_air in Wetdeposition(t) 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
+
+    d = param_to_var(d, :cloudFrac, :ρ_air)
+    ConnectorSystem([
+            d.cloudFrac ~ g.A3cld₊CLOUD,
+            d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+        ], d, g)
+end
+
+end
\ No newline at end of file
diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl
index 00c3b792..f7420fca 100644
--- a/ext/GasChemExt.jl
+++ b/ext/GasChemExt.jl
@@ -2,48 +2,31 @@ module GasChemExt
 
 using AtmosphericDeposition, GasChem, EarthSciMLBase, Unitful, ModelingToolkit
 
-export register_couplings_ext
-
-# Use a global flag to track initialization and ensure that the coupling between SuperFast and NEI Emission is only initialized once
-const couplings_registered_ext = Ref(false)
-
-function register_couplings_ext()
-    if couplings_registered_ext[]
-        println("Couplings have already been registered.")
-        return
-    end
-    
-    println("Registering couplings in ext")
-    @parameters t [unit = u"s"]
-
-    register_coupling(SuperFast(t), Wetdeposition(t)) do c, e
-        operator_compose(convert(ODESystem, c), e, Dict(
-            c.SO2 => e.SO2,
-            c.NO2 => e.NO2,
-            c.O3 => e.O3,
-            c.CH4 => e.CH4,
-            c.CO => e.CO,
-            c.DMS => e.DMS,
-            c.ISOP => e.ISOP))
-    end
-
-    register_coupling(SuperFast(t), DrydepositionG(t)) do c, e
-        operator_compose(convert(ODESystem, c), e, Dict(
-            c.SO2 => e.SO2,
-            c.NO2 => e.NO2,
-            c.O3 => e.O3,
-            c.NO => e.NO,
-            c.H2O2 => e.H2O2,
-            c.CH2O => e.CH2O))
-    end
-    
-    couplings_registered_ext[] = true
-    println("Coupling registry after registration: ", EarthSciMLBase.coupling_registry)
+function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.DrydepositionGCoupler)
+    c, d = c.sys, d.sys
+
+    operator_compose(convert(ODESystem, c), d, Dict(
+        c.SO2 => d.SO2,
+        c.NO2 => d.NO2,
+        c.O3 => d.O3,
+        c.NO => d.NO,
+        c.H2O2 => d.H2O2,
+        c.CH2O => d.CH2O,
+    ))
 end
 
-function __init__()
-    println("Initializing EmissionExt module")
-    register_couplings_ext()
+function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.WetdepositionCoupler)
+    c, d = c.sys, d.sys
+
+    operator_compose(convert(ODESystem, c), d, Dict(
+        c.SO2 => d.SO2,
+        c.NO2 => d.NO2,
+        c.O3 => d.O3,
+        c.CH4 => d.CH4,
+        c.CO => d.CO,
+        c.DMS => d.DMS,
+        c.ISOP => d.ISOP,
+    ))
 end
 
 end
\ No newline at end of file
diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
index 9698f332..b54789ae 100644
--- a/src/dry_deposition.jl
+++ b/src/dry_deposition.jl
@@ -220,12 +220,9 @@ 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]
 
-# Add unit "ppb" to Unitful 
-module MyUnits
-using Unitful
-@unit ppb "ppb" Number 1 / 1000000000 false
+struct DrydepositionGCoupler
+    sys
 end
-Unitful.register(MyUnits)
 
 """
 Description: This is a box model used to calculate the gas species concentration rate changed by dry deposition.
@@ -236,28 +233,32 @@ Build Drydeposition model (gas)
 	d = DrydepositionG(t)
 ```
 """
-function DrydepositionG(t)
+function DrydepositionG(t; name=:DrydepositionG)
     iSeason = 1
     iLandUse = 10
     rain = false
     dew = false
-    @parameters z = 50 [unit = u"m", description = "top of the surface layer"]
-    @parameters z₀ = 0.04 [unit = u"m", description = "roughness lenght"]
-    @parameters u_star = 0.44 [unit = u"m/s", description = "friction velocity"]
-    @parameters L = 0 [unit = u"m", description = "Monin-Obukhov length"]
-    @parameters ρA = 1.2 [unit = u"kg*m^-3", description = "air density"]
-    @parameters G = 300 [unit = u"W*m^-2", description = "solar irradiation"]
-    @parameters T = 298 [unit = u"K", description = "temperature"]
-    @parameters θ = 0 [description = "slope of the local terrain, in unit radians"]
+    params = @parameters(
+        z = 50, [unit = u"m", description = "top of the surface layer"],
+        z₀ = 0.04, [unit = u"m", description = "roughness lenght"],
+        u_star = 0.44, [unit = u"m/s", description = "friction velocity"],
+        L = 0, [unit = u"m", description = "Monin-Obukhov length"],
+        ρA = 1.2, [unit = u"kg*m^-3", description = "air density"],
+        G = 300, [unit = u"W*m^-2", description = "solar irradiation"],
+        T = 298, [unit = u"K", description = "temperature"],
+        θ = 0, [description = "slope of the local terrain, in unit radians"],
+    )
 
     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"]
+    vars = @variables( 
+        SO2(t), [unit = u"nmol/mol"],
+        O3(t),[unit = u"nmol/mol"],
+        NO2(t), [unit = u"nmol/mol"],
+        NO(t), [unit = u"nmol/mol"],
+        H2O2(t), [unit = u"nmol/mol"],
+        CH2O(t), [unit = u"nmol/mol"],
+    )
 
     eqs = [
         D(SO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, iSeason, iLandUse, rain, dew, true, false) / z * SO2
@@ -268,6 +269,7 @@ function DrydepositionG(t)
         D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O
     ]
 
-    ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG)
+    ODESystem(eqs, t, vars, params; name=name,
+        metadata=Dict(:coupletype => DrydepositionGCoupler))
 end
 
diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl
index 2efd507e..d143195c 100644
--- a/src/wet_deposition.jl
+++ b/src/wet_deposition.jl
@@ -41,12 +41,9 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz)
 end
 wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0]
 
-# Add unit "ppb" to Unitful 
-module myUnits
-using Unitful
-@unit ppb "ppb" Number 1 / 1000000000 false
+struct WetdepositionCoupler
+    sys
 end
-Unitful.register(myUnits)
 
 """
 Description: This is a box model used to calculate wet deposition based on formulas at EMEP model.
@@ -57,21 +54,25 @@ Build Wetdeposition model
 	wd = Wetdeposition(t)
 ```
 """
-function Wetdeposition(t)
-    @parameters cloudFrac = 0.5 [description = "fraction of grid cell covered by clouds"]
-    @parameters qrain = 0.5 [description = "rain mixing ratio"]
-    @parameters ρ_air = 1.204 [unit = u"kg*m^-3", description = "air density"]
-    @parameters Δz = 200 [unit = u"m", description = "fall distance"]
+function Wetdeposition(t; name=:Wetdeposition)
+    params = @parameters(
+        cloudFrac = 0.5, [description = "fraction of grid cell covered by clouds"],
+        qrain = 0.5, [description = "rain mixing ratio"],
+        ρ_air = 1.204, [unit = u"kg*m^-3", description = "air density"],
+        Δz = 200, [unit = u"m", description = "fall distance"],
+    )
 
     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"]
+    vars = @variables(
+        SO2(t), [unit = u"nmol/mol"],
+        O3(t), [unit = u"nmol/mol"],
+        NO2(t), [unit = u"nmol/mol"],
+        CH4(t), [unit = u"nmol/mol"],
+        CO(t), [unit = u"nmol/mol"],
+        DMS(t), [unit = u"nmol/mol"],
+        ISOP(t), [unit = u"nmol/mol"],
+    )
 
     eqs = [
         D(SO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2
@@ -83,5 +84,6 @@ function Wetdeposition(t)
         D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP
     ]
 
-    ODESystem(eqs, t, [SO2, O3, NO2, CH4, CO, DMS, ISOP], [cloudFrac, qrain, ρ_air, Δz]; name=:WetDeposition)
+    ODESystem(eqs, t, vars, params; name=name,
+        metadata=Dict(:coupletype => WetdepositionCoupler))
 end
diff --git a/test/connector_test.jl b/test/connector_test.jl
index 5c813e74..7ee7f918 100644
--- a/test/connector_test.jl
+++ b/test/connector_test.jl
@@ -1,16 +1,37 @@
 using AtmosphericDeposition
-using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase
+using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData
 
-@testset "Connector" begin
-    ModelingToolkit.check_units(eqs...) = nothing
+@testset "GasChemExt" begin
     start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1))
-    @parameters t
+    @parameters t [unit = u"s"]
     composed_ode = couple(SuperFast(t), FastJX(t), DrydepositionG(t), Wetdeposition(t))
     tspan = (start, start+3600*24*3)
     sys = structural_simplify(get_mtk(composed_ode))
     @test length(states(sys)) ≈ 18
 
     eqs = string(equations(sys))
-    wanteqs = ["Differential(t)(superfast₊O3(t)) ~ superfast₊DrydepositionG_ddt_O3ˍt(t) + superfast₊WetDeposition_ddt_O3ˍt(t)"]
+    wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"]
     @test contains(string(eqs), wanteqs[1])
 end
+
+@testset "EarthSciDataExt" begin
+    @parameters t [unit = u"s"]
+
+    @parameters lat = 40
+    @parameters lon = -97
+    @parameters lev = 1
+    geosfp = GEOSFP("4x5", t)
+
+    model = couple(SuperFast(t), FastJX(t), geosfp, Wetdeposition(t), DrydepositionG(t))
+
+    sys = structural_simplify(get_mtk(model))
+    @test length(states(sys)) ≈ 18
+
+    eqs = string(equations(get_mtk(model)))
+    wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)"
+    @test contains(eqs, wanteq)
+    wanteq = "Wetdeposition₊cloudFrac(t) ~ GEOSFP₊A3cld₊CLOUD(t)"
+    @test contains(eqs, wanteq)
+    wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P * PaPerhPa/(GEOSFP₊I3₊T*R)*kgperg*MW_air"
+    @test contains(eqs, wanteq)
+end
\ No newline at end of file

From 46379366642c9cc22c12fee02c1d31ff4ae7e8b1 Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Fri, 6 Sep 2024 18:05:45 -0500
Subject: [PATCH 02/11] Chang iSeason & iLandUse as parameters and add
 connections with geosfp to get qrain and L

---
 Project.toml                 |  3 ++-
 docs/src/wesley1989.md       | 18 ++++++++-----
 src/AtmosphericDeposition.jl |  1 +
 src/dry_deposition.jl        | 12 ++++-----
 src/wesley1989.jl            | 51 ++++++++++++++++++++----------------
 test/drydep_test.jl          |  6 ++---
 6 files changed, 53 insertions(+), 38 deletions(-)

diff --git a/Project.toml b/Project.toml
index f70640de..977787e1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["EarthSciML authors and contributors"]
 version = "0.2.1"
 
 [deps]
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -20,8 +21,8 @@ EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
 GasChem = "58070593-4751-4c87-a5d1-63807d11d76c"
 
 [extensions]
-GasChemExt = "GasChem"
 EarthSciDataExt = "EarthSciData"
+GasChemExt = "GasChem"
 
 [compat]
 DifferentialEquations = "7"
diff --git a/docs/src/wesley1989.md b/docs/src/wesley1989.md
index 95ffc6bf..c5404633 100644
--- a/docs/src/wesley1989.md
+++ b/docs/src/wesley1989.md
@@ -39,13 +39,19 @@ a. For season index ```iSeason```, there're five seasonal categories
     4. Winter, snow on ground
     5. Transitional 
 
-b. For land use index ```iLandUse```, there're five land use categories
+b. For land use index ```iLandUse```, there're eleven land use categories
 
-    1. Evergreen–needleleaf trees
-    2. Deciduous broadleaf trees
-    3. Grass
-    4. Desert
-    5. Shrubs and interrupted woodlands 
+    1. Urban land 
+    2. agricultural land
+    3. range land
+    4. deciduous forest
+    5. coniferous forest 
+    6. mixed forest including wetland 
+    7. water, both salt and fresh 
+    8. barren land, mostly desert
+    9. nonforested wetland
+    10. mixed agricultural and range land
+    11. rocky open areas with low-growing shrubs
 
 c. For gasData ```gasData```, we include the following species: 
 
diff --git a/src/AtmosphericDeposition.jl b/src/AtmosphericDeposition.jl
index 8414bf38..a7ccfe4f 100644
--- a/src/AtmosphericDeposition.jl
+++ b/src/AtmosphericDeposition.jl
@@ -5,6 +5,7 @@ using StaticArrays
 using ModelingToolkit
 using EarthSciMLBase
 using IfElse
+using DataInterpolations
 
 include("wesley1989.jl")
 include("dry_deposition.jl")
diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
index b54789ae..3d1234ff 100644
--- a/src/dry_deposition.jl
+++ b/src/dry_deposition.jl
@@ -185,13 +185,13 @@ 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, iLandUse, 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 = WesleySurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit
+    Rc = WesleySurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit
     return 1 / (Ra + Rb + Rc)
 end
 
@@ -202,7 +202,7 @@ Dp is particle diameter [m], Ts is surface air temperature [K], P is pressure [P
 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)
+function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason, iLandUse)
     Ra = ra(z, z₀, u_star, L)
     μ = mu(Ts)
     Cc = cc(Dp, Ts, P, μ)
@@ -210,7 +210,7 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason:
     if iLandUse == 4 # dessert
         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, (obtain_value(iSeason, iLandUse,A_table) * 10^(-3)) * unit_m)
     end
     D = dParticle(Ts, P, Dp, Cc, μ)
     Sc = sc(μ, ρA, D)
@@ -234,11 +234,11 @@ Build Drydeposition model (gas)
 ```
 """
 function DrydepositionG(t; name=:DrydepositionG)
-    iSeason = 1
-    iLandUse = 10
     rain = false
     dew = false
     params = @parameters(
+        iSeason = 1, [description = "Index for season"],
+        iLandUse = 10, [description = "Index for land-use"],
         z = 50, [unit = u"m", description = "top of the surface layer"],
         z₀ = 0.04, [unit = u"m", description = "roughness lenght"],
         u_star = 0.44, [unit = u"m/s", description = "friction velocity"],
diff --git a/src/wesley1989.jl b/src/wesley1989.jl
index 62f64e9e..11fc7468 100644
--- a/src/wesley1989.jl
+++ b/src/wesley1989.jl
@@ -84,10 +84,17 @@ 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
 
+# Obtain values from matrix using symbolic parameter iSeason and iLandUse
+function obtain_value(iSeason, iLandUse, matrix)
+    index=(iLandUse-1)*5+iSeason
+    interpolate_r_i = DataInterpolations.LinearInterpolation(vec(matrix),1:55)
+    interpolate_r_i(index)
+end
+
 # 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)
+function r_s(G, Ts, iSeason, iLandUse, 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))))
+    rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts))))
     # Adjust for dew and rain (from "Effects of dew and rain" section).
     if rainOrDew
         rs *= 3
@@ -111,7 +118,7 @@ function r_smx(r_s, Dh2oPerDx::T, r_mx::T) where {T<:AbstractFloat}
 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}
+function r_lux(Hstar::T, fo::T, iSeason, iLandUse, 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
@@ -122,10 +129,10 @@ function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bo
             end
         elseif isO3
             # equation 11
-            rlux = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse]))
+            rlux = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu)))
         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))
+            rluO = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) #equation 11
+            rlux = 1.0 / (1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu) / (1.0e-5 * Hstar + fo))
                           + 1.0e-7 * Hstar + fo / rluO) # equation 14, modified to match Walmsley eq. 5g
         end
     elseif rain && (iSeason != 4)
@@ -134,31 +141,31 @@ function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bo
                 rlux = 50 #equation 13 and a half
             else
                 # equation 12
-                rlux = 1.0 / (1.0 / 5000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse]))
+                rlux = 1.0 / (1.0 / 5000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu)))
             end
         elseif isO3
             # equation 13
-            rlux = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse]))
+            rlux = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu)))
         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
+            rluO = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) # equation 13
+            rlux = 1.0 / (1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu) / (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)
+        rlux = obtain_value(iSeason, iLandUse,r_lu) / (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])
+function r_clx(Hstar::T, fo::T, iSeason, iLandUse) where {T<:AbstractFloat}
+    return 1.0 / (Hstar / (1.0e5 * obtain_value(iSeason, iLandUse,r_clS)) +
+                  fo / obtain_value(iSeason, iLandUse,r_clO))
 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])
+function r_gsx(Hstar::T, fo::T, iSeason, iLandUse) where {T<:AbstractFloat}
+    return 1.0 / (Hstar / (1.0e5 * obtain_value(iSeason, iLandUse,r_gsS)) +
+                  fo / obtain_value(iSeason, iLandUse,r_gsO))
 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).
@@ -174,17 +181,17 @@ function WesleySurfaceResistance(gasData::GasData, G, Ts, θ,
     rclx = 0.0
     rgsx = 0.0
     if isSO2
-        rclx = r_clS[iSeason, iLandUse]
-        rgsx = r_gsS[iSeason, iLandUse]
+        rclx = obtain_value(iSeason, iLandUse,r_clS)
+        rgsx = obtain_value(iSeason, iLandUse,r_gsS)
     elseif isO3
-        rclx = r_clO[iSeason, iLandUse]
-        rgsx = r_gsO[iSeason, iLandUse]
+        rclx = obtain_value(iSeason, iLandUse,r_clO)
+        rgsx = obtain_value(iSeason, iLandUse,r_gsO)
     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]
+    rac = obtain_value(iSeason, iLandUse,r_ac)
 
     # 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
diff --git a/test/drydep_test.jl b/test/drydep_test.jl
index 517e3df4..935d319e 100644
--- a/test/drydep_test.jl
+++ b/test/drydep_test.jl
@@ -13,7 +13,7 @@ begin
     @parameters ρParticle [unit = u"kg*m^-3"]
     @parameters ρA [unit = u"kg*m^-3"]
     @parameters G [unit = u"W*m^-2"]
-    @parameters θ
+    @parameters θ, iLandUse, iSeason
 end
 
 @testset "mfp" begin
@@ -28,7 +28,7 @@ end
     #@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, AtmosphericDeposition.So2Data, G, T, θ, 1, 1, false, false, true, false)) == u"m/s"
+    @test ModelingToolkit.get_unit(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.So2Data, G, T, θ, iSeason, iLandUse, 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
 
@@ -73,7 +73,7 @@ end
 
 @testset "DryDepGas" begin
     vd_true = 0.03 # m/s
-    @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.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, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, defaults...)) - vd_true) / vd_true < 0.33
 end
 
 @testset "DryDepParticle" begin

From 69aa7a91947d6a8e42b074489a36091a67674dff Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Fri, 6 Sep 2024 18:06:46 -0500
Subject: [PATCH 03/11] Add connections with geosfp to calculate L and qrain

---
 ext/EarthSciDataExt.jl | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl
index 08babae6..40632acc 100644
--- a/ext/EarthSciDataExt.jl
+++ b/ext/EarthSciDataExt.jl
@@ -10,13 +10,18 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler,
         MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
         kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+        vK = 0.4, [description = "von Karman's constant"],
+        Cp = 1000, [unit = u"W*s/kg/K", description="specific heat at constant pressure"],
+        gg = 9.81, [unit = u"m*s^-2", description="gravitational acceleration"],
     )
 
     # ρA in DrydepositionG(t) 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
 
-    d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA)
+    # Monin-Obhukov length = -Air density * Cp * T(surface air) * Ustar^3/（vK   * g  * Sensible Heat flux）
+
+    d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA, :L)
     ConnectorSystem([
             d.T ~ g.I3₊T,
             d.z ~ g.A1₊PBLH,
@@ -24,6 +29,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler,
             d.u_star ~ g.A1₊USTAR,
             d.G ~ g. A1₊SWGDN,
             d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+            d.L ~ -g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX),
         ], d, g)
 end
 
@@ -35,16 +41,21 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g
         MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
         kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+        Vdr = 5.0, [unit = u"m/s", description="droplet velocity"],
     )
 
     # ρ_air in Wetdeposition(t) 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
 
-    d = param_to_var(d, :cloudFrac, :ρ_air)
+    # From EMEP algorithm: P = QRAIN * Vdr * ρgas => QRAIN = P / Vdr / ρgas
+    # kg*m-2*s-1/(m/s)/(kg/m3)
+
+    d = param_to_var(d, :cloudFrac, :ρ_air, :qrain)
     ConnectorSystem([
             d.cloudFrac ~ g.A3cld₊CLOUD,
             d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+            d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air),
         ], d, g)
 end
 

From e8d87bfcb2ea5812ae05b6a4d37a5d03ce06fb3c Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Thu, 26 Sep 2024 15:54:06 -0500
Subject: [PATCH 04/11] Update Deposition to be compatible with MTKv9

---
 Project.toml                 | 13 +++++--------
 ext/EarthSciDataExt.jl       | 20 ++++++++------------
 ext/GasChemExt.jl            |  2 +-
 src/AtmosphericDeposition.jl |  6 ++++--
 src/dry_deposition.jl        | 26 ++++++++++++++------------
 src/wesley1989.jl            |  5 +++--
 src/wet_deposition.jl        | 16 ++++++++--------
 test/connector_test.jl       | 32 +++++++++++++++-----------------
 test/drydep_test.jl          | 22 +++++++++-------------
 test/wetdep_test.jl          |  2 +-
 10 files changed, 68 insertions(+), 76 deletions(-)

diff --git a/Project.toml b/Project.toml
index 977787e1..bcf6e927 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,15 +6,13 @@ version = "0.2.1"
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 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"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
 EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
@@ -26,13 +24,12 @@ GasChemExt = "GasChem"
 
 [compat]
 DifferentialEquations = "7"
-EarthSciMLBase = "0.15"
-GasChem = "0.6"
-IfElse = "0.1"
-ModelingToolkit = "8"
+EarthSciData = "0.9"
+EarthSciMLBase = "0.16"
+GasChem = "0.7"
+ModelingToolkit = "9"
 SafeTestsets = "0.1"
 StaticArrays = "1"
-Unitful = "1"
 julia = "1"
 
 [extras]
diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl
index 40632acc..e7bc628c 100644
--- a/ext/EarthSciDataExt.jl
+++ b/ext/EarthSciDataExt.jl
@@ -1,14 +1,12 @@
 module EarthSciDataExt
 
-using AtmosphericDeposition, EarthSciData, EarthSciMLBase, Unitful, ModelingToolkit
+using AtmosphericDeposition, EarthSciData, EarthSciMLBase, DynamicQuantities, ModelingToolkit
 
 function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, g::EarthSciData.GEOSFPCoupler)
     d, g = d.sys, g.sys
 
     @constants(
-        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
-        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
-        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
         vK = 0.4, [description = "von Karman's constant"],
         Cp = 1000, [unit = u"W*s/kg/K", description="specific heat at constant pressure"],
@@ -27,9 +25,9 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler,
             d.z ~ g.A1₊PBLH,
             d.z₀ ~ g.A1₊Z0M,
             d.u_star ~ g.A1₊USTAR,
-            d.G ~ g. A1₊SWGDN,
-            d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
-            d.L ~ -g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX),
+            d.G ~ g.A1₊SWGDN,
+            d.ρA ~ g.P/(g.I3₊T*R)*MW_air,
+            d.L ~ -g.P/(g.I3₊T*R)*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX),
         ], d, g)
 end
 
@@ -37,9 +35,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g
     d, g = d.sys, g.sys
 
     @constants(
-        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
-        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
-        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
         Vdr = 5.0, [unit = u"m/s", description="droplet velocity"],
     )
@@ -54,8 +50,8 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g
     d = param_to_var(d, :cloudFrac, :ρ_air, :qrain)
     ConnectorSystem([
             d.cloudFrac ~ g.A3cld₊CLOUD,
-            d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
-            d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air),
+            d.ρ_air ~ g.P/(g.I3₊T*R)*MW_air,
+            d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P/(g.I3₊T*R)*MW_air),
         ], d, g)
 end
 
diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl
index f7420fca..be976bf6 100644
--- a/ext/GasChemExt.jl
+++ b/ext/GasChemExt.jl
@@ -1,6 +1,6 @@
 module GasChemExt
 
-using AtmosphericDeposition, GasChem, EarthSciMLBase, Unitful, ModelingToolkit
+using AtmosphericDeposition, GasChem, EarthSciMLBase, DynamicQuantities, ModelingToolkit
 
 function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.DrydepositionGCoupler)
     c, d = c.sys, d.sys
diff --git a/src/AtmosphericDeposition.jl b/src/AtmosphericDeposition.jl
index a7ccfe4f..98100af7 100644
--- a/src/AtmosphericDeposition.jl
+++ b/src/AtmosphericDeposition.jl
@@ -1,11 +1,13 @@
 module AtmosphericDeposition
 
-using Unitful
+using DynamicQuantities
 using StaticArrays
 using ModelingToolkit
 using EarthSciMLBase
-using IfElse
 using DataInterpolations
+using ModelingToolkit:t
+
+@register_unit ppb 1
 
 include("wesley1989.jl")
 include("dry_deposition.jl")
diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
index 3d1234ff..1f18c4ad 100644
--- a/src/dry_deposition.jl
+++ b/src/dry_deposition.jl
@@ -14,12 +14,14 @@ 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₀))
+    ζ = ifelse((L / unit_m == 0), 0, z / L)
+    ζ₀ = ifelse((L / unit_m == 0), 0, z₀ / L)
+    rₐ_1 = ifelse((0 < ζ), 1 / (κ * u_star) * (log(z / z₀) + 4.7 * (ζ - ζ₀)), 1 / (κ * u_star) * log(z / z₀))
+    rₐ_1 = ifelse((ζ < 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)
+    rₐ = ifelse((-1 < ζ), 1 / (κ * u_star) * [log(z / z₀) + log(((η₀^2 + 1) * (η₀ + 1)^2) / ((η^2 + 1) * (η + 1)^2)) + 2 * (atan(η) - atan(η₀))][1], rₐ_1)
+    rₐ = ifelse((ζ < 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
 
@@ -58,7 +60,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((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
 
@@ -233,7 +235,7 @@ Build Drydeposition model (gas)
 	d = DrydepositionG(t)
 ```
 """
-function DrydepositionG(t; name=:DrydepositionG)
+function DrydepositionG(; name=:DrydepositionG)
     rain = false
     dew = false
     params = @parameters(
@@ -252,12 +254,12 @@ function DrydepositionG(t; name=:DrydepositionG)
     D = Differential(t)
 
     vars = @variables( 
-        SO2(t), [unit = u"nmol/mol"],
-        O3(t),[unit = u"nmol/mol"],
-        NO2(t), [unit = u"nmol/mol"],
-        NO(t), [unit = u"nmol/mol"],
-        H2O2(t), [unit = u"nmol/mol"],
-        CH2O(t), [unit = u"nmol/mol"],
+        SO2(t), [unit = u"ppb"],
+        O3(t),[unit = u"ppb"],
+        NO2(t), [unit = u"ppb"],
+        NO(t), [unit = u"ppb"],
+        H2O2(t), [unit = u"ppb"],
+        CH2O(t), [unit = u"ppb"],
     )
 
     eqs = [
diff --git a/src/wesley1989.jl b/src/wesley1989.jl
index 11fc7468..e2661a36 100644
--- a/src/wesley1989.jl
+++ b/src/wesley1989.jl
@@ -94,7 +94,8 @@ end
 # 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, iLandUse, rainOrDew::Bool)
     rs = 0.0
-    rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts))))
+    rs = ifelse((Ts >= 39.9), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts))))
+    rs = ifelse((Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts))))
     # Adjust for dew and rain (from "Effects of dew and rain" section).
     if rainOrDew
         rs *= 3
@@ -194,7 +195,7 @@ function WesleySurfaceResistance(gasData::GasData, G, Ts, θ,
     rac = obtain_value(iSeason, iLandUse,r_ac)
 
     # 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
+    correction = ifelse((Ts < 0.0), 1000.0 * exp(-Ts - 4), 0) # [s m-1] #mark
     rlux += correction
     rclx += correction
     rgsx += correction
diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl
index d143195c..a19a8859 100644
--- a/src/wet_deposition.jl
+++ b/src/wet_deposition.jl
@@ -54,7 +54,7 @@ Build Wetdeposition model
 	wd = Wetdeposition(t)
 ```
 """
-function Wetdeposition(t; name=:Wetdeposition)
+function Wetdeposition(; name=:Wetdeposition)
     params = @parameters(
         cloudFrac = 0.5, [description = "fraction of grid cell covered by clouds"],
         qrain = 0.5, [description = "rain mixing ratio"],
@@ -65,13 +65,13 @@ function Wetdeposition(t; name=:Wetdeposition)
     D = Differential(t)
 
     vars = @variables(
-        SO2(t), [unit = u"nmol/mol"],
-        O3(t), [unit = u"nmol/mol"],
-        NO2(t), [unit = u"nmol/mol"],
-        CH4(t), [unit = u"nmol/mol"],
-        CO(t), [unit = u"nmol/mol"],
-        DMS(t), [unit = u"nmol/mol"],
-        ISOP(t), [unit = u"nmol/mol"],
+        SO2(t), [unit = u"ppb"],
+        O3(t), [unit = u"ppb"],
+        NO2(t), [unit = u"ppb"],
+        CH4(t), [unit = u"ppb"],
+        CO(t), [unit = u"ppb"],
+        DMS(t), [unit = u"ppb"],
+        ISOP(t), [unit = u"ppb"],
     )
 
     eqs = [
diff --git a/test/connector_test.jl b/test/connector_test.jl
index 7ee7f918..264ce166 100644
--- a/test/connector_test.jl
+++ b/test/connector_test.jl
@@ -1,13 +1,13 @@
 using AtmosphericDeposition
-using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData
+using Test, DynamicQuantities, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData
+using ModelingToolkit:t
 
 @testset "GasChemExt" begin
     start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1))
-    @parameters t [unit = u"s"]
-    composed_ode = couple(SuperFast(t), FastJX(t), DrydepositionG(t), Wetdeposition(t))
-    tspan = (start, start+3600*24*3)
-    sys = structural_simplify(get_mtk(composed_ode))
-    @test length(states(sys)) ≈ 18
+    composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition())
+    combined_mtk = convert(ODESystem, composed_ode)
+    sys = structural_simplify(combined_mtk)
+    @test length(unknowns(sys)) ≈ 18
 
     eqs = string(equations(sys))
     wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"]
@@ -15,23 +15,21 @@ using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciDa
 end
 
 @testset "EarthSciDataExt" begin
-    @parameters t [unit = u"s"]
-
-    @parameters lat = 40
-    @parameters lon = -97
+    @parameters lat = deg2rad(40.0f0) [unit=u"rad"]
+    @parameters lon = deg2rad(-97.0f0) [unit=u"rad"]
     @parameters lev = 1
-    geosfp = GEOSFP("4x5", t)
-
-    model = couple(SuperFast(t), FastJX(t), geosfp, Wetdeposition(t), DrydepositionG(t))
+    geosfp = GEOSFP("4x5")
+    
+    model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG())
 
-    sys = structural_simplify(get_mtk(model))
-    @test length(states(sys)) ≈ 18
+    sys = structural_simplify(convert(ODESystem, model))
+    @test length(unknowns(sys)) ≈ 18
 
-    eqs = string(equations(get_mtk(model)))
+    eqs = string(equations(convert(ODESystem, model)))
     wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)"
     @test contains(eqs, wanteq)
     wanteq = "Wetdeposition₊cloudFrac(t) ~ GEOSFP₊A3cld₊CLOUD(t)"
     @test contains(eqs, wanteq)
-    wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P * PaPerhPa/(GEOSFP₊I3₊T*R)*kgperg*MW_air"
+    wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P/(GEOSFP₊I3₊T*R)*kgperg*MW_air"
     @test contains(eqs, wanteq)
 end
\ No newline at end of file
diff --git a/test/drydep_test.jl b/test/drydep_test.jl
index 935d319e..2594be5c 100644
--- a/test/drydep_test.jl
+++ b/test/drydep_test.jl
@@ -1,5 +1,5 @@
 using AtmosphericDeposition
-using Test, Unitful, ModelingToolkit
+using Test, DynamicQuantities, ModelingToolkit
 
 begin
     @parameters T [unit = u"K"]
@@ -17,19 +17,15 @@ begin
 end
 
 @testset "mfp" begin
-    @test substitute(mfp(T, P, μ), Dict(T => 298, P => 101300, μ => 1.8e-5, defaults...)) ≈ 6.512893276888993e-8
+    @test substitute(mfp(T, P, μ), Dict(T => 298, P => 101300, μ => 1.8e-5, AtmosphericDeposition.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 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, AtmosphericDeposition.So2Data, G, T, θ, iSeason, iLandUse, 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
@@ -37,7 +33,7 @@ end
     μ_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...)))
+        push!(μ_test, substitute(mu(T), Dict(T => T_[i], AtmosphericDeposition.defaults...)))
     end
     for i in 1:6
         @test (μ_test[i] - μ_list[i]) / μ_list[i] < 0.01
@@ -49,7 +45,7 @@ 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, AtmosphericDeposition.defaults...)), AtmosphericDeposition.defaults...)))
     end
     for i in 1:16
         @test (Cc_test[i] - Cc_list[i]) / Cc_list[i] < 0.03
@@ -63,8 +59,8 @@ end
     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...)))
+        push!(Cc_list, substitute(cc(Dp, T, P, μ), Dict(Dp => Dp_[i] * 1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, AtmosphericDeposition.defaults...)), AtmosphericDeposition.defaults...)))
+        push!(Vs_test, substitute(vs(Dp, ρParticle, Cc, μ), Dict(Dp => Dp_[i] * 1e-6, ρParticle => 1000, Cc => Cc_list[i], μ => 1.836522217711828e-5, AtmosphericDeposition.defaults...)))
     end
     for i in 1:4
         @test (Vs_test[i] - Vs_list[i]) / Vs_list[i] < 1
@@ -73,7 +69,7 @@ end
 
 @testset "DryDepGas" begin
     vd_true = 0.03 # m/s
-    @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, defaults...)) - vd_true) / vd_true < 0.33
+    @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, AtmosphericDeposition.defaults...)) - vd_true) / vd_true < 0.33
 end
 
 @testset "DryDepParticle" begin
@@ -81,7 +77,7 @@ end
     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], AtmosphericDeposition.defaults...)))
     end
     for i in 1:4
         @test vd_list[i] - vd_true[i] < 0.015
diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl
index 8157aeb6..aece7d7c 100644
--- a/test/wetdep_test.jl
+++ b/test/wetdep_test.jl
@@ -1,5 +1,5 @@
 using AtmosphericDeposition
-using Test, Unitful, ModelingToolkit
+using Test, DynamicQuantities, ModelingToolkit
 
 @parameters cloudFrac = 0.5
 @parameters qrain = 0.5

From 98193f8cff557f9485fa0c237f87145b9a1d2ccb Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Mon, 30 Sep 2024 16:47:23 -0500
Subject: [PATCH 05/11] Update documentation andinclude species HNO3, H2O2,
 CH2O

---
 docs/Project.toml          |  2 +-
 docs/src/dry_deposition.md | 26 +++++++++++++-------------
 docs/src/emep.md           | 16 ++++++++--------
 ext/EarthSciDataExt.jl     |  4 ++--
 ext/GasChemExt.jl          |  2 ++
 src/dry_deposition.jl      |  2 ++
 src/wet_deposition.jl      |  6 ++++++
 test/connector_test.jl     |  7 +++++--
 8 files changed, 39 insertions(+), 26 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 9bfbf718..c75f724b 100644
--- a/docs/Project.toml
+++ b/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"
\ No newline at end of file
diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md
index fc2bcfec..15fbc70a 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -3,7 +3,7 @@
 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
@@ -11,16 +11,16 @@ 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)
 
diff --git a/docs/src/emep.md b/docs/src/emep.md
index 4a9205c2..463b3027 100644
--- a/docs/src/emep.md
+++ b/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
 ```
diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl
index e7bc628c..e5d19ef3 100644
--- a/ext/EarthSciDataExt.jl
+++ b/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
 
diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl
index be976bf6..6add2acf 100644
--- a/ext/GasChemExt.jl
+++ b/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
 
diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
index 1f18c4ad..7adff712 100644
--- a/src/dry_deposition.jl
+++ b/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,
diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl
index a19a8859..83af024e 100644
--- a/src/wet_deposition.jl
+++ b/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,
diff --git a/test/connector_test.jl b/test/connector_test.jl
index 264ce166..20badd7b 100644
--- a/test/connector_test.jl
+++ b/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)"

From aa9566d23c713d37ebbc45e69ea7f5fe8f954805 Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Mon, 30 Sep 2024 16:59:04 -0500
Subject: [PATCH 06/11] Update Documentation

---
 docs/src/dry_deposition.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md
index 15fbc70a..26390908 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -51,14 +51,14 @@ The parameters in the model are:
 ```julia
 parameters(sys) # [iLandUse, z, z₀, u_star, L, ρA, G, iSeason, T, θ]
 ```
-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].
+where ```iLandUse``` is the index for the landuse, ```iSeason``` is the index for the season, ```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 = [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]
+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)
 

From 382c6b164c86cc777858fbee7f2b05c0f6c3bfde Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Mon, 30 Sep 2024 17:27:54 -0500
Subject: [PATCH 07/11] Correct the examples in documentation

---
 docs/src/dry_deposition.md |  6 +++---
 docs/src/emep.md           | 14 +++++++-------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md
index 26390908..8b449319 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -55,10 +55,10 @@ where ```iLandUse``` is the index for the landuse, ```iSeason``` is the index fo
 
 Let's run some simulation with different value for parameter ```z```. 
 ```@example 1
-@unpack O3 = sys
+@unpack O3,z = sys
 
-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]
+p1 = [z=>50]
+p2 = [z=>10]
 sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
 sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
diff --git a/docs/src/emep.md b/docs/src/emep.md
index 463b3027..865688e4 100644
--- a/docs/src/emep.md
+++ b/docs/src/emep.md
@@ -49,16 +49,16 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig
 ## Parameters
 The parameters in the model are:
 ```julia @example 1
-parameters(sys) # [cloudFrac, qrain, ρ_air, Δz]
+parameters(sys) #[ρ_air, qrain, Δz, cloudFrac]
 ```
 where ```cloudFrac``` is fraction of grid cell covered by clouds, ```qrain``` is rain mixing ratio, ```ρ_air``` is air density [kg/m3], and ```Δz``` is fall distance [m].
 
 Let's run some simulation with different value for parameter ```cloudFrac```. 
 ```@example 1
-@unpack O3 = sys
+@unpack O3,cloudFrac,qrain = sys
 
-p1 = [0.3,0.5,1.204,200]
-p2 = [0.6,0.5,1.204,200]
+p1 = [cloudFrac=>0.3]
+p2 = [cloudFrac=>0.6]
 sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
 sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
@@ -68,12 +68,12 @@ From the plot we could see that with larger cloud fraction, the wet deposition r
 
 Let's run some simulation with different value for parameter ```qrain``` 
 ```@example 1
-p3 = [0.5,0.3,1.204,200]
-p4 = [0.5,0.6,1.204,200]
+p3 = [qrain=>0.3]
+p4 = [qrain=>0.6]
 sol3 = solve(ODEProblem(sys, u0, tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0)
 sol4 = solve(ODEProblem(sys, u0, tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
-plot([sol3[O3],sol4[O3]], label = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
+plot([sol3[O3],sol4[O3]], label = ["qrain=0.3" "qrain=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```
 The graph indicates that an increase in the rain mixing ratio leads to a corresponding rise in the rate of wet deposition.
 

From d058416012925fe6954452837db5a0a461121a71 Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Tue, 1 Oct 2024 13:34:49 -0500
Subject: [PATCH 08/11] Change deposition species

---
 docs/src/dry_deposition.md | 10 +++++-----
 docs/src/emep.md           |  3 +--
 ext/GasChemExt.jl          |  5 -----
 src/dry_deposition.jl      | 14 ++++++--------
 src/wet_deposition.jl      | 20 ++++++--------------
 5 files changed, 18 insertions(+), 34 deletions(-)

diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md
index 8b449319..38959c98 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -19,9 +19,9 @@ 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,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
+tspan = (0.0, 3600*24)
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
+sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters 
 ```
 
 ```@setup 1
@@ -36,8 +36,8 @@ model = DrydepositionG()
 
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24) 
-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
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
+sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters 
 ```
 
 which we can plot as
diff --git a/docs/src/emep.md b/docs/src/emep.md
index 865688e4..d1711e90 100644
--- a/docs/src/emep.md
+++ b/docs/src/emep.md
@@ -35,8 +35,7 @@ model = Wetdeposition()
 
 sys = structural_simplify(model)
 tspan = (0.0, 3600*24)
-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, [])
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
 sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```
 
diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl
index 6add2acf..0e9f96a3 100644
--- a/ext/GasChemExt.jl
+++ b/ext/GasChemExt.jl
@@ -9,7 +9,6 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos
         c.SO2 => d.SO2,
         c.NO2 => d.NO2,
         c.O3 => d.O3,
-        c.NO => d.NO,
         c.H2O2 => d.H2O2,
         c.CH2O => d.CH2O,
     ))
@@ -22,10 +21,6 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos
         c.SO2 => d.SO2,
         c.NO2 => d.NO2,
         c.O3 => d.O3,
-        c.CH4 => d.CH4,
-        c.CO => d.CO,
-        c.DMS => d.DMS,
-        c.ISOP => d.ISOP,
         c.H2O2 => d.H2O2,
         c.CH2O => d.CH2O,
     ))
diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
index 7adff712..cc4bdda2 100644
--- a/src/dry_deposition.jl
+++ b/src/dry_deposition.jl
@@ -254,20 +254,18 @@ function DrydepositionG(; name=:DrydepositionG)
     D = Differential(t)
 
     vars = @variables( 
-        SO2(t), [unit = u"ppb"],
-        O3(t),[unit = u"ppb"],
-        NO2(t), [unit = u"ppb"],
-        NO(t), [unit = u"ppb"],
-        H2O2(t), [unit = u"ppb"],
-        CH2O(t), [unit = u"ppb"],
-        HNO3(t), [unit = u"ppb"],
+        SO2(t) = 2, [unit = u"ppb"],
+        O3(t) = 10,[unit = u"ppb"],
+        NO2(t) = 10, [unit = u"ppb"],
+        H2O2(t) = 2.34, [unit = u"ppb"],
+        CH2O(t) = 0.15, [unit = u"ppb"],
+        HNO3(t) = 10, [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
         D(HNO3) ~ -DryDepGas(z, z₀, u_star, L, ρA, Hno3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * HNO3
diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl
index 83af024e..6f69695a 100644
--- a/src/wet_deposition.jl
+++ b/src/wet_deposition.jl
@@ -65,26 +65,18 @@ function Wetdeposition(; name=:Wetdeposition)
     D = Differential(t)
 
     vars = @variables(
-        SO2(t), [unit = u"ppb"],
-        O3(t), [unit = u"ppb"],
-        NO2(t), [unit = u"ppb"],
-        CH4(t), [unit = u"ppb"],
-        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"],
+        SO2(t) = 2, [unit = u"ppb"],
+        O3(t) = 10, [unit = u"ppb"],
+        NO2(t) = 10, [unit = u"ppb"],
+        H2O2(t) = 2.34, [unit = u"ppb"],
+        HNO3(t) = 10, [unit = u"ppb"],
+        CH2O(t) = 0.15, [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
         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

From d45d5b7e94ed6db6e6bd32b27541fa8157d17dd2 Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Tue, 1 Oct 2024 13:52:56 -0500
Subject: [PATCH 09/11] Fix typo

---
 docs/src/dry_deposition.md | 4 ++--
 docs/src/emep.md           | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md
index 38959c98..31e58ad6 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -59,8 +59,8 @@ Let's run some simulation with different value for parameter ```z```.
 
 p1 = [z=>50]
 p2 = [z=>10]
-sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
-sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol1 = solve(ODEProblem(sys, [], tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol2 = solve(ODEProblem(sys, [], tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
 plot([sol1[O3],sol2[O3]], label = ["z=50m" "z=10m"], title = "Change of O3 concentration due to dry deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```
diff --git a/docs/src/emep.md b/docs/src/emep.md
index d1711e90..81f7af50 100644
--- a/docs/src/emep.md
+++ b/docs/src/emep.md
@@ -58,8 +58,8 @@ Let's run some simulation with different value for parameter ```cloudFrac```.
 
 p1 = [cloudFrac=>0.3]
 p2 = [cloudFrac=>0.6]
-sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
-sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol1 = solve(ODEProblem(sys, [], tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol2 = solve(ODEProblem(sys, [], tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
 plot([sol1[O3],sol2[O3]], label = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```
@@ -69,8 +69,8 @@ Let's run some simulation with different value for parameter ```qrain```
 ```@example 1
 p3 = [qrain=>0.3]
 p4 = [qrain=>0.6]
-sol3 = solve(ODEProblem(sys, u0, tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0)
-sol4 = solve(ODEProblem(sys, u0, tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol3 = solve(ODEProblem(sys, [], tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol4 = solve(ODEProblem(sys, [], tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
 plot([sol3[O3],sol4[O3]], label = ["qrain=0.3" "qrain=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```

From ebccc8e8acb247aa7799aefcc5acdb652b87f7de Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Fri, 1 Nov 2024 16:57:54 -0500
Subject: [PATCH 10/11] Update packages version

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index bcf6e927..4970a66a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,8 +24,8 @@ GasChemExt = "GasChem"
 
 [compat]
 DifferentialEquations = "7"
-EarthSciData = "0.9"
-EarthSciMLBase = "0.16"
+EarthSciData = "0.10"
+EarthSciMLBase = "0.19"
 GasChem = "0.7"
 ModelingToolkit = "9"
 SafeTestsets = "0.1"

From a2af6b52ae22e44b62daf461da9b6da7f1b57b4d Mon Sep 17 00:00:00 2001
From: jialinl6 <jialinl6@illinois.edu>
Date: Fri, 1 Nov 2024 17:36:55 -0500
Subject: [PATCH 11/11] Change geos-fp

---
 test/connector_test.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/test/connector_test.jl b/test/connector_test.jl
index 20badd7b..90752cf4 100644
--- a/test/connector_test.jl
+++ b/test/connector_test.jl
@@ -2,6 +2,11 @@ using AtmosphericDeposition
 using Test, DynamicQuantities, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData
 using ModelingToolkit:t
 
+domain = DomainInfo(DateTime(2022, 1, 1), DateTime(2022, 1, 3);
+    latrange=deg2rad(-85.0f0):deg2rad(2):deg2rad(85.0f0),
+    lonrange=deg2rad(-180.0f0):deg2rad(2.5):deg2rad(175.0f0),
+    levrange=1:10, dtype=Float64)
+
 @testset "GasChemExt" begin
     start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1))
     composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition())
@@ -20,7 +25,8 @@ end
     @parameters lat = deg2rad(40.0f0) [unit=u"rad"]
     @parameters lon = deg2rad(-97.0f0) [unit=u"rad"]
     @parameters lev = 1
-    geosfp = GEOSFP("4x5")
+
+    geosfp = GEOSFP("4x5", domain)
     
     model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG())