diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c491b10 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.DS_Store +.vscode +output/ diff --git a/Manifest.toml b/Manifest.toml index b1f5d1f..c4a3262 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -12,6 +12,9 @@ git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" version = "3.3.1" +[[ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" + [[ArnoldiMethod]] deps = ["LinearAlgebra", "Random", "StaticArrays"] git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9" @@ -19,10 +22,7 @@ uuid = "ec485272-7323-5ecc-a04f-4719b315124d" version = "0.1.0" [[Artifacts]] -deps = ["Pkg"] -git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744" uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -version = "1.3.0" [[AxisAlgorithms]] deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] @@ -41,9 +41,9 @@ version = "1.0.1" [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "5d64be50ea9b43a89b476be773e125cef03c7cd5" +git-tree-sha1 = "dcc25ff085cf548bc8befad5ce048391a7c07d40" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.10.1" +version = "0.10.11" [[Classes]] deps = ["DataStructures", "InteractiveUtils", "MacroTools"] @@ -76,15 +76,13 @@ version = "1.0.2" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "e4e2b39db08f967cc1360951f01e8a75ec441cab" +git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.30.0" +version = "3.31.0" [[CompilerSupportLibraries_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70" +deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "0.3.4+0" [[Compose]] deps = ["Base64", "Colors", "DataStructures", "Dates", "IterTools", "JSON", "LinearAlgebra", "Measures", "Printf", "Random", "Requires", "Statistics", "UUIDs"] @@ -94,9 +92,9 @@ version = "0.9.2" [[ConstructionBase]] deps = ["LinearAlgebra"] -git-tree-sha1 = "1dc43957fb9a1574fa1b7a449e101bd1fd3a9fb7" +git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.2.1" +version = "1.3.0" [[Crayons]] git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d" @@ -104,15 +102,15 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.0.4" [[DataAPI]] -git-tree-sha1 = "dfb3b7e89e395be1e25c2ad6d7690dc29cc53b1d" +git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.6.0" +version = "1.7.0" [[DataFrames]] deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "66ee4fe515a9294a8836ef18eea7239c6ac3db5e" +git-tree-sha1 = "1dadfca11c0e08e03ab15b63aaeda55266754bad" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "1.1.1" +version = "1.2.0" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -150,16 +148,20 @@ uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" version = "0.24.18" [[DocStringExtensions]] -deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "9d4f64f79012636741cf01133158a54b24924c32" +deps = ["LibGit2"] +git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.4" +version = "0.8.5" [[DoubleFloats]] -deps = ["GenericLinearAlgebra", "GenericSchur", "LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"] -git-tree-sha1 = "cfc5657a37c2881a728d76bd14ad808ca096d601" +deps = ["GenericLinearAlgebra", "LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"] +git-tree-sha1 = "1c962cf7e75c09a5f1fbf504df7d6a06447a1129" uuid = "497a8b3b-efae-58df-a0af-a86822472b78" -version = "1.1.22" +version = "1.1.23" + +[[Downloads]] +deps = ["ArgTools", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" [[Electron]] deps = ["Base64", "FilePaths", "JSON", "Pkg", "Sockets", "URIParser", "UUIDs"] @@ -168,33 +170,27 @@ uuid = "a1bb12fb-d4d1-54b4-b10a-ee7951ef7ad3" version = "3.1.2" [[ExprTools]] -git-tree-sha1 = "10407a39b87f29d47ebaca8edbc75d7c302ff93e" +git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" -version = "0.1.3" - -[[EzXML]] -deps = ["Printf", "XML2_jll"] -git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268" -uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615" -version = "1.1.0" +version = "0.1.6" [[FFTW]] -deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"] -git-tree-sha1 = "1b48dbde42f307e48685fa9213d8b9f8c0d87594" +deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "f985af3b9f4e278b1d24434cbb546d6092fca661" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.3.2" +version = "1.4.3" [[FFTW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "5a0d4b6a22a34d17d53543bd124f4b08ed78e8b0" +git-tree-sha1 = "3676abafff7e4ff07bbd2c42b3d8201f31653dcc" uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" -version = "3.3.9+7" +version = "3.3.9+8" [[FileIO]] deps = ["Pkg", "Requires", "UUIDs"] -git-tree-sha1 = "ca36405da56db2f4730b29cd2ad2bf5869baea3c" +git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf" uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" -version = "1.9.1" +version = "1.10.1" [[FilePaths]] deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"] @@ -210,9 +206,9 @@ version = "0.9.10" [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "31939159aeb8ffad1d4d8ee44d07f8558273120a" +git-tree-sha1 = "693210145367e7685d8604aee33d9bfb85db8b31" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.11.7" +version = "0.11.9" [[FixedPointNumbers]] deps = ["Statistics"] @@ -236,12 +232,6 @@ git-tree-sha1 = "ff291c1827030ffaacaf53e3c83ed92d4d5e6fb6" uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" version = "0.2.5" -[[GenericSchur]] -deps = ["LinearAlgebra", "Printf"] -git-tree-sha1 = "372e48d7f3ced17fdc888a841bcce77be417ce57" -uuid = "c145ed77-6b09-5dd9-b285-bf645a82121e" -version = "0.5.0" - [[GlobalSensitivityAnalysis]] deps = ["DataFrames", "DataStructures", "Distributed", "Distributions", "KernelDensity", "NumericalIntegration", "ProgressMeter", "Random", "Sobol", "Statistics", "StatsBase", "VegaLite"] git-tree-sha1 = "3fb26b1bf1ee9f897341ee1a6ff186da514e8666" @@ -255,10 +245,10 @@ uuid = "a2cc645c-3eea-5389-862e-a155d0052231" version = "0.4.4" [[HTTP]] -deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508" +deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"] +git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.19" +version = "0.9.12" [[Inflate]] git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" @@ -282,10 +272,10 @@ deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[Interpolations]] -deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] -git-tree-sha1 = "1e0e51692a3a77f1eeb51bf741bdd0439ed210e7" +deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "1470c80592cf1f0a35566ee5e93c5f8221ebc33a" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.13.2" +version = "0.13.3" [[Intervals]] deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"] @@ -317,9 +307,9 @@ version = "1.0.0" [[JLD2]] deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "8ee0507881f9fbd2d6063c8befe1df4d95966d78" +git-tree-sha1 = "4813826871754cf52607e76ad37acb36ccf52719" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.8" +version = "0.4.11" [[JLLWrappers]] deps = ["Preferences"] @@ -346,24 +336,28 @@ uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" version = "0.6.3" [[LazyArtifacts]] -deps = ["Pkg"] -git-tree-sha1 = "4bb5499a1fc437342ea9ab7e319ede5a457c0968" +deps = ["Artifacts", "Pkg"] uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" -version = "1.3.0" + +[[LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" + +[[LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" [[LibGit2]] -deps = ["Printf"] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +[[LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" + [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -[[Libiconv_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "8e924324b2e9275a51407a4e06deb3455b1e359f" -uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.0+7" - [[LightGraphs]] deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6" @@ -376,9 +370,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[LogExpFunctions]] deps = ["DocStringExtensions", "LinearAlgebra"] -git-tree-sha1 = "1ba664552f1ef15325e68dc4c05c3ef8c2d5d885" +git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.2.4" +version = "0.2.5" [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -406,10 +400,8 @@ uuid = "739be429-bea8-5141-9913-cc70e7f3736d" version = "1.0.3" [[MbedTLS_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0eef589dd1c26a3ac9d753fe1a8bcad63f956fa6" +deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.16.8+1" [[Measures]] git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f" @@ -424,9 +416,9 @@ version = "0.6.7" [[Mimi]] deps = ["CSVFiles", "Classes", "DataFrames", "DataStructures", "Dates", "DelimitedFiles", "Distributions", "Electron", "FileIO", "FilePaths", "GlobalSensitivityAnalysis", "GraphPlot", "IterTools", "IteratorInterfaceExtensions", "JSON", "LightGraphs", "LinearAlgebra", "Logging", "MacroTools", "MetaGraphs", "NamedArrays", "Pkg", "ProgressMeter", "Random", "Serialization", "Statistics", "StatsBase", "StringBuilders", "TableTraits", "VegaLite"] -git-tree-sha1 = "b41925a57b0d0f59654546ed6fc5463de5596cdc" +git-tree-sha1 = "265766488a3c7c4e45a2dbeaea0c2de9e02208f8" uuid = "e4e893b0-ee5e-52ea-8111-44b3bdec128c" -version = "1.2.3" +version = "1.3.0" [[Missings]] deps = ["DataAPI"] @@ -443,11 +435,14 @@ git-tree-sha1 = "916b850daad0d46b8c71f65f719c49957e9513ed" uuid = "78c3b35d-d492-501b-9361-3d52fe80e533" version = "0.7.1" +[[MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" + [[MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "ad9b2bce6021631e0e20706d361972343a03e642" +git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01" uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -version = "0.2.19" +version = "0.2.20" [[NamedArrays]] deps = ["Combinatorics", "DataStructures", "DelimitedFiles", "InvertedIndices", "LinearAlgebra", "Random", "Requires", "SparseArrays", "Statistics"] @@ -455,6 +450,9 @@ git-tree-sha1 = "9ba8ddb0c06a08b1bad81b7120d13288e5d766fa" uuid = "86f7a689-2022-50b4-a561-43c23ac3c673" version = "0.9.5" +[[NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" + [[NodeJS]] deps = ["Pkg"] git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a" @@ -474,15 +472,15 @@ version = "0.3.3" [[OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "1381a7142eefd4cd12f052a4d2d790fe21bd1d55" +git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.9.2" +version = "1.10.2" [[OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3" +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.3+4" +version = "0.5.5+0" [[OrderedCollections]] git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" @@ -491,9 +489,9 @@ version = "1.4.1" [[PDMats]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "f82a0e71f222199de8e9eb9a09977bd0767d52a0" +git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.0" +version = "0.11.1" [[Parsers]] deps = ["Dates"] @@ -502,14 +500,14 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" version = "1.1.0" [[Pkg]] -deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[Polynomials]] deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"] -git-tree-sha1 = "3685cb1e2ccbbb9973684774f956f5c6e4673bc9" +git-tree-sha1 = "0bbfdcd8cda81b8144de4be8a67f5717e959a005" uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" -version = "2.0.12" +version = "2.0.14" [[PooledArrays]] deps = ["DataAPI", "Future"] @@ -525,9 +523,9 @@ version = "1.2.2" [[PrettyTables]] deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"] -git-tree-sha1 = "b60494adf99652d220cdef46f8a32232182cc22d" +git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6" uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" -version = "1.0.1" +version = "1.1.0" [[Printf]] deps = ["Unicode"] @@ -552,7 +550,7 @@ uuid = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" version = "0.5.5" [[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets"] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[Random]] @@ -570,9 +568,9 @@ uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" version = "1.1.1" [[Reexport]] -git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5" +git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220" uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.0.0" +version = "1.1.0" [[Requires]] deps = ["UUIDs"] @@ -582,15 +580,15 @@ version = "1.1.3" [[Rmath]] deps = ["Random", "Rmath_jll"] -git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370" +git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f" uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" -version = "0.6.1" +version = "0.7.0" [[Rmath_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "1b7bf41258f6c5c9c31df8c1ba34c1fc88674957" +git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7" uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" -version = "0.2.2+2" +version = "0.3.0+0" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -610,9 +608,9 @@ uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" [[SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] -git-tree-sha1 = "daf7aec3fe3acb2131388f93a4c409b8c7f62226" +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" -version = "0.9.3" +version = "0.9.4" [[Sobol]] deps = ["DelimitedFiles", "Random"] @@ -625,9 +623,9 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc" [[SortingAlgorithms]] deps = ["DataStructures"] -git-tree-sha1 = "2ec1962eba973f383239da22e75218565c390a96" +git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508" uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.0.0" +version = "1.0.1" [[SparseArrays]] deps = ["LinearAlgebra", "Random"] @@ -641,9 +639,9 @@ version = "1.5.1" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "42378d3bab8b4f57aa1ca443821b752850592668" +git-tree-sha1 = "a43a7b58a6e7dc933b2fa2e0ca653ccf8bb8fd0e" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.2.2" +version = "1.2.6" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -677,9 +675,7 @@ uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [[TOML]] deps = ["Dates"] -git-tree-sha1 = "44aaac2d2aec4a850302f9aa69127c74f0c3787e" uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" [[TableShowUtils]] deps = ["DataValues", "Dates", "JSON", "Markdown", "Test"] @@ -701,12 +697,16 @@ version = "1.0.1" [[Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"] -git-tree-sha1 = "aa30f8bb63f9ff3f8303a06c604c8500a69aa791" +git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.4.3" +version = "1.4.4" + +[[Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" [[Test]] -deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[TextParse]] @@ -716,10 +716,10 @@ uuid = "e0df1984-e451-5cb5-8b61-797a481e67e3" version = "1.0.1" [[TimeZones]] -deps = ["Dates", "EzXML", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"] -git-tree-sha1 = "960099aed321e05ac649c90d583d59c9309faee1" +deps = ["Dates", "Future", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"] +git-tree-sha1 = "81753f400872e5074768c9a77d4c44e70d409ef0" uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53" -version = "1.5.5" +version = "1.5.6" [[TranscodingStreams]] deps = ["Random", "Test"] @@ -733,6 +733,11 @@ git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d" uuid = "30578b45-9adc-5946-b283-645ec420af67" version = "0.4.1" +[[URIs]] +git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355" +uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" +version = "1.3.0" + [[UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" @@ -764,12 +769,6 @@ git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c" uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" version = "0.5.3" -[[XML2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "be0db24f70aae7e2b89f2f3092e93b8606d659a6" -uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.9.10+3" - [[ZipFile]] deps = ["Libdl", "Printf", "Zlib_jll"] git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7" @@ -777,7 +776,13 @@ uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" version = "0.9.3" [[Zlib_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6" +deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.11+18" + +[[nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" + +[[p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" diff --git a/Project.toml b/Project.toml index 54f4559..6bde4f1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,20 @@ +name = "MimiFAIRv2" +uui = "db6a9ac4-572f-40b4-b774-4040310d200a" +version = "0.1.0-DEV" + [deps] CSVFiles = "5d742f6a-9f54-50ce-8119-2520741973ca" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Mimi = "e4e893b0-ee5e-52ea-8111-44b3bdec128c" + +[compat] +CSVFiles = "1.0" +DataFrames = "1.0" +Mimi = "1.3" +julia = "1.4" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index a32efa1..7a75566 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,35 @@ -# MimiFAIRv2.0.jl +# MimiFAIRv2.jl -This is a work-in-progress respository for a Julia-Mimi implementation of the FAIRv2.0 simple climate model. The model description paper can be found at [FaIRv2.0.0: A Generalized Impulse Response Model for Climate Uncertainty and Future Scenario Exploration](https://gmd.copernicus.org/articles/14/3007/2021/gmd-14-3007-2021.html). +This is a work-in-progress repository for a Julia-Mimi implementation of the FAIRv2.0 simple climate model. The model description paper can be found at [FaIRv2.0.0: A Generalized Impulse Response Model for Climate Uncertainty and Future Scenario Exploration](https://gmd.copernicus.org/articles/14/3007/2021/gmd-14-3007-2021.html). -## Running Mimi-FAIRv2.0 -To run the model, execute the following code: +## Preparing the Software Environment + +To add the package to your current environment, run the following command at the julia package REPL: + +```julia +pkg> add https://github.com/FrankErrickson/MimiFAIRv2.jl.git +``` + +```julia +include("src/MimiFAIRv2.jl") # load the MimiFAIRv2 module +``` + +You probably also want to install the Mimi package into your julia environment, so that you can use some of the tools in there: + +```julia +pkg> add Mimi +``` + +## Running the Model + +The model uses the Mimi framework and it is highly recommended to read the Mimi documentation first to understand the code structure. The basic way to access a copy of the default MimiFAIRv2 model and explore the resuts is the following: ```julia -# Load the model code. -include("src/MimiFAIRv2.0.jl") +using Mimi +using MimiFAIRv2 -# Create an instance of Mimi-FAIRv2.0. -m = get_model() +# Create an instance of MimiFAIRv2. +m = MimiFAIRv2.get_model() # Run the model. run(m) @@ -33,4 +52,4 @@ The `get_model()` function currently has the following keyword arguments: \ \ -![Python vs. Julia temperature comparison](https://github.com/FrankErrickson/MimiFAIRv2.0.jl/blob/main/data/python_replication_data/Python_Mimi_FAIR2_temperature_comparison.png) +![Python vs. Julia temperature comparison](https://github.com/FrankErrickson/MimiFAIRv2.jl/blob/main/data/python_replication_data/Python_Mimi_FAIR2_temperature_comparison.png) \ No newline at end of file diff --git a/src/MimiFAIRv2.0.jl b/src/MimiFAIRv2.0.jl deleted file mode 100644 index 300bdf3..0000000 --- a/src/MimiFAIRv2.0.jl +++ /dev/null @@ -1,291 +0,0 @@ -# Load required packages. -using CSVFiles, DataFrames, Mimi - -# Load helper functions and Mimi-FAIR v2.0 model commponent files. -include("helper_functions.jl") -include("components/co2_cycle.jl") -include("components/ch4_cycle.jl") -include("components/n2o_cycle.jl") -include("components/montreal_gas_cycles.jl") -include("components/flourinated_gas_cycles.jl") -include("components/aerosol_plus_gas_cycles.jl") -include("components/radiative_forcing.jl") -include("components/temperature.jl") - - -function get_model(;emissions_forcing_scenario::String="ssp585", start_year::Int=1750, end_year::Int=2500, TCR::Float64=1.79, RWF::Float64=0.552, F2x::Float64=3.759) - - # --------------------------------------------- - # --------------------------------------------- - # Set Up Data and Parameter Values - # --------------------------------------------- - # --------------------------------------------- - - # Load emissions and forcing data and crop to appropriate model years (scenarios span 1750-2500 by default). - scenario_indices = indexin(start_year:end_year, 1750:2500) - forcing_data = DataFrame(load(joinpath(@__DIR__, "..", "data", "rcmip_"*emissions_forcing_scenario*"_effective_radiative_forcing_1750_to_2500.csv"), skiplines_begin=6))[scenario_indices,:] - emissions_data = DataFrame(load(joinpath(@__DIR__, "..", "data", "rcmip_"*emissions_forcing_scenario*"_emissions_1750_to_2500.csv"), skiplines_begin=6))[scenario_indices,:] - - # Load FAIR default gas cycle (gas) and indirect radiative forcing (irf_p) parameters. - gas_p = DataFrame(load(joinpath(@__DIR__, "..", "data", "default_gas_cycle_parameters.csv"), skiplines_begin=6)) - irf_p = DataFrame(load(joinpath(@__DIR__, "..", "data", "default_indirect_radiative_forcing_parameters.csv"), skiplines_begin=7)) - - # Load initial conditions for 1750 under the RCMIP emissions & forcing scenario. - init_gas_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7)) - init_thermal_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_thermal_conditions_1750.csv"), skiplines_begin=7)) - - # Calculate exogenous effective radiative forcing values (sum of land use, solar, and volcanic forcings). - exogenous_forcing = sum(Array(forcing_data[:, [:land_use,:volcanic,:solar]]), dims=2)[:,1] - - # Isolate specific gas parameters for convenience. - co2_p = filter(:gas_name => ==("carbon_dioxide"), gas_p) - ch4_p = filter(:gas_name => ==("methane"), gas_p) - n2o_p = filter(:gas_name => ==("nitrous_oxide"), gas_p) - montreal_gas_p = filter(:gas_group => ==("montreal"), gas_p) - flourinated_gas_p = filter(:gas_group => ==("flourinated"), gas_p) - aerosol_plus_gas_p = filter(:gas_group => ==("aerosol_plus"), gas_p) - - # Sort arrays of parameters for multiple gases so they are listed alphabetically. - sort!(flourinated_gas_p, :gas_name) - sort!(montreal_gas_p, :gas_name) - sort!(aerosol_plus_gas_p, :gas_name) - - #Isolate Montreal indirect forcing parameters and sort alphabetically. - montreal_irf_p = filter(:gas_group => ==("montreal"), irf_p) - sort!(montreal_irf_p, :gas_name) - - # Isolate initial model conditions (initial conditions currently default to year 1750). - co2_init = filter(:gas_name => ==("carbon_dioxide"), init_gas_vals) - ch4_init = filter(:gas_name => ==("methane"), init_gas_vals) - n2o_init = filter(:gas_name => ==("nitrous_oxide"), init_gas_vals) - montreal_init = filter(:gas_group => ==("montreal"), init_gas_vals) - flourinated_init = filter(:gas_group => ==("flourinated"), init_gas_vals) - aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals) - - # Sort arrays of initial conditions for multiple gases so they are listed alphabetically. - sort!(montreal_init, :gas_name) - sort!(flourinated_init, :gas_name) - sort!(aerosol_plus_init, :gas_name) - - # Extract emissions arrays for multi-gas groupings. - montreal_emissions = emissions_data[:, Symbol.(montreal_init.gas_name)] - flourinated_emissions = emissions_data[:, Symbol.(flourinated_init.gas_name)] - aerosol_plus_emissions = emissions_data[:, Symbol.(aerosol_plus_init.gas_name)] - - # Create arrays for 'a' parameter groups. - co2_a = vec(Array(co2_p[:,[:a1,:a2,:a3,:a4]])) - ch4_a = vec(Array(ch4_p[:,[:a1,:a2,:a3,:a4]])) - n2o_a = vec(Array(n2o_p[:,[:a1,:a2,:a3,:a4]])) - montreal_a = Array(montreal_gas_p[:,[:a1,:a2,:a3,:a4]]) - flourinated_a = Array(flourinated_gas_p[:,[:a1,:a2,:a3,:a4]]) - aerosol_plus_a = Array(aerosol_plus_gas_p[:,[:a1,:a2,:a3,:a4]]) - - # Create arrays for 'τ' parameter groups. - co2_τ = vec(Array(co2_p[:,[:tau1,:tau2,:tau3,:tau4]])) - ch4_τ = vec(Array(ch4_p[:,[:tau1,:tau2,:tau3,:tau4]])) - n2o_τ = vec(Array(n2o_p[:,[:tau1,:tau2,:tau3,:tau4]])) - montreal_τ = Array(montreal_gas_p[:,[:tau1,:tau2,:tau3,:tau4]]) - flourinated_τ = Array(flourinated_gas_p[:,[:tau1,:tau2,:tau3,:tau4]]) - aerosol_plus_τ = Array(aerosol_plus_gas_p[:,[:tau1,:tau2,:tau3,:tau4]]) - - # Calculate constants to approximate numerical solution for state-dependent timescale adjustment factor from Millar et al. (2017). - g0_co2, g1_co2 = calculate_g0_g1(co2_a, co2_τ) - g0_ch4, g1_ch4 = calculate_g0_g1(ch4_a, ch4_τ) - g0_n2o, g1_n2o = calculate_g0_g1(n2o_a, n2o_τ) - g0_montreal, g1_montreal = calculate_g0_g1(montreal_a, montreal_τ) - g0_flourinated, g1_flourinated = calculate_g0_g1(flourinated_a, flourinated_τ) - g0_aerosol_plus, g1_aerosol_plus = calculate_g0_g1(aerosol_plus_a, aerosol_plus_τ) - - # Calculate default thermal parameter values (defaults: transient climate response=1.79, realized warming fraction=0.552, forcing from a doubling of CO₂=3.759) - thermal_p = get_thermal_parameter_defaults(TCR=1.79, RWF=0.552, F2x=3.759) - - # Calculate thermal decay factors, defined as exp(-1/d). - thermal_decay_factors = exp.(-1.0 ./ vec(Array(thermal_p[:,[:d1,:d2,:d3]]))) - - - # --------------------------------------------- - # --------------------------------------------- - # Initialize Mimi model. - # --------------------------------------------- - # --------------------------------------------- - - # Create a Mimi model. - m = Model() - - # Set time and gas-grouping indices. - set_dimension!(m, :time, start_year:end_year) - set_dimension!(m, :montreal_gases, montreal_gas_p[:,:gas_name]) - set_dimension!(m, :flourinated_gases, flourinated_gas_p[:,:gas_name]) - set_dimension!(m, :aerosol_plus_gases, aerosol_plus_gas_p[:,:gas_name]) - - # --------------------------------------------- - # Add components to model - # --------------------------------------------- - - add_comp!(m, co2_cycle) - add_comp!(m, ch4_cycle) - add_comp!(m, n2o_cycle) - add_comp!(m, montreal_cycles) - add_comp!(m, flourinated_cycles) - add_comp!(m, aerosol_plus_cycles) - add_comp!(m, radiative_forcing) - add_comp!(m, temperature) - - # --------------------------------------------- - # Set component-specific parameters - # --------------------------------------------- - - # ---- Carbon Cycle ---- # - set_param!(m, :co2_cycle, :co2_0, co2_init.concentration[1]) - set_param!(m, :co2_cycle, :r0_co2, co2_p.r0[1]) - set_param!(m, :co2_cycle, :rU_co2, co2_p.rC[1]) - set_param!(m, :co2_cycle, :rT_co2, co2_p.rT[1]) - set_param!(m, :co2_cycle, :rA_co2, co2_p.rA[1]) - set_param!(m, :co2_cycle, :GU_co2_0, co2_init.cumulative_uptake[1]) - set_param!(m, :co2_cycle, :g0_co2, g0_co2) - set_param!(m, :co2_cycle, :g1_co2, g1_co2) - set_param!(m, :co2_cycle, :R0_co2, vec(Array(co2_init[:,[:pool1,:pool2,:pool3,:pool4]]))) - set_param!(m, :co2_cycle, :emiss2conc_co2, co2_p.emis2conc[1]) - set_param!(m, :co2_cycle, :a_co2, co2_a) - set_param!(m, :co2_cycle, :τ_co2, co2_τ) - set_param!(m, :co2_cycle, :E_co2, emissions_data.carbon_dioxide) - - # ---- Methane Cycle ---- # - set_param!(m, :ch4_cycle, :ch4_0, ch4_init.concentration[1]) - set_param!(m, :ch4_cycle, :r0_ch4, ch4_p.r0[1]) - set_param!(m, :ch4_cycle, :rU_ch4, ch4_p.rC[1]) - set_param!(m, :ch4_cycle, :rT_ch4, ch4_p.rT[1]) - set_param!(m, :ch4_cycle, :rA_ch4, ch4_p.rA[1]) - set_param!(m, :ch4_cycle, :GU_ch4_0, ch4_init.cumulative_uptake[1]) - set_param!(m, :ch4_cycle, :g0_ch4, g0_ch4) - set_param!(m, :ch4_cycle, :g1_ch4, g1_ch4) - set_param!(m, :ch4_cycle, :R0_ch4, vec(Array(ch4_init[:,[:pool1,:pool2,:pool3,:pool4]]))) - set_param!(m, :ch4_cycle, :emiss2conc_ch4, ch4_p.emis2conc[1]) - set_param!(m, :ch4_cycle, :a_ch4, ch4_a) - set_param!(m, :ch4_cycle, :τ_ch4, ch4_τ) - set_param!(m, :ch4_cycle, :E_ch4, emissions_data.methane) - - # ---- Nitrous Oxide Cycle ---- # - set_param!(m, :n2o_cycle, :n2o_0, n2o_init.concentration[1]) - set_param!(m, :n2o_cycle, :r0_n2o, n2o_p.r0[1]) - set_param!(m, :n2o_cycle, :rU_n2o, n2o_p.rC[1]) - set_param!(m, :n2o_cycle, :rT_n2o, n2o_p.rT[1]) - set_param!(m, :n2o_cycle, :rA_n2o, n2o_p.rA[1]) - set_param!(m, :n2o_cycle, :GU_n2o_0, n2o_init.cumulative_uptake[1]) - set_param!(m, :n2o_cycle, :g0_n2o, g0_n2o) - set_param!(m, :n2o_cycle, :g1_n2o, g1_n2o) - set_param!(m, :n2o_cycle, :R0_n2o, vec(Array(n2o_init[:,[:pool1,:pool2,:pool3,:pool4]]))) - set_param!(m, :n2o_cycle, :emiss2conc_n2o, n2o_p.emis2conc[1]) - set_param!(m, :n2o_cycle, :a_n2o, n2o_a) - set_param!(m, :n2o_cycle, :τ_n2o, n2o_τ) - set_param!(m, :n2o_cycle, :E_n2o, emissions_data.nitrous_oxide) - - # ---- Flourinated Gas Cycles ---- # - set_param!(m, :flourinated_cycles, :flourinated_0, flourinated_init[:, :concentration]) - set_param!(m, :flourinated_cycles, :r0_flourinated, flourinated_gas_p[:,:r0]) - set_param!(m, :flourinated_cycles, :rU_flourinated, flourinated_gas_p[:,:rC]) - set_param!(m, :flourinated_cycles, :rT_flourinated, flourinated_gas_p[:,:rT]) - set_param!(m, :flourinated_cycles, :rA_flourinated, flourinated_gas_p[:,:rA]) - set_param!(m, :flourinated_cycles, :GU_flourinated_0, flourinated_init[:, :cumulative_uptake]) - set_param!(m, :flourinated_cycles, :g0_flourinated, g0_flourinated) - set_param!(m, :flourinated_cycles, :g1_flourinated, g1_flourinated) - set_param!(m, :flourinated_cycles, :R0_flourinated, Array(flourinated_init[:,[:pool1,:pool2,:pool3,:pool4]])) - set_param!(m, :flourinated_cycles, :emiss2conc_flourinated, flourinated_gas_p[:,:emis2conc]) - set_param!(m, :flourinated_cycles, :a_flourinated, flourinated_a) - set_param!(m, :flourinated_cycles, :τ_flourinated, flourinated_τ) - set_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions)) - - # ---- Montreal Protocol Gas Cycles ---- # - set_param!(m, :montreal_cycles, :montreal_0, montreal_init[:, :concentration]) - set_param!(m, :montreal_cycles, :r0_montreal, montreal_gas_p[:,:r0]) - set_param!(m, :montreal_cycles, :rU_montreal, montreal_gas_p[:,:rC]) - set_param!(m, :montreal_cycles, :rT_montreal, montreal_gas_p[:,:rT]) - set_param!(m, :montreal_cycles, :rA_montreal, montreal_gas_p[:,:rA]) - set_param!(m, :montreal_cycles, :GU_montreal_0, montreal_init[:, :cumulative_uptake]) - set_param!(m, :montreal_cycles, :g0_montreal, g0_montreal) - set_param!(m, :montreal_cycles, :g1_montreal, g1_montreal) - set_param!(m, :montreal_cycles, :R0_montreal, Array(montreal_init[:,[:pool1,:pool2,:pool3,:pool4]])) - set_param!(m, :montreal_cycles, :emiss2conc_montreal, montreal_gas_p[:,:emis2conc]) - set_param!(m, :montreal_cycles, :a_montreal, montreal_a) - set_param!(m, :montreal_cycles, :τ_montreal, montreal_τ) - set_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions)) - - # ---- Tropospheric Ozone Precursors, Aerosols, & Reactive Gas Cycles (Aerosol+) ---- # - set_param!(m, :aerosol_plus_cycles, :aerosol_plus_0, aerosol_plus_init[:, :concentration]) - set_param!(m, :aerosol_plus_cycles, :r0_aerosol_plus, aerosol_plus_gas_p[:,:r0]) - set_param!(m, :aerosol_plus_cycles, :rU_aerosol_plus, aerosol_plus_gas_p[:,:rC]) - set_param!(m, :aerosol_plus_cycles, :rT_aerosol_plus, aerosol_plus_gas_p[:,:rT]) - set_param!(m, :aerosol_plus_cycles, :rA_aerosol_plus, aerosol_plus_gas_p[:,:rA]) - set_param!(m, :aerosol_plus_cycles, :GU_aerosol_plus_0, aerosol_plus_init[:, :cumulative_uptake]) - set_param!(m, :aerosol_plus_cycles, :g0_aerosol_plus, g0_aerosol_plus) - set_param!(m, :aerosol_plus_cycles, :g1_aerosol_plus, g1_aerosol_plus) - set_param!(m, :aerosol_plus_cycles, :R0_aerosol_plus, Array(aerosol_plus_init[:,[:pool1,:pool2,:pool3,:pool4]])) - set_param!(m, :aerosol_plus_cycles, :emiss2conc_aerosol_plus, aerosol_plus_gas_p[:,:emis2conc]) - set_param!(m, :aerosol_plus_cycles, :a_aerosol_plus, aerosol_plus_a) - set_param!(m, :aerosol_plus_cycles, :τ_aerosol_plus, aerosol_plus_τ) - set_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions)) - - # ---- Radiative Forcing ---- # - set_param!(m, :radiative_forcing, :co2_f, vec(Array(co2_p[:, [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :ch4_f, vec(Array(ch4_p[:, [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :ch4_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="methane|o3"), [:f1, :f2, :f3]])) ) - set_param!(m, :radiative_forcing, :ch4_h2o_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="methane|strat_h2o"), [:f1, :f2, :f3]])) ) - set_param!(m, :radiative_forcing, :n2o_f, vec(Array(n2o_p[:, [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :n2o_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nitrous_oxide|o3"), [:f1, :f2, :f3]])) ) - set_param!(m, :radiative_forcing, :montreal_f, Array(montreal_gas_p[:,[:f1, :f2, :f3]])) - set_param!(m, :radiative_forcing, :montreal_ind_f, Array(montreal_irf_p[:,[:f1, :f2, :f3]])) - set_param!(m, :radiative_forcing, :flourinated_f, Array(flourinated_gas_p[:,[:f1, :f2, :f3]])) - set_param!(m, :radiative_forcing, :bc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="bc"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :bc_snow_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="bc|bc_on_snow"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :bc_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="bc|aci"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :co_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="co"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :co_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="co|o3"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nh3_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nh3"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nmvoc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nmvoc"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nmvoc_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nmvoc|o3"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nox_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nox"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nox_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nox|o3"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nox_avi_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nox_avi"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :nox_avi_contrails_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nox_avi|contrails"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :oc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="oc"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :oc_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="oc|aci"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :so2_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="so2"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :so2_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="so2|aci"), [:f1, :f2, :f3]]))) - set_param!(m, :radiative_forcing, :exogenous_RF, exogenous_forcing) - - # ---- Global Temperature Anomaly ---- # - set_param!(m, :temperature, :Tj_0, vec(Array(init_thermal_vals[:,[:thermal_box_1,:thermal_box_2,:thermal_box_3]]))) - set_param!(m, :temperature, :T_0, init_thermal_vals.global_temp_anomaly[1]) - set_param!(m, :temperature, :q, vec(Array(thermal_p[:,[:q1,:q2,:q3]]))) - set_param!(m, :temperature, :decay_factor, thermal_decay_factors) - - # --- Set Parameters Common to Multiple Components ---- # - set_param!(m, :co2_pi, co2_p.PI_conc[1]) - set_param!(m, :ch4_pi, ch4_p.PI_conc[1]) - set_param!(m, :n2o_pi, n2o_p.PI_conc[1]) - set_param!(m, :montreal_pi, montreal_gas_p[:,:PI_conc]) - set_param!(m, :flourinated_pi, flourinated_gas_p[:,:PI_conc]) - set_param!(m, :aerosol_plus_pi, aerosol_plus_gas_p[:,:PI_conc]) - - # --------------------------------------------- - # Create Connections Between Mimi Components - # --------------------------------------------- - - # Syntax is :component_needing_a_parameter_input => :name_of_that_parameter, :component_calculating_required_values => :name_of_variable_output - connect_param!(m, :montreal_cycles => :Tj, :temperature => :Tj) - connect_param!(m, :flourinated_cycles => :Tj, :temperature => :Tj) - connect_param!(m, :aerosol_plus_cycles => :Tj, :temperature => :Tj) - connect_param!(m, :co2_cycle => :Tj, :temperature => :Tj) - connect_param!(m, :ch4_cycle => :Tj, :temperature => :Tj) - connect_param!(m, :n2o_cycle => :Tj, :temperature => :Tj) - connect_param!(m, :radiative_forcing => :co2_conc, :co2_cycle => :co2) - connect_param!(m, :radiative_forcing => :ch4_conc, :ch4_cycle => :ch4) - connect_param!(m, :radiative_forcing => :n2o_conc, :n2o_cycle => :n2o) - connect_param!(m, :radiative_forcing => :montreal_conc, :montreal_cycles => :montreal_conc) - connect_param!(m, :radiative_forcing => :flourinated_conc, :flourinated_cycles => :flourinated_conc) - connect_param!(m, :radiative_forcing => :aerosol_plus_conc, :aerosol_plus_cycles => :aerosol_plus_conc) - connect_param!(m, :temperature => :F, :radiative_forcing => :total_RF) - - # Return model. - return m -end diff --git a/src/MimiFAIRv2.jl b/src/MimiFAIRv2.jl new file mode 100644 index 0000000..b913b29 --- /dev/null +++ b/src/MimiFAIRv2.jl @@ -0,0 +1,334 @@ +module MimiFAIRv2 + +# Load required packages. +using CSVFiles, DataFrames, Mimi + +# Load helper functions and MimiFAIRv2 model commponent files. +include("helper_functions.jl") +include("components/co2_cycle.jl") +include("components/ch4_cycle.jl") +include("components/n2o_cycle.jl") +include("components/montreal_gas_cycles.jl") +include("components/flourinated_gas_cycles.jl") +include("components/aerosol_plus_gas_cycles.jl") +include("components/radiative_forcing.jl") +include("components/temperature.jl") + +""" + get_model(;emissions_forcing_scenario::String="ssp585", start_year::Int=1750, + end_year::Int=2500, TCR::Float64=1.79, RWF::Float64=0.552, + F2x::Float64=3.759) + +Return a constructed model with default settings for emissions forcing scenario +(default to ssp585), start year (default to 1750), end year (default to 2500), +transient climate response (default to 1.79), realized warming fraction (default +to 5.22), and forcing from a doubling of CO₂ (default to 3.759). +""" +function get_model(;emissions_forcing_scenario::String="ssp585", start_year::Int=1750, end_year::Int=2500, TCR::Float64=1.79, RWF::Float64=0.552, F2x::Float64=3.759) + + # TODO - should this warn or error? + if start_year !== 1750 + @warn "Model should not be set to start with a year differing from 1750." + end + + # --------------------------------------------- + # --------------------------------------------- + # Set Up Data and Parameter Values + # --------------------------------------------- + # --------------------------------------------- + + # Load emissions and forcing data and crop to appropriate model years (scenarios span 1750-2500 by default). + scenario_indices = indexin(start_year:end_year, 1750:2500) + forcing_data = DataFrame(load(joinpath(@__DIR__, "..", "data", "rcmip_"*emissions_forcing_scenario*"_effective_radiative_forcing_1750_to_2500.csv"), skiplines_begin=6))[scenario_indices,:] + emissions_data = DataFrame(load(joinpath(@__DIR__, "..", "data", "rcmip_"*emissions_forcing_scenario*"_emissions_1750_to_2500.csv"), skiplines_begin=6))[scenario_indices,:] + + # Load FAIR default gas cycle (gas) and indirect radiative forcing (irf_p) parameters. + gas_p = DataFrame(load(joinpath(@__DIR__, "..", "data", "default_gas_cycle_parameters.csv"), skiplines_begin=6)) + irf_p = DataFrame(load(joinpath(@__DIR__, "..", "data", "default_indirect_radiative_forcing_parameters.csv"), skiplines_begin=7)) + + # Load initial conditions for 1750 under the RCMIP emissions & forcing scenario. + init_gas_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7)) + init_thermal_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_thermal_conditions_1750.csv"), skiplines_begin=7)) + + # Calculate exogenous effective radiative forcing values (sum of land use, solar, and volcanic forcings). + exogenous_forcing = sum(Array(forcing_data[:, [:land_use,:volcanic,:solar]]), dims=2)[:,1] + + # Isolate specific gas parameters for convenience. + co2_p = filter(:gas_name => ==("carbon_dioxide"), gas_p) + ch4_p = filter(:gas_name => ==("methane"), gas_p) + n2o_p = filter(:gas_name => ==("nitrous_oxide"), gas_p) + montreal_gas_p = filter(:gas_group => ==("montreal"), gas_p) + flourinated_gas_p = filter(:gas_group => ==("flourinated"), gas_p) + aerosol_plus_gas_p = filter(:gas_group => ==("aerosol_plus"), gas_p) + + # Sort arrays of parameters for multiple gases so they are listed alphabetically. + sort!(flourinated_gas_p, :gas_name) + sort!(montreal_gas_p, :gas_name) + sort!(aerosol_plus_gas_p, :gas_name) + + #Isolate Montreal indirect forcing parameters and sort alphabetically. + montreal_irf_p = filter(:gas_group => ==("montreal"), irf_p) + sort!(montreal_irf_p, :gas_name) + + # Isolate initial model conditions (initial conditions currently default to year 1750). + co2_init = filter(:gas_name => ==("carbon_dioxide"), init_gas_vals) + ch4_init = filter(:gas_name => ==("methane"), init_gas_vals) + n2o_init = filter(:gas_name => ==("nitrous_oxide"), init_gas_vals) + montreal_init = filter(:gas_group => ==("montreal"), init_gas_vals) + flourinated_init = filter(:gas_group => ==("flourinated"), init_gas_vals) + aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals) + + # Sort arrays of initial conditions for multiple gases so they are listed alphabetically. + sort!(montreal_init, :gas_name) + sort!(flourinated_init, :gas_name) + sort!(aerosol_plus_init, :gas_name) + + # Extract emissions arrays for multi-gas groupings. + montreal_emissions = emissions_data[:, Symbol.(montreal_init.gas_name)] + flourinated_emissions = emissions_data[:, Symbol.(flourinated_init.gas_name)] + aerosol_plus_emissions = emissions_data[:, Symbol.(aerosol_plus_init.gas_name)] + + # Create helper arrays of indices to use for pulling parameter groups. + a_idxs = [:a1,:a2,:a3,:a4] + τ_idxs = [:tau1,:tau2,:tau3,:tau4] + + # Create arrays for 'a' parameter groups. + co2_a = vec(Array(co2_p[:, a_idxs])) + ch4_a = vec(Array(ch4_p[:, a_idxs])) + n2o_a = vec(Array(n2o_p[:, a_idxs])) + montreal_a = Array(montreal_gas_p[:, a_idxs]) + flourinated_a = Array(flourinated_gas_p[:, a_idxs]) + aerosol_plus_a = Array(aerosol_plus_gas_p[:, a_idxs]) + + # Create arrays for 'τ' parameter groups. + co2_τ = vec(Array(co2_p[:, τ_idxs])) + ch4_τ = vec(Array(ch4_p[:, τ_idxs])) + n2o_τ = vec(Array(n2o_p[:, τ_idxs])) + montreal_τ = Array(montreal_gas_p[:, τ_idxs]) + flourinated_τ = Array(flourinated_gas_p[:, τ_idxs]) + aerosol_plus_τ = Array(aerosol_plus_gas_p[:, τ_idxs]) + + # Calculate constants to approximate numerical solution for state-dependent timescale adjustment factor from Millar et al. (2017). + g0_co2, g1_co2 = calculate_g0_g1(co2_a, co2_τ) + g0_ch4, g1_ch4 = calculate_g0_g1(ch4_a, ch4_τ) + g0_n2o, g1_n2o = calculate_g0_g1(n2o_a, n2o_τ) + g0_montreal, g1_montreal = calculate_g0_g1(montreal_a, montreal_τ) + g0_flourinated, g1_flourinated = calculate_g0_g1(flourinated_a, flourinated_τ) + g0_aerosol_plus, g1_aerosol_plus = calculate_g0_g1(aerosol_plus_a, aerosol_plus_τ) + + # Calculate default thermal parameter values + # Defaults pulled from get_model function defaults are are as follows: + # transient climate response = 1.79 + # realized warming fraction = 0.552 + # forcing from a doubling of CO₂ = 3.759 + thermal_p = get_thermal_parameter_defaults(TCR = TCR, RWF= RWF, F2x = F2x) + + # Calculate thermal decay factors, defined as exp(-1/d). + thermal_decay_factors = exp.(-1.0 ./ vec(Array(thermal_p[:,[:d1,:d2,:d3]]))) + + + # --------------------------------------------- + # --------------------------------------------- + # Initialize Mimi model. + # --------------------------------------------- + # --------------------------------------------- + + # Create a Mimi model. + m = Model() + + # Set time and gas-grouping indices. + set_dimension!(m, :time, start_year:end_year) + set_dimension!(m, :montreal_gases, montreal_gas_p[:,:gas_name]) + set_dimension!(m, :flourinated_gases, flourinated_gas_p[:,:gas_name]) + set_dimension!(m, :aerosol_plus_gases, aerosol_plus_gas_p[:,:gas_name]) + + # --------------------------------------------- + # Add components to model + # --------------------------------------------- + + add_comp!(m, co2_cycle) + add_comp!(m, ch4_cycle) + add_comp!(m, n2o_cycle) + add_comp!(m, montreal_cycles) + add_comp!(m, flourinated_cycles) + add_comp!(m, aerosol_plus_cycles) + add_comp!(m, radiative_forcing) + add_comp!(m, temperature) + + # --------------------------------------------- + # Set component-specific parameters + # --------------------------------------------- + + # ---- Carbon Cycle ---- # + update_param!(m, :co2_cycle, :co2_0, co2_init.concentration[1]) + update_param!(m, :co2_cycle, :r0_co2, co2_p.r0[1]) + update_param!(m, :co2_cycle, :rU_co2, co2_p.rC[1]) + update_param!(m, :co2_cycle, :rT_co2, co2_p.rT[1]) + update_param!(m, :co2_cycle, :rA_co2, co2_p.rA[1]) + update_param!(m, :co2_cycle, :GU_co2_0, co2_init.cumulative_uptake[1]) + update_param!(m, :co2_cycle, :g0_co2, g0_co2) + update_param!(m, :co2_cycle, :g1_co2, g1_co2) + update_param!(m, :co2_cycle, :R0_co2, vec(Array(co2_init[:,[:pool1,:pool2,:pool3,:pool4]]))) + update_param!(m, :co2_cycle, :emiss2conc_co2, co2_p.emis2conc[1]) + update_param!(m, :co2_cycle, :a_co2, co2_a) + update_param!(m, :co2_cycle, :τ_co2, co2_τ) + update_param!(m, :co2_cycle, :E_co2, emissions_data.carbon_dioxide) + + # ---- Methane Cycle ---- # + update_param!(m, :ch4_cycle, :ch4_0, ch4_init.concentration[1]) + update_param!(m, :ch4_cycle, :r0_ch4, ch4_p.r0[1]) + update_param!(m, :ch4_cycle, :rU_ch4, ch4_p.rC[1]) + update_param!(m, :ch4_cycle, :rT_ch4, ch4_p.rT[1]) + update_param!(m, :ch4_cycle, :rA_ch4, ch4_p.rA[1]) + update_param!(m, :ch4_cycle, :GU_ch4_0, ch4_init.cumulative_uptake[1]) + update_param!(m, :ch4_cycle, :g0_ch4, g0_ch4) + update_param!(m, :ch4_cycle, :g1_ch4, g1_ch4) + update_param!(m, :ch4_cycle, :R0_ch4, vec(Array(ch4_init[:,[:pool1,:pool2,:pool3,:pool4]]))) + update_param!(m, :ch4_cycle, :emiss2conc_ch4, ch4_p.emis2conc[1]) + update_param!(m, :ch4_cycle, :a_ch4, ch4_a) + update_param!(m, :ch4_cycle, :τ_ch4, ch4_τ) + update_param!(m, :ch4_cycle, :E_ch4, emissions_data.methane) + + # ---- Nitrous Oxide Cycle ---- # + update_param!(m, :n2o_cycle, :n2o_0, n2o_init.concentration[1]) + update_param!(m, :n2o_cycle, :r0_n2o, n2o_p.r0[1]) + update_param!(m, :n2o_cycle, :rU_n2o, n2o_p.rC[1]) + update_param!(m, :n2o_cycle, :rT_n2o, n2o_p.rT[1]) + update_param!(m, :n2o_cycle, :rA_n2o, n2o_p.rA[1]) + update_param!(m, :n2o_cycle, :GU_n2o_0, n2o_init.cumulative_uptake[1]) + update_param!(m, :n2o_cycle, :g0_n2o, g0_n2o) + update_param!(m, :n2o_cycle, :g1_n2o, g1_n2o) + update_param!(m, :n2o_cycle, :R0_n2o, vec(Array(n2o_init[:,[:pool1,:pool2,:pool3,:pool4]]))) + update_param!(m, :n2o_cycle, :emiss2conc_n2o, n2o_p.emis2conc[1]) + update_param!(m, :n2o_cycle, :a_n2o, n2o_a) + update_param!(m, :n2o_cycle, :τ_n2o, n2o_τ) + update_param!(m, :n2o_cycle, :E_n2o, emissions_data.nitrous_oxide) + + # ---- Flourinated Gas Cycles ---- # + update_param!(m, :flourinated_cycles, :flourinated_0, flourinated_init[:, :concentration]) + update_param!(m, :flourinated_cycles, :r0_flourinated, flourinated_gas_p[:,:r0]) + update_param!(m, :flourinated_cycles, :rU_flourinated, flourinated_gas_p[:,:rC]) + update_param!(m, :flourinated_cycles, :rT_flourinated, flourinated_gas_p[:,:rT]) + update_param!(m, :flourinated_cycles, :rA_flourinated, flourinated_gas_p[:,:rA]) + update_param!(m, :flourinated_cycles, :GU_flourinated_0, flourinated_init[:, :cumulative_uptake]) + update_param!(m, :flourinated_cycles, :g0_flourinated, g0_flourinated) + update_param!(m, :flourinated_cycles, :g1_flourinated, g1_flourinated) + update_param!(m, :flourinated_cycles, :R0_flourinated, Array(flourinated_init[:,[:pool1,:pool2,:pool3,:pool4]])) + update_param!(m, :flourinated_cycles, :emiss2conc_flourinated, flourinated_gas_p[:,:emis2conc]) + update_param!(m, :flourinated_cycles, :a_flourinated, flourinated_a) + update_param!(m, :flourinated_cycles, :τ_flourinated, flourinated_τ) + update_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions)) + + # ---- Montreal Protocol Gas Cycles ---- # + update_param!(m, :montreal_cycles, :montreal_0, montreal_init[:, :concentration]) + update_param!(m, :montreal_cycles, :r0_montreal, montreal_gas_p[:,:r0]) + update_param!(m, :montreal_cycles, :rU_montreal, montreal_gas_p[:,:rC]) + update_param!(m, :montreal_cycles, :rT_montreal, montreal_gas_p[:,:rT]) + update_param!(m, :montreal_cycles, :rA_montreal, montreal_gas_p[:,:rA]) + update_param!(m, :montreal_cycles, :GU_montreal_0, montreal_init[:, :cumulative_uptake]) + update_param!(m, :montreal_cycles, :g0_montreal, g0_montreal) + update_param!(m, :montreal_cycles, :g1_montreal, g1_montreal) + update_param!(m, :montreal_cycles, :R0_montreal, Array(montreal_init[:,[:pool1,:pool2,:pool3,:pool4]])) + update_param!(m, :montreal_cycles, :emiss2conc_montreal, montreal_gas_p[:,:emis2conc]) + update_param!(m, :montreal_cycles, :a_montreal, montreal_a) + update_param!(m, :montreal_cycles, :τ_montreal, montreal_τ) + update_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions)) + + # ---- Tropospheric Ozone Precursors, Aerosols, & Reactive Gas Cycles (Aerosol+) ---- # + update_param!(m, :aerosol_plus_cycles, :aerosol_plus_0, aerosol_plus_init[:, :concentration]) + update_param!(m, :aerosol_plus_cycles, :r0_aerosol_plus, aerosol_plus_gas_p[:,:r0]) + update_param!(m, :aerosol_plus_cycles, :rU_aerosol_plus, aerosol_plus_gas_p[:,:rC]) + update_param!(m, :aerosol_plus_cycles, :rT_aerosol_plus, aerosol_plus_gas_p[:,:rT]) + update_param!(m, :aerosol_plus_cycles, :rA_aerosol_plus, aerosol_plus_gas_p[:,:rA]) + update_param!(m, :aerosol_plus_cycles, :GU_aerosol_plus_0, aerosol_plus_init[:, :cumulative_uptake]) + update_param!(m, :aerosol_plus_cycles, :g0_aerosol_plus, g0_aerosol_plus) + update_param!(m, :aerosol_plus_cycles, :g1_aerosol_plus, g1_aerosol_plus) + update_param!(m, :aerosol_plus_cycles, :R0_aerosol_plus, Array(aerosol_plus_init[:,[:pool1,:pool2,:pool3,:pool4]])) + update_param!(m, :aerosol_plus_cycles, :emiss2conc_aerosol_plus, aerosol_plus_gas_p[:,:emis2conc]) + update_param!(m, :aerosol_plus_cycles, :a_aerosol_plus, aerosol_plus_a) + update_param!(m, :aerosol_plus_cycles, :τ_aerosol_plus, aerosol_plus_τ) + update_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions)) + + # ---- Radiative Forcing ---- # + update_param!(m, :radiative_forcing, :co2_f, vec(Array(co2_p[:, [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :ch4_f, vec(Array(ch4_p[:, [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :ch4_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="methane|o3"), [:f1, :f2, :f3]])) ) + update_param!(m, :radiative_forcing, :ch4_h2o_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="methane|strat_h2o"), [:f1, :f2, :f3]])) ) + update_param!(m, :radiative_forcing, :n2o_f, vec(Array(n2o_p[:, [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :n2o_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nitrous_oxide|o3"), [:f1, :f2, :f3]])) ) + update_param!(m, :radiative_forcing, :montreal_f, Array(montreal_gas_p[:,[:f1, :f2, :f3]])) + update_param!(m, :radiative_forcing, :montreal_ind_f, Array(montreal_irf_p[:,[:f1, :f2, :f3]])) + update_param!(m, :radiative_forcing, :flourinated_f, Array(flourinated_gas_p[:,[:f1, :f2, :f3]])) + update_param!(m, :radiative_forcing, :bc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="bc"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :bc_snow_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="bc|bc_on_snow"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :bc_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="bc|aci"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :co_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="co"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :co_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="co|o3"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nh3_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nh3"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nmvoc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nmvoc"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nmvoc_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nmvoc|o3"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nox_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nox"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nox_o3_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nox|o3"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nox_avi_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="nox_avi"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :nox_avi_contrails_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="nox_avi|contrails"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :oc_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="oc"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :oc_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="oc|aci"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :so2_f, vec(Array(aerosol_plus_gas_p[(aerosol_plus_gas_p.gas_name.=="so2"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :so2_aci_f, vec(Array(irf_p[(irf_p.indirect_forcing_effect.=="so2|aci"), [:f1, :f2, :f3]]))) + update_param!(m, :radiative_forcing, :exogenous_RF, exogenous_forcing) + + # ---- Global Temperature Anomaly ---- # + update_param!(m, :temperature, :Tj_0, vec(Array(init_thermal_vals[:,[:thermal_box_1,:thermal_box_2,:thermal_box_3]]))) + update_param!(m, :temperature, :T_0, init_thermal_vals.global_temp_anomaly[1]) + update_param!(m, :temperature, :q, vec(Array(thermal_p[:,[:q1,:q2,:q3]]))) + update_param!(m, :temperature, :decay_factor, thermal_decay_factors) + + # --- Set Parameters Common to Multiple Components ---- # + add_shared_param!(m, :shared_co2_pi, co2_p.PI_conc[1]) + connect_param!(m, :co2_cycle, :co2_pi, :shared_co2_pi) + connect_param!(m, :radiative_forcing, :co2_pi, :shared_co2_pi) + + add_shared_param!(m, :shared_ch4_pi, ch4_p.PI_conc[1]) + connect_param!(m, :ch4_cycle, :ch4_pi, :shared_ch4_pi) + connect_param!(m, :radiative_forcing, :ch4_pi, :shared_ch4_pi) + + add_shared_param!(m, :shared_n2o_pi, n2o_p.PI_conc[1]) + connect_param!(m, :n2o_cycle, :n2o_pi, :shared_n2o_pi) + connect_param!(m, :radiative_forcing, :n2o_pi, :shared_n2o_pi) + + add_shared_param!(m, :shared_montreal_pi, montreal_gas_p[:,:PI_conc], dims = [:montreal_gases]) + connect_param!(m, :montreal_cycles, :montreal_pi, :shared_montreal_pi) + connect_param!(m, :radiative_forcing, :montreal_pi, :shared_montreal_pi) + + add_shared_param!(m, :shared_flourinated_pi, flourinated_gas_p[:,:PI_conc], dims = [:flourinated_gases]) + connect_param!(m, :flourinated_cycles, :flourinated_pi, :shared_flourinated_pi) + connect_param!(m, :radiative_forcing, :flourinated_pi, :shared_flourinated_pi) + + add_shared_param!(m, :shared_aerosol_plus_pi, aerosol_plus_gas_p[:,:PI_conc], dims = [:aerosol_plus_gases]) + connect_param!(m, :aerosol_plus_cycles, :aerosol_plus_pi, :shared_aerosol_plus_pi) + connect_param!(m, :radiative_forcing, :aerosol_plus_pi, :shared_aerosol_plus_pi) + + # --------------------------------------------- + # Create Connections Between Mimi Components + # --------------------------------------------- + + # Syntax is :component_needing_a_parameter_input => :name_of_that_parameter, :component_calculating_required_values => :name_of_variable_output + connect_param!(m, :montreal_cycles => :Tj, :temperature => :Tj) + connect_param!(m, :flourinated_cycles => :Tj, :temperature => :Tj) + connect_param!(m, :aerosol_plus_cycles => :Tj, :temperature => :Tj) + connect_param!(m, :co2_cycle => :Tj, :temperature => :Tj) + connect_param!(m, :ch4_cycle => :Tj, :temperature => :Tj) + connect_param!(m, :n2o_cycle => :Tj, :temperature => :Tj) + connect_param!(m, :radiative_forcing => :co2_conc, :co2_cycle => :co2) + connect_param!(m, :radiative_forcing => :ch4_conc, :ch4_cycle => :ch4) + connect_param!(m, :radiative_forcing => :n2o_conc, :n2o_cycle => :n2o) + connect_param!(m, :radiative_forcing => :montreal_conc, :montreal_cycles => :montreal_conc) + connect_param!(m, :radiative_forcing => :flourinated_conc, :flourinated_cycles => :flourinated_conc) + connect_param!(m, :radiative_forcing => :aerosol_plus_conc, :aerosol_plus_cycles => :aerosol_plus_conc) + connect_param!(m, :temperature => :F, :radiative_forcing => :total_RF) + + # Return model. + return m +end + +end # module diff --git a/src/components/aerosol_plus_gas_cycles.jl b/src/components/aerosol_plus_gas_cycles.jl index 1505d26..e5ff6d7 100644 --- a/src/components/aerosol_plus_gas_cycles.jl +++ b/src/components/aerosol_plus_gas_cycles.jl @@ -4,31 +4,31 @@ @defcomp aerosol_plus_cycles begin - aerosol_plus_gases = Index() # Index for tropospheric ozone precursors, aerosols and reactive gas - - aerosol_plus_0 = Parameter(index=[aerosol_plus_gases]) # Aerosol+ gas concentration in initial model period (ppb). - aerosol_plus_pi = Parameter(index=[aerosol_plus_gases]) # Pre-industrial Aerosol+ concentrations (ppb). - g0_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 Aerosol+ gas cycle parameterization at α=1. - g1_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 Aerosol+ gas cycle parameterization at α=1. - emiss2conc_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Conversion factor between emissions and concentrations. - GU_aerosol_plus_0 = Parameter(index=[aerosol_plus_gases]) # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Strength of pre-industrial uptake from atmosphere. - rA_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). - rT_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). - rU_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Fraction of emissions entering iᵗʰ Aerosol+ gas reservior. - R0_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + aerosol_plus_gases = Index() # Index for tropospheric ozone precursors, aerosols and reactive gas + + aerosol_plus_0 = Parameter(index=[aerosol_plus_gases]) # Aerosol+ gas concentration in initial model period (ppb). + aerosol_plus_pi = Parameter(index=[aerosol_plus_gases]) # Pre-industrial Aerosol+ concentrations (ppb). + g0_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 Aerosol+ gas cycle parameterization at α=1. + g1_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 Aerosol+ gas cycle parameterization at α=1. + emiss2conc_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Conversion factor between emissions and concentrations. + GU_aerosol_plus_0 = Parameter(index=[aerosol_plus_gases]) # Initial model period value for cumulative uptake of agent (unit of E⁻¹). + r0_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Strength of pre-industrial uptake from atmosphere. + rA_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). + rT_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). + rU_aerosol_plus = Parameter(index=[aerosol_plus_gases]) # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). + τ_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Fraction of emissions entering iᵗʰ Aerosol+ gas reservior. + R0_aerosol_plus = Parameter(index=[aerosol_plus_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). E_aerosol_plus = Parameter(index=[time, aerosol_plus_gases]) # Annual Aerosol+ emissions (Tg yr⁻¹). - Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). + Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). α_aerosol_plus = Variable(index=[time, aerosol_plus_gases]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. aerosol_plus_conc = Variable(index=[time, aerosol_plus_gases]) # Total atmospheric Aerosol+ concentrations (ppb). GA_aerosol_plus = Variable(index=[time, aerosol_plus_gases]) # Atmospheric burden of agent above pre-industrial levels (unit of E). GU_aerosol_plus = Variable(index=[time, aerosol_plus_gases]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). iIRFT100_aerosol_plus = Variable(index=[time, aerosol_plus_gases]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_aerosol_plus = Variable(index=[time, aerosol_plus_gases, 4]) - R_aerosol_plus = Variable(index=[time, aerosol_plus_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + decay_rate_aerosol_plus = Variable(index=[time, aerosol_plus_gases, 4]) # Decay rates + R_aerosol_plus = Variable(index=[time, aerosol_plus_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/ch4_cycle.jl b/src/components/ch4_cycle.jl index cae4ae6..83fcd12 100644 --- a/src/components/ch4_cycle.jl +++ b/src/components/ch4_cycle.jl @@ -4,29 +4,29 @@ @defcomp ch4_cycle begin - ch4_0 = Parameter() # Methane concentration in initial model period (ppb). - ch4_pi = Parameter() # Pre-industrial methane concentrations (ppb). - g0_ch4 = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 methane cycle parameterization at α=1. - g1_ch4 = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 methane cycle parameterization at α=1. - emiss2conc_ch4 = Parameter() # Conversion factor between emissions and concentrations. - GU_ch4_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_ch4 = Parameter() # Strength of pre-industrial uptake from atmosphere. - rA_ch4 = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). - rT_ch4 = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). - rU_ch4 = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_ch4 = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_ch4 = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ methane reservior. - R0_ch4 = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). - E_ch4 = Parameter(index=[time]) # Annual methane emissions (TgCH₄ yr⁻¹). - Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). - - α_ch4 = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. - ch4 = Variable(index=[time]) # Total atmospheric methane concentrations (ppb). - GA_ch4 = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). - GU_ch4 = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). - iIRFT100_ch4 = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_ch4 = Variable(index=[time,4]) - R_ch4 = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + ch4_0 = Parameter() # Methane concentration in initial model period (ppb). + ch4_pi = Parameter() # Pre-industrial methane concentrations (ppb). + g0_ch4 = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 methane cycle parameterization at α=1. + g1_ch4 = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 methane cycle parameterization at α=1. + emiss2conc_ch4 = Parameter() # Conversion factor between emissions and concentrations. + GU_ch4_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). + r0_ch4 = Parameter() # Strength of pre-industrial uptake from atmosphere. + rA_ch4 = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). + rT_ch4 = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). + rU_ch4 = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). + τ_ch4 = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_ch4 = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ methane reservior. + R0_ch4 = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + E_ch4 = Parameter(index=[time]) # Annual methane emissions (TgCH₄ yr⁻¹). + Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). + + α_ch4 = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. + ch4 = Variable(index=[time]) # Total atmospheric methane concentrations (ppb). + GA_ch4 = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). + GU_ch4 = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). + iIRFT100_ch4 = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). + decay_rate_ch4 = Variable(index=[time,4]) # Decay rates + R_ch4 = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/co2_cycle.jl b/src/components/co2_cycle.jl index 477ecda..7416d47 100644 --- a/src/components/co2_cycle.jl +++ b/src/components/co2_cycle.jl @@ -4,29 +4,29 @@ @defcomp co2_cycle begin - co2_0 = Parameter() # Carbon dioxide concentration in initial model period (ppm). - co2_pi = Parameter() # Pre-industrial carbon dioxide concentrations (ppm). - g0_co2 = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 carbon cycle parameterization at α=1. - g1_co2 = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 carbon cycle parameterization at α=1. - emiss2conc_co2 = Parameter() # Conversion factor between emissions (GtC) and concentrations (ppm). - GU_co2_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_co2 = Parameter() # Strength of pre-industrial uptake from atmosphere. - rA_co2 = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). - rT_co2 = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). - rU_co2 = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_co2 = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_co2 = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ carbon reservior. - R0_co2 = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). - E_co2 = Parameter(index=[time]) # Annual carbon dioxide emissions (GtC yr⁻¹). - Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). - - α_co2 = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. - co2 = Variable(index=[time]) # Total atmospheric carbon dioxide concentrations (ppm). - GA_co2 = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). - GU_co2 = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). - iIRFT100_co2 = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_co2 = Variable(index=[time,4]) - R_co2 = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + co2_0 = Parameter() # Carbon dioxide concentration in initial model period (ppm). + co2_pi = Parameter() # Pre-industrial carbon dioxide concentrations (ppm). + g0_co2 = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 carbon cycle parameterization at α=1. + g1_co2 = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 carbon cycle parameterization at α=1. + emiss2conc_co2 = Parameter() # Conversion factor between emissions (GtC) and concentrations (ppm). + GU_co2_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). + r0_co2 = Parameter() # Strength of pre-industrial uptake from atmosphere. + rA_co2 = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). + rT_co2 = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). + rU_co2 = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). + τ_co2 = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_co2 = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ carbon reservior. + R0_co2 = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + E_co2 = Parameter(index=[time]) # Annual carbon dioxide emissions (GtC yr⁻¹). + Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). + + α_co2 = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. + co2 = Variable(index=[time]) # Total atmospheric carbon dioxide concentrations (ppm). + GA_co2 = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). + GU_co2 = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). + iIRFT100_co2 = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). + decay_rate_co2 = Variable(index=[time,4]) # Decay rates + R_co2 = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/flourinated_gas_cycles.jl b/src/components/flourinated_gas_cycles.jl index db93aa7..1769941 100644 --- a/src/components/flourinated_gas_cycles.jl +++ b/src/components/flourinated_gas_cycles.jl @@ -4,31 +4,31 @@ @defcomp flourinated_cycles begin - flourinated_gases = Index() # Index for flourinated gases controlled under the Kyoto Protocol. - - flourinated_0 = Parameter(index=[flourinated_gases]) # flourinated gas concentration in initial model period (ppb). - flourinated_pi = Parameter(index=[flourinated_gases]) # Pre-industrial flourinated gas concentrations (ppb). - g0_flourinated = Parameter(index=[flourinated_gases]) # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 flourinated gas cycle parameterization at α=1. - g1_flourinated = Parameter(index=[flourinated_gases]) # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 flourinated gas cycle parameterization at α=1. - emiss2conc_flourinated = Parameter(index=[flourinated_gases]) # Conversion factor between emissions and concentrations. - GU_flourinated_0 = Parameter(index=[flourinated_gases]) # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_flourinated = Parameter(index=[flourinated_gases]) # Strength of pre-industrial uptake from atmosphere. - rA_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). - rT_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). - rU_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_flourinated = Parameter(index=[flourinated_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_flourinated = Parameter(index=[flourinated_gases, 4]) # Fraction of emissions entering iᵗʰ flourinated gas reservior. - R0_flourinated = Parameter(index=[flourinated_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). - E_flourinated = Parameter(index=[time, flourinated_gases]) # Annual flourinated gas emissions (Tg yr⁻¹). - Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). - - α_flourinated = Variable(index=[time, flourinated_gases]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. - flourinated_conc = Variable(index=[time, flourinated_gases]) # Total atmospheric flourinated gas concentrations (ppb). - GA_flourinated = Variable(index=[time, flourinated_gases]) # Atmospheric burden of agent above pre-industrial levels (unit of E). - GU_flourinated = Variable(index=[time, flourinated_gases]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). - iIRFT100_flourinated = Variable(index=[time, flourinated_gases]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_flourinated = Variable(index=[time, flourinated_gases, 4]) - R_flourinated = Variable(index=[time, flourinated_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + flourinated_gases = Index() # Index for flourinated gases controlled under the Kyoto Protocol. + + flourinated_0 = Parameter(index=[flourinated_gases]) # flourinated gas concentration in initial model period (ppb). + flourinated_pi = Parameter(index=[flourinated_gases]) # Pre-industrial flourinated gas concentrations (ppb). + g0_flourinated = Parameter(index=[flourinated_gases]) # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 flourinated gas cycle parameterization at α=1. + g1_flourinated = Parameter(index=[flourinated_gases]) # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 flourinated gas cycle parameterization at α=1. + emiss2conc_flourinated = Parameter(index=[flourinated_gases]) # Conversion factor between emissions and concentrations. + GU_flourinated_0 = Parameter(index=[flourinated_gases]) # Initial model period value for cumulative uptake of agent (unit of E⁻¹). + r0_flourinated = Parameter(index=[flourinated_gases]) # Strength of pre-industrial uptake from atmosphere. + rA_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). + rT_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). + rU_flourinated = Parameter(index=[flourinated_gases]) # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). + τ_flourinated = Parameter(index=[flourinated_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_flourinated = Parameter(index=[flourinated_gases, 4]) # Fraction of emissions entering iᵗʰ flourinated gas reservior. + R0_flourinated = Parameter(index=[flourinated_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + E_flourinated = Parameter(index=[time, flourinated_gases]) # Annual flourinated gas emissions (Tg yr⁻¹). + Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). + + α_flourinated = Variable(index=[time, flourinated_gases]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. + flourinated_conc = Variable(index=[time, flourinated_gases]) # Total atmospheric flourinated gas concentrations (ppb). + GA_flourinated = Variable(index=[time, flourinated_gases]) # Atmospheric burden of agent above pre-industrial levels (unit of E). + GU_flourinated = Variable(index=[time, flourinated_gases]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). + iIRFT100_flourinated = Variable(index=[time, flourinated_gases]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). + decay_rate_flourinated = Variable(index=[time, flourinated_gases, 4]) # Decay rates + R_flourinated = Variable(index=[time, flourinated_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/montreal_gas_cycles.jl b/src/components/montreal_gas_cycles.jl index c555190..d0ab6af 100644 --- a/src/components/montreal_gas_cycles.jl +++ b/src/components/montreal_gas_cycles.jl @@ -4,31 +4,31 @@ @defcomp montreal_cycles begin - montreal_gases = Index() # Index for Ozone Depleting Substances controlled under the Montreal Protocol. + montreal_gases = Index() # Index for Ozone Depleting Substances controlled under the Montreal Protocol. montreal_0 = Parameter(index=[montreal_gases]) # Montreal gas concentration in initial model period (ppb). - montreal_pi = Parameter(index=[montreal_gases]) # Pre-industrial Montreal gas concentrations (ppb). + montreal_pi = Parameter(index=[montreal_gases]) # Pre-industrial Montreal gas concentrations (ppb). g0_montreal = Parameter(index=[montreal_gases]) # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 Montreal gas cycle parameterization at α=1. g1_montreal = Parameter(index=[montreal_gases]) # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 Montreal gas cycle parameterization at α=1. - emiss2conc_montreal = Parameter(index=[montreal_gases]) # Conversion factor between emissions and concentrations. + emiss2conc_montreal = Parameter(index=[montreal_gases]) # Conversion factor between emissions and concentrations. GU_montreal_0 = Parameter(index=[montreal_gases]) # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_montreal = Parameter(index=[montreal_gases]) # Strength of pre-industrial uptake from atmosphere. + r0_montreal = Parameter(index=[montreal_gases]) # Strength of pre-industrial uptake from atmosphere. rA_montreal = Parameter(index=[montreal_gases]) # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). rT_montreal = Parameter(index=[montreal_gases]) # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). rU_montreal = Parameter(index=[montreal_gases]) # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_montreal = Parameter(index=[montreal_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_montreal = Parameter(index=[montreal_gases, 4]) # Fraction of emissions entering iᵗʰ Montreal gas reservior. - R0_montreal = Parameter(index=[montreal_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). - E_montreal = Parameter(index=[time, montreal_gases]) # Annual Montreal gas emissions (Tg yr⁻¹). - Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). + τ_montreal = Parameter(index=[montreal_gases, 4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_montreal = Parameter(index=[montreal_gases, 4]) # Fraction of emissions entering iᵗʰ Montreal gas reservior. + R0_montreal = Parameter(index=[montreal_gases, 4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + E_montreal = Parameter(index=[time, montreal_gases]) # Annual Montreal gas emissions (Tg yr⁻¹). + Tj = Parameter(index=[time, 3]) # Temperature change for three thermal pools (K). α_montreal = Variable(index=[time, montreal_gases]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. montreal_conc = Variable(index=[time, montreal_gases]) # Total atmospheric Montreal gas concentrations (ppb). GA_montreal = Variable(index=[time, montreal_gases]) # Atmospheric burden of agent above pre-industrial levels (unit of E). GU_montreal = Variable(index=[time, montreal_gases]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). iIRFT100_montreal = Variable(index=[time, montreal_gases]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_montreal = Variable(index=[time, montreal_gases, 4]) - R_montreal = Variable(index=[time, montreal_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + decay_rate_montreal = Variable(index=[time, montreal_gases, 4]) # Decay rates + R_montreal = Variable(index=[time, montreal_gases, 4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/n2o_cycle.jl b/src/components/n2o_cycle.jl index 28a29c7..60b2547 100644 --- a/src/components/n2o_cycle.jl +++ b/src/components/n2o_cycle.jl @@ -4,29 +4,29 @@ @defcomp n2o_cycle begin - n2o_0 = Parameter() # Nitrous oxide concentration in initial model period (ppb). - n2o_pi = Parameter() # Pre-industrial nitrous oxide concentrations (ppb). - g0_n2o = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 nitrous oxide cycle parameterization at α=1. - g1_n2o = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 nitrous oxide cycle parameterization at α=1. - emiss2conc_n2o = Parameter() # Conversion factor between emissions and concentrations. - GU_n2o_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). - r0_n2o = Parameter() # Strength of pre-industrial uptake from atmosphere. - rA_n2o = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). - rT_n2o = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). - rU_n2o = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). - τ_n2o = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). - a_n2o = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ nitrous oxide reservior. - R0_n2o = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). - E_n2o = Parameter(index=[time]) # Annual nitrous oxide emissions (TgN yr⁻¹). - Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). - - α_n2o = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. - n2o = Variable(index=[time]) # Total atmospheric nitrous oxide concentrations (ppb). - GA_n2o = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). - GU_n2o = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). - iIRFT100_n2o = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). - decay_rate_n2o = Variable(index=[time,4]) - R_n2o = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). + n2o_0 = Parameter() # Nitrous oxide concentration in initial model period (ppb). + n2o_pi = Parameter() # Pre-industrial nitrous oxide concentrations (ppb). + g0_n2o = Parameter() # Constant to set approximation of value for α equal to Millar et al. (2017) numerical solution for iIRF100 nitrous oxide cycle parameterization at α=1. + g1_n2o = Parameter() # Constant to set approximation of gradient for α equal to Millar et al. (2017) numerical solution for iIRF100 nitrous oxide cycle parameterization at α=1. + emiss2conc_n2o = Parameter() # Conversion factor between emissions and concentrations. + GU_n2o_0 = Parameter() # Initial model period value for cumulative uptake of agent (unit of E⁻¹). + r0_n2o = Parameter() # Strength of pre-industrial uptake from atmosphere. + rA_n2o = Parameter() # Sensitivity of uptake from atmosphere to current atmospheric burden of agent (unit of E⁻¹). + rT_n2o = Parameter() # Sensitivity of uptake from atmosphere to model temperature change since initialization (K⁻¹). + rU_n2o = Parameter() # Sensitivity of uptake from atmosphere to cumulative uptake of agent since model initialization (unit of E⁻¹). + τ_n2o = Parameter(index=[4]) # Atmospheric lifetime of gas in iᵗʰ reservior (years). + a_n2o = Parameter(index=[4]) # Fraction of emissions entering iᵗʰ nitrous oxide reservior. + R0_n2o = Parameter(index=[4]) # Initial model period value for quantity of agent in iᵗʰ atmospheric reservior (unit of E). + E_n2o = Parameter(index=[time]) # Annual nitrous oxide emissions (TgN yr⁻¹). + Tj = Parameter(index=[time,3]) # Temperature change for three thermal pools (K). + + α_n2o = Variable(index=[time]) # State-dependent multiplicative adjustment coefficient of reservior lifetimes. + n2o = Variable(index=[time]) # Total atmospheric nitrous oxide concentrations (ppb). + GA_n2o = Variable(index=[time]) # Atmospheric burden of agent above pre-industrial levels (unit of E). + GU_n2o = Variable(index=[time]) # Cumulative uptake of agent since model initialization (unit of E⁻¹). + iIRFT100_n2o = Variable(index=[time]) # 100-year integrated impulse response function (the average airborne fraction over a 100-year period). + decay_rate_n2o = Variable(index=[time,4]) # Decay rates + R_n2o = Variable(index=[time,4]) # Quantity of agent in iᵗʰ atmospheric reservior (unit of E). function run_timestep(p, v, d, t) diff --git a/src/components/radiative_forcing.jl b/src/components/radiative_forcing.jl index 676f704..5ee7755 100644 --- a/src/components/radiative_forcing.jl +++ b/src/components/radiative_forcing.jl @@ -4,64 +4,80 @@ @defcomp radiative_forcing begin - montreal_gases = Index() # Index for Ozone Depleting Substances controlled under the Montreal Protocol. - aerosol_plus_gases = Index() # Index for tropospheric ozone precursors, aerosols and reactive gas - flourinated_gases = Index() # Index for flourinated gases controlled under the Kyoto Protocol. - - # CO2, CH4, and N2O parametes - co2_conc = Parameter(index=[time]) - ch4_conc = Parameter(index=[time]) - n2o_conc = Parameter(index=[time]) - co2_pi = Parameter() - ch4_pi = Parameter() - n2o_pi = Parameter() - co2_f = Parameter(index=[3]) # f coefficient for - ch4_f = Parameter(index=[3]) # f coefficient for - ch4_o3_f = Parameter(index=[3]) # f coefficient for - ch4_h2o_f = Parameter(index=[3]) # f coefficient for - n2o_f = Parameter(index=[3]) # f coefficient for - n2o_o3_f = Parameter(index=[3]) # f coefficient for + montreal_gases = Index() # Index for Ozone Depleting Substances controlled under the Montreal Protocol. + aerosol_plus_gases = Index() # Index for tropospheric ozone precursors, aerosols and reactive gas + flourinated_gases = Index() # Index for flourinated gases controlled under the Kyoto Protocol. + + ## + ## Parameters + ## + + # CO2, CH4, and N2O parameters + co2_conc = Parameter(index=[time]) + ch4_conc = Parameter(index=[time]) + n2o_conc = Parameter(index=[time]) + co2_pi = Parameter() + ch4_pi = Parameter() + n2o_pi = Parameter() + co2_f = Parameter(index=[3]) # f coefficient for co2 + ch4_f = Parameter(index=[3]) # f coefficient for ch4 + ch4_o3_f = Parameter(index=[3]) # f coefficient for ch4_o3 + ch4_h2o_f = Parameter(index=[3]) # f coefficient for ch4_h2o + n2o_f = Parameter(index=[3]) # f coefficient for n2o_f + n2o_o3_f = Parameter(index=[3]) # f coefficient for n2o_o3 # Aerosol+ parameters - aerosol_plus_conc = Parameter(index=[time, aerosol_plus_gases]) - aerosol_plus_pi = Parameter(index=[aerosol_plus_gases]) - bc_f = Parameter(index=[3]) # f coefficient for - bc_snow_f = Parameter(index=[3]) # f coefficient for - bc_aci_f = Parameter(index=[3]) # f coefficient for - co_f = Parameter(index=[3]) # f coefficient for - co_o3_f = Parameter(index=[3]) # f coefficient for - nh3_f = Parameter(index=[3]) # f coefficient for - nmvoc_f = Parameter(index=[3]) # f coefficient for - nmvoc_o3_f = Parameter(index=[3]) # f coefficient for - nox_f = Parameter(index=[3]) # f coefficient for - nox_o3_f = Parameter(index=[3]) # f coefficient for - nox_avi_f = Parameter(index=[3]) # f coefficient for - nox_avi_contrails_f = Parameter(index=[3]) # f coefficient for - oc_f = Parameter(index=[3]) # f coefficient for - oc_aci_f = Parameter(index=[3]) # f coefficient for - so2_f = Parameter(index=[3]) # f coefficient for - so2_aci_f = Parameter(index=[3]) # f coefficient for + aerosol_plus_conc = Parameter(index=[time, aerosol_plus_gases]) + aerosol_plus_pi = Parameter(index=[aerosol_plus_gases]) + bc_f = Parameter(index=[3]) # f coefficient for bc + bc_snow_f = Parameter(index=[3]) # f coefficient for bc_snow + bc_aci_f = Parameter(index=[3]) # f coefficient for bc_aci + co_f = Parameter(index=[3]) # f coefficient for co + co_o3_f = Parameter(index=[3]) # f coefficient for co_o3 + nh3_f = Parameter(index=[3]) # f coefficient for nh3 + nmvoc_f = Parameter(index=[3]) # f coefficient for nmvoc + nmvoc_o3_f = Parameter(index=[3]) # f coefficient for nmvoc_o3 + nox_f = Parameter(index=[3]) # f coefficient for nox + nox_o3_f = Parameter(index=[3]) # f coefficient for nox_o3 + nox_avi_f = Parameter(index=[3]) # f coefficient for nox_avi + nox_avi_contrails_f = Parameter(index=[3]) # f coefficient for nox_avi_contrails + oc_f = Parameter(index=[3]) # f coefficient for oc + oc_aci_f = Parameter(index=[3]) # f coefficient for oc_aci + so2_f = Parameter(index=[3]) # f coefficient for so2 + so2_aci_f = Parameter(index=[3]) # f coefficient for so2_aci # Montreal gas parameters - montreal_conc = Parameter(index=[time, montreal_gases]) - montreal_pi = Parameter(index=[montreal_gases]) - montreal_f = Parameter(index=[montreal_gases, 3]) # f coefficient - montreal_ind_f = Parameter(index=[montreal_gases, 3]) # f coefficient for indirect forcing effect on O3 + montreal_conc = Parameter(index=[time, montreal_gases]) + montreal_pi = Parameter(index=[montreal_gases]) + montreal_f = Parameter(index=[montreal_gases, 3]) # f coefficient + montreal_ind_f = Parameter(index=[montreal_gases, 3]) # f coefficient for indirect forcing effect on O3 # Flourinated gas parameters - flourinated_conc = Parameter(index=[time, flourinated_gases]) - flourinated_pi = Parameter(index=[flourinated_gases]) - flourinated_f = Parameter(index=[flourinated_gases, 3]) # f coefficient + flourinated_conc = Parameter(index=[time, flourinated_gases]) + flourinated_pi = Parameter(index=[flourinated_gases]) + flourinated_f = Parameter(index=[flourinated_gases, 3]) # f coefficient - exogenous_RF = Parameter(index=[time]) + # Other parameters + exogenous_RF = Parameter(index=[time]) + ## + ## Variables + ## + + # CO2, CH4, and N2O variables + co2_RF = Variable(index=[time]) + ch4_RF = Variable(index=[time]) + ch4_o3_RF = Variable(index=[time]) + ch4_h2o_RF = Variable(index=[time]) + n2o_RF = Variable(index=[time]) + n2o_o3_RF = Variable(index=[time]) + + # Aerosol+ variables bc_RF = Variable(index=[time]) bc_snow_RF = Variable(index=[time]) bc_aci_RF = Variable(index=[time]) - co_RF = Variable(index=[time]) co_o3_RF = Variable(index=[time]) - nh3_RF = Variable(index=[time]) nmvoc_RF = Variable(index=[time]) nmvoc_o3_RF = Variable(index=[time]) @@ -74,20 +90,14 @@ so2_RF = Variable(index=[time]) so2_aci_RF = Variable(index=[time]) - co2_RF = Variable(index=[time]) - - ch4_RF = Variable(index=[time]) - ch4_o3_RF = Variable(index=[time]) - ch4_h2o_RF = Variable(index=[time]) - - n2o_RF = Variable(index=[time]) - n2o_o3_RF = Variable(index=[time]) - + # Montreal gas variables montreal_RF = Variable(index=[time, montreal_gases]) montreal_ind_RF = Variable(index=[time, montreal_gases]) - + + # Flourinated gas variables flourinated_RF = Variable(index=[time, flourinated_gases]) + # Other variables total_RF = Variable(index=[time]) function run_timestep(p, v, d, t) diff --git a/src/components/temperature.jl b/src/components/temperature.jl index 57b3605..208deb3 100644 --- a/src/components/temperature.jl +++ b/src/components/temperature.jl @@ -6,10 +6,10 @@ #d = Parameter(index=[3]) # Thermal response timescales: [1] Thermal equilibration of deep ocean & [2] Thermal admustment of upper ocean (years). decay_factor = Parameter(index=[3]) # Thermal response decay factor, calculated as exp(-1/d) where d represents thermal response timescale of jth thermal box. - q = Parameter(index=[3]) # Raditive forcing coefficient: [1] Thermal equilibration of deep ocean & [2] Thermal admustment of upper ocean (K W⁻¹m²). - F = Parameter(index=[time]) # Total radiative forcing (Wm⁻²). - Tj_0 = Parameter(index=[3]) - T_0 = Parameter() + q = Parameter(index=[3]) # Raditive forcing coefficient: [1] Thermal equilibration of deep ocean & [2] Thermal admustment of upper ocean (K W⁻¹m²). + F = Parameter(index=[time]) # Total radiative forcing (Wm⁻²). + Tj_0 = Parameter(index=[3]) + T_0 = Parameter() T = Variable(index=[time]) # Global mean surface temperature anomaly (K). Tj = Variable(index=[time,3]) # Temperature change for three thermal pools (K). diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 7da2811..41cb731 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -1,6 +1,6 @@ # #---------------------------------------------------------------------------------------------------------------------- # #---------------------------------------------------------------------------------------------------------------------- -# This file contains functions and other snippets of code that are used in various calculations for Mimi-FAIRv2.0. +# This file contains functions and other snippets of code that are used in various calculations for MimiFAIRv2. # #---------------------------------------------------------------------------------------------------------------------- # #---------------------------------------------------------------------------------------------------------------------- @@ -21,6 +21,20 @@ # F2x: Radiative forcing from a doubling of carbon dioxide concentrations (Wm⁻²). #---------------------------------------------------------------------------------------------------------------------- +""" + get_thermal_parameter_defaults(;TCR::Float64=1.79, RWF::Float64=0.552, F2x::Float64=3.759) + +Return a dataframe of the FAIRv2.0.0-alpha default climate parameters using the +optional keyword arguments for transient climate response (K) (TCR; defaults to 1.79), +realized warming fraction in ratio of TCR/ECS (RWF; defaults to 0.552), and radiative +forcing from a double of carbon dioxide concentrations (Wm⁻²) (F2x; defaults to +3.759). + +The response timescales d1-3 (and the shortest-timescale coefficient, q1) are set +to the central estimate of a CMIP6 inferred distribution constrained with observational +warming. The constraint does not significantly affect the central estimates of +the prior (ie. raw CMIP6 inference) distribution. +""" function get_thermal_parameter_defaults(;TCR::Float64=1.79, RWF::Float64=0.552, F2x::Float64=3.759) # Set default values for d1, d2, d3, and q1 parameters. @@ -63,6 +77,15 @@ end #---------------------------------------------------------------------------------------------------------------------- # Function version for a single gas (where a and τ are vectors). +""" + calculate_g0_g1(a::Array{Float64,1}, τ::Array{Float64,1}) + +For a SINGLE gas, calculate and return constants for estimating the state-dependent +adjustment coefficient of a reservior's lifetime (α) such that it approximates +the Millar et al. (2017) numerical solution for a iIRF100 carbon cycle parameterization +at α=1. The two required arguments are a, the fraction of emissions entering iᵗʰ +atmospheric pool, and τ, the atmospheric lifetime of gas in iᵗʰ pool. +""" function calculate_g0_g1(a::Array{Float64,1}, τ::Array{Float64,1}) g1 = sum(a .* τ .* (1.0 .- (1.0 .+ 100.0 ./ τ) .* exp.(-100.0 ./ τ))) g0 = exp(-1.0 * sum(a .* τ .* (1.0 .- exp.(-100.0 ./ τ))) / g1) @@ -70,6 +93,15 @@ function calculate_g0_g1(a::Array{Float64,1}, τ::Array{Float64,1}) end # Function version for multiple gases (where a and τ are 2-d arrays). +""" + calculate_g0_g1(a::Array{Float64,1}, τ::Array{Float64,1}) + +For MULITPLE gases, calculate and return constants for estimating the state-dependent +adjustment coefficient of a reservior's lifetime (α) such that it approximates +the Millar et al. (2017) numerical solution for a iIRF100 carbon cycle parameterization +at α=1. The two required arguments are a, the fraction of emissions entering iᵗʰ +atmospheric pool, and τ, the atmospheric lifetime of gas in iᵗʰ pool. +""" function calculate_g0_g1(a::Array{Float64,2}, τ::Array{Float64,2}) g1 = vec(sum(a .* τ .* (1.0 .- (1.0 .+ 100.0 ./ τ) .* exp.(-100.0 ./ τ)), dims=2)) g0 = vec(exp.(-1.0 .* sum(a .* τ .* (1.0 .- exp.(-100.0 ./ τ)), dims=2) ./ g1)) @@ -91,6 +123,18 @@ end # C_pi: Pre-industrial concentration of forcing agent in atmosphere. #---------------------------------------------------------------------------------------------------------------------- +""" + calculate_RF(f, C, C_pi) + +Calculate and return the effective radiative forcing based on the atmospheric +concentration of a greenhouse gas or other forcing agent following Equation (6) +in Leach et al. (2021). + +The required arguments are as follows: + f: A vector of three concentration-forcing coefficients (1 = logarithmic, 2 = linear, 3 = square root) + C: Concentration of forcing agent in atmosphere + C_pi: Pre-industrial concentration of forcing agent in atmosphere +""" function calculate_RF(f, C, C_pi) if C <= 0.0 diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..da693d4 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,184 @@ +using DataFrames +using CSVFiles +using Test +using MimiFAIRv2 +using Mimi + +using MimiFAIRv2: get_model # load `get_model` function to avoid need for `MimiFAIRv2.` prefix + +@testset "UI" begin + + # run basic steps presented in README + m = get_model() + run(m) + fair_temps = m[:temperature, :T] + + # Run model with altered keyword args and make sure results reflect the + # appropriate changes, or at least a change at all + # - start and end years + # - emissions scenario + # - TCR, RWF, and F2x + + # end year + m = get_model(end_year = 2100) + run(m) + fair_temps = m[:temperature, :T] + @test length(fair_temps) == length(1750:2100) + + # start year - should throw a warning if use a start year other than 1750 + @test_logs (:warn, "Model should not be set to start with a year differing from 1750.") get_model(start_year = 1760) + m = get_model(start_year = 1760) + run(m) + fair_temps = m[:temperature, :T] + @test length(fair_temps) == length(1760:2500) + + # emissions scenario + m1 = get_model(emissions_forcing_scenario = "ssp126") + m2 = get_model(emissions_forcing_scenario = "ssp585") + run(m1) + run(m2) + m1_fair_temps = m1[:temperature, :T] + m2_fair_temps = m2[:temperature, :T] + # Mimi.plot(m1, :temperature, :T) + # Mimi.plot(m2, :temperature, :T) + @test m1_fair_temps !== m2_fair_temps + @test maximum(m1_fair_temps) < maximum(m2_fair_temps) + + # TCR - transient climate response (default to 1.79) + m1 = get_model() + m2 = get_model(TCR = 2.0) + run(m1) + run(m2) + m1_fair_temps = m1[:temperature, :T] + m2_fair_temps = m2[:temperature, :T] + # Mimi.plot(m1, :temperature, :T) + # Mimi.plot(m2, :temperature, :T) + @test m1_fair_temps !== m2_fair_temps + @test maximum(m1_fair_temps) < maximum(m2_fair_temps) + + # RWF - realized warming fraction (default to 0.552) + m1 = get_model() + m2 = get_model(RWF = 0.5) + run(m1) + run(m2) + m1_fair_temps = m1[:temperature, :T] + m2_fair_temps = m2[:temperature, :T] + # Mimi.plot(m1, :temperature, :T) + # Mimi.plot(m2, :temperature, :T) + @test m1_fair_temps !== m2_fair_temps + @test maximum(m1_fair_temps) < maximum(m2_fair_temps) + + # F2x - forcing from a doubling of CO₂ (default to 3.759) + m1 = get_model() + m2 = get_model(F2x = 3.5) + run(m1) + run(m2) + m1_fair_temps = m1[:temperature, :T] + m2_fair_temps = m2[:temperature, :T] + # Mimi.plot(m1, :temperature, :T) + # Mimi.plot(m2, :temperature, :T) + @test m1_fair_temps !== m2_fair_temps + @test maximum(m1_fair_temps) < maximum(m2_fair_temps) + +end + +@testset "Python Comparison - With Python Replication Data" begin + + # run model for each possible SSP and compare variable values between julia + # and python versons + + for ssp in ["ssp119", "ssp126", "ssp245", "ssp370", "ssp585"] + + #Load initial conditions (just need parameter names) and extract gas names. + init_gas_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7)) + montreal_init = filter(:gas_group => ==("montreal"), init_gas_vals) + flourinated_init = filter(:gas_group => ==("flourinated"), init_gas_vals) + aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals) + + # Sort arrays of initial conditions for multiple gases so they are listed alphabetically. + sort!(montreal_init, :gas_name) + sort!(flourinated_init, :gas_name) + sort!(aerosol_plus_init, :gas_name) + + # Load replication emissions and forcing from Python. + python_emiss = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_emissions.csv"))) + python_exog_RF = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_forcing.csv")))[:,:External] + + # Extract emissions arrays for multi-gas groupings. + montreal_emissions = python_emiss[:, Symbol.(montreal_init.gas_name)] + flourinated_emissions = python_emiss[:, Symbol.(flourinated_init.gas_name)] + aerosol_plus_emissions = python_emiss[:, Symbol.(aerosol_plus_init.gas_name)] + + # Get a model instance for that SSP. + m = MimiFAIRv2.get_model(emissions_forcing_scenario = ssp, start_year=1750, end_year=2100) + + # Update exogenous forcing and emissions data. + update_param!(m, :radiative_forcing, :exogenous_RF, python_exog_RF) + update_param!(m, :co2_cycle, :E_co2, python_emiss.carbon_dioxide) + update_param!(m, :ch4_cycle, :E_ch4, python_emiss.methane) + update_param!(m, :n2o_cycle, :E_n2o, python_emiss.nitrous_oxide) + update_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions)) + update_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions)) + update_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions)) + + run(m) + + # Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe. + julia_results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T]) + + # compare temperatures + python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_temperature.csv"))) |> DataFrame + @test maximum(abs, julia_results[!, :temp] .- python_results[!, :default]) ≈ 0. atol = 1e-3 + + # compare radiative forcing + python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_forcing.csv"))) |> DataFrame + select!(python_results, "Forcing component", "Total") + @test maximum(abs, julia_results[!, :rf] .- python_results[!, :Total]) ≈ 0. atol = 1e-3 + + # compare emissions + conv = select!(DataFrame(load(joinpath(@__DIR__, "..", "data", "default_gas_cycle_parameters.csv"), skiplines_begin=6)), :gas_name, :emis2conc) + filter(:gas_name => ==("carbon_dioxide"), conv) + + python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_concentrations.csv"))) |> DataFrame + select!(python_results, "Gas name", "carbon_dioxide") + @test maximum(abs, julia_results[!, :co2] .- python_results[!, :carbon_dioxide]) ≈ 0. atol = 1e-3 + + end +end + +@testset "Python Comparison - With Default Model" begin + + # run model for each possible SSP and compare variable values between julia + # and python versons + + for ssp in ["ssp119", "ssp126", "ssp245", "ssp370", "ssp585"] + + # Get a model instance for that SSP. + m = MimiFAIRv2.get_model(emissions_forcing_scenario = ssp, start_year=1750, end_year=2100) + run(m) + + # Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe. + julia_results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T]) + + python_results = DataFrame(:year => collect(1750:2100)) + + # add temperature + res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_temperature.csv"))) |> DataFrame + insertcols!(python_results, 2, :temp => res[!, :default]) + + # compare radiative forcing + res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_forcing.csv"))) |> DataFrame + insertcols!(python_results, 2, :rf => res[!, :Total]) + + # compare emissions + res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_concentrations.csv"))) |> DataFrame + insertcols!(python_results, 2, :co2 => res[!, :carbon_dioxide]) + + for name in [:co2, :rf, :temp] + maxdiff = maximum(abs, julia_results[!, name] .- python_results[!, name]) + if maxdiff > 1e-3 + @warn("Maximum absolute difference for $(ssp)'s $name is $maxdiff") + end + end + end +end diff --git a/test/test_plot_comparison_code.jl b/test/test_plot_comparison_code.jl new file mode 100644 index 0000000..3e79841 --- /dev/null +++ b/test/test_plot_comparison_code.jl @@ -0,0 +1,152 @@ +using DataFrames +using CSVFiles +using Mimi + +include(joinpath(@__DIR__, "..", "src/MimiFAIRv2.jl")) + +#Load initial conditions (just need parameter names) and extract gas names. +init_gas_vals = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7)) +montreal_init = filter(:gas_group => ==("montreal"), init_gas_vals) +flourinated_init = filter(:gas_group => ==("flourinated"), init_gas_vals) +aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals) + +# Sort arrays of initial conditions for multiple gases so they are listed alphabetically. +sort!(montreal_init, :gas_name) +sort!(flourinated_init, :gas_name) +sort!(aerosol_plus_init, :gas_name) + +# Set SSP (options = ['ssp119','ssp126','ssp245','ssp370','ssp585']) +ssp = "ssp370" + +# Load replication emissions and forcing from Python. +python_emiss = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_emissions.csv"))) +python_exog_RF = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_forcing.csv")))[:,:External] + +# Extract emissions arrays for multi-gas groupings. +montreal_emissions = python_emiss[:, Symbol.(montreal_init.gas_name)] +flourinated_emissions = python_emiss[:, Symbol.(flourinated_init.gas_name)] +aerosol_plus_emissions = python_emiss[:, Symbol.(aerosol_plus_init.gas_name)] + +# Get a model instance for that SSP. +m = MimiFAIRv2.get_model(emissions_forcing_scenario=ssp, start_year=1750, end_year=2100) + +# Update exogenous forcing and emissions data. +update_param!(m, :radiative_forcing, :exogenous_RF, python_exog_RF) +update_param!(m, :co2_cycle, :E_co2, python_emiss.carbon_dioxide) +update_param!(m, :ch4_cycle, :E_ch4, python_emiss.methane) +update_param!(m, :n2o_cycle, :E_n2o, python_emiss.nitrous_oxide) +update_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions)) +update_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions)) +update_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions)) + +run(m) + +# Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe. +results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T]) + +# Save results. +save("Mimi_FAIR_replication_"*ssp*".csv", results) + + + +#------------------------------------------------------ +#------------------------------------------------------ +#------------------------------------------------------ +# PLOTTING CODE IN R +#------------------------------------------------------ +#------------------------------------------------------ +#------------------------------------------------------ + +library(ggplot2) + +setwd("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl") + +# Load Mimi results +mimi_ssp119 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp119.csv") +mimi_ssp126 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp126.csv") +mimi_ssp245 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp245.csv") +mimi_ssp370 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp370.csv") +mimi_ssp585 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp585.csv") + +# Load Python temperature results +python_ssp119 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp119_temperature.csv") +python_ssp126 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp126_temperature.csv") +python_ssp245 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp245_temperature.csv") +python_ssp370 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp370_temperature.csv") +python_ssp585 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp585_temperature.csv") + +# Put data into a dataframe. +mimi_data = data.frame(year=1750:2100, ssp119=mimi_ssp119[,3], ssp126=mimi_ssp126[,3], ssp245=mimi_ssp245[,3], ssp370=mimi_ssp370[,3], ssp585=mimi_ssp585[,3]) +python_data = data.frame(year=1750:2100, ssp119=python_ssp119[,2], ssp126=python_ssp126[,2], ssp245=python_ssp245[,2], ssp370=python_ssp370[,2], ssp585=python_ssp585[,2]) + +# Crop it to all start in the year 2000 (then do historical in black). +index_2000 = which(1750:2100 == 2000) + +short_mimi_data = mimi_data[(index_2000+1):351, ] +short_python_data = python_data[(index_2000+1):351, ] + +# Get historicla data +historic_mimi_data = mimi_data[1:(index_2000+1), ] +historic_python_data = python_data[1:(index_2000+1), ] + +# Thin python data to every 5-years and to shortened time frame +thin_indices = seq(1,length(short_python_data[,1]), by=3) +thin_indices_historic = seq(1,length(historic_python_data[,1]), by=1) + +python_data_thinned = short_python_data[thin_indices, ] +historic_python_data_thinned = historic_python_data[thin_indices_historic, ] + +ssp119_col = "#6ace85" +ssp126_col = "#7f9be7" +ssp245_col = "#9c4ed2" +ssp370_col = "darkorange" +#ssp370_col = "#9c4ed2" +ssp585_col = "#ea3675" + + +p = ggplot() + +p = p + geom_hline(yintercept=0, color="red", linetype="22") + +p = p + geom_line(data=historic_mimi_data, aes(x=year, y=ssp119), color="black") +p = p + geom_point(data=historic_python_data_thinned, aes(x=year, y=ssp119), shape=21, size=1) + +p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp119), color=ssp119_col) +p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp119), fill=ssp119_col, shape=21, size=1) + +p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp126), color=ssp126_col) +p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp126), fill=ssp126_col, shape=21, size=1) + +p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp245), color=ssp245_col) +p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp245), fill=ssp245_col, shape=21, size=1) + +p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp370), color=ssp370_col) +p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp370), fill=ssp370_col, shape=21, size=1) + +p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp585), color=ssp585_col) +p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp585), fill=ssp585_col, shape=21, size=1) + +p = p + xlab("Year") +p = p + ylab("Global Surface Temperature Anomaly (K)") +#p = p + labs(title = "Python vs. Julia-Mimi versions of FaIR v2.0") +p = p + labs(title = "Python vs. Julia-Mimi versions of FaIR v2.0", subtitle = "RCMIP Emission & Forcing Scenarios (1750-2100)") + +p = p + theme(panel.background = element_rect(fill = "transparent"), + plot.background = element_rect(fill = "transparent", color="NA"), + panel.grid.minor = element_blank(), + panel.grid.major = element_line(color="gray90", size=0.25), + #axis.line.y = element_blank(), + axis.line = element_line(color="black", size=0.25), + axis.ticks = element_line(color="black", size=0.25), + axis.text = element_text(size=9, colour="black"), + axis.title = element_text(size=9, colour="black"), + legend.position="none", + #plot.title = element_blank(), + #axis.title.y = element_blank(), + axis.ticks.length=unit(.1, "cm")) + #axis.text.y = element_blank(), + #axis.ticks.y = element_blank()) + +ggsave(p, file="Python_Mimi_FAIR2_temperature_comparison.jpg", device="jpeg", type="cairo", width=200, height=130, unit="mm", dpi=200) +ggsave(p, file="Python_Mimi_FAIR2_temperature_comparison.pdf", device="pdf", width=200, height=130, unit="mm", dpi=200, useDingbats=FALSE) +