Merge branch 'EarthSciML:main' into main

Project.toml

@@ -11,18 +11,22 @@ 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]
+EarthSciDataExt = "EarthSciData"
 GasChemExt = "GasChem"
 
 [compat]
 DifferentialEquations = "7"
-EarthSciMLBase = "0.10"
+EarthSciData = "0.9"
+EarthSciMLBase = "0.15"
 GasChem = "0.6"
 IfElse = "0.1"
 ModelingToolkit = "8"
@@ -33,9 +37,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"]

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

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

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
 

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

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