diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 000000000..3d7aa22e4 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,1622 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.10.10" +manifest_format = "2.0" +project_hash = "fc5090dba0d641efaadaa43b2214980eadcc9eb9" + +[[deps.AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.5.0" + + [deps.AbstractFFTs.extensions] + AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + AbstractFFTsTestExt = "Test" + + [deps.AbstractFFTs.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.Accessors]] +deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"] +git-tree-sha1 = "3b86719127f50670efe356bc11073d84b4ed7a5d" +uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +version = "0.1.42" + + [deps.Accessors.extensions] + AxisKeysExt = "AxisKeys" + IntervalSetsExt = "IntervalSets" + LinearAlgebraExt = "LinearAlgebra" + StaticArraysExt = "StaticArrays" + StructArraysExt = "StructArrays" + TestExt = "Test" + UnitfulExt = "Unitful" + + [deps.Accessors.weakdeps] + AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" + Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[[deps.Adapt]] +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "7e35fca2bdfba44d797c53dfe63a51fabf39bfc0" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "4.4.0" +weakdeps = ["SparseArrays", "StaticArrays"] + + [deps.Adapt.extensions] + AdaptSparseArraysExt = "SparseArrays" + AdaptStaticArraysExt = "StaticArrays" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "d57bd3762d308bded22c3b82d033bff85f6195c6" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.4.0" + +[[deps.ArrayInterface]] +deps = ["Adapt", "LinearAlgebra"] +git-tree-sha1 = "dbd8c3bbbdbb5c2778f85f4422c39960eac65a42" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "7.20.0" + + [deps.ArrayInterface.extensions] + ArrayInterfaceBandedMatricesExt = "BandedMatrices" + ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" + ArrayInterfaceCUDAExt = "CUDA" + ArrayInterfaceCUDSSExt = "CUDSS" + ArrayInterfaceChainRulesCoreExt = "ChainRulesCore" + ArrayInterfaceChainRulesExt = "ChainRules" + ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" + ArrayInterfaceMetalExt = "Metal" + ArrayInterfaceReverseDiffExt = "ReverseDiff" + ArrayInterfaceSparseArraysExt = "SparseArrays" + ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" + ArrayInterfaceTrackerExt = "Tracker" + + [deps.ArrayInterface.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" + ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + Metal = "dde4c033-4e86-420c-a63e-0dd931031962" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Atomix]] +deps = ["UnsafeAtomics"] +git-tree-sha1 = "29bb0eb6f578a587a49da16564705968667f5fa8" +uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458" +version = "1.1.2" + + [deps.Atomix.extensions] + AtomixCUDAExt = "CUDA" + AtomixMetalExt = "Metal" + AtomixOpenCLExt = "OpenCL" + AtomixoneAPIExt = "oneAPI" + + [deps.Atomix.weakdeps] + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + Metal = "dde4c033-4e86-420c-a63e-0dd931031962" + OpenCL = "08131aa3-fb12-5dee-8b74-c09406e224a2" + oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" + +[[deps.BFloat16s]] +deps = ["LinearAlgebra", "Printf", "Random"] +git-tree-sha1 = "3b642331600250f592719140c60cf12372b82d66" +uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b" +version = "0.5.1" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.BitFlags]] +git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d" +uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" +version = "0.1.9" + +[[deps.BitTwiddlingConvenienceFunctions]] +deps = ["Static"] +git-tree-sha1 = "f21cfd4950cb9f0587d5067e69405ad2acd27b87" +uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" +version = "0.1.6" + +[[deps.Blosc_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "535c80f1c0847a4c967ea945fca21becc9de1522" +uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9" +version = "1.21.7+0" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1b96ea4a01afe0ea4090c5c8039690672dd13f2e" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.9+0" + +[[deps.CEnum]] +git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.5.0" + +[[deps.CFTime]] +deps = ["Dates", "Printf"] +git-tree-sha1 = "2ed76cf4ac70526e3df565435d65e7c7b5c7a77a" +uuid = "179af706-886a-5703-950a-314cd64e0468" +version = "0.2.4" + +[[deps.CPUSummary]] +deps = ["CpuId", "IfElse", "PrecompileTools", "Preferences", "Static"] +git-tree-sha1 = "f3a21d7fc84ba618a779d1ed2fcca2e682865bab" +uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" +version = "0.2.7" + +[[deps.CUDA]] +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Compiler_jll", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "GPUToolbox", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics", "demumble_jll"] +git-tree-sha1 = "224c078c40f7cbd98f393d1722953f3a7f24a3ca" +pinned = true +uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" +version = "5.8.5" + + [deps.CUDA.extensions] + ChainRulesCoreExt = "ChainRulesCore" + EnzymeCoreExt = "EnzymeCore" + SparseMatricesCSRExt = "SparseMatricesCSR" + SpecialFunctionsExt = "SpecialFunctions" + + [deps.CUDA.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" + SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[[deps.CUDA_Compiler_jll]] +deps = ["Artifacts", "CUDA_Driver_jll", "CUDA_Runtime_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "8c4f340dd6501a93c4b99b690797772e4a203099" +uuid = "d1e2174e-dfdc-576e-b43e-73b79eb1aca8" +version = "0.2.1+0" + +[[deps.CUDA_Driver_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e6a1d9f5518122c186fd27786b61d2053cfa1b0c" +uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc" +version = "13.0.1+0" + +[[deps.CUDA_Runtime_Discovery]] +deps = ["Libdl"] +git-tree-sha1 = "f9a521f52d236fe49f1028d69e549e7f2644bb72" +uuid = "1af6417a-86b4-443c-805f-a4643ffb695f" +version = "1.0.0" + +[[deps.CUDA_Runtime_jll]] +deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "e24c6de116c0735c37e83b8bc05ed60d4d359693" +uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" +version = "0.19.1+0" + +[[deps.CatViews]] +deps = ["Random", "Test"] +git-tree-sha1 = "23d1f1e10d4e24374112fcf800ac981d14a54b24" +uuid = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31" +version = "1.0.0" + +[[deps.ChunkCodecCore]] +git-tree-sha1 = "51f4c10ee01bda57371e977931de39ee0f0cdb3e" +uuid = "0b6fb165-00bc-4d37-ab8b-79f91016dbe1" +version = "1.0.0" + +[[deps.ChunkCodecLibZlib]] +deps = ["ChunkCodecCore", "Zlib_jll"] +git-tree-sha1 = "cee8104904c53d39eb94fd06cbe60cb5acde7177" +uuid = "4c0bbee4-addc-4d73-81a0-b6caacae83c8" +version = "1.0.0" + +[[deps.ChunkCodecLibZstd]] +deps = ["ChunkCodecCore", "Zstd_jll"] +git-tree-sha1 = "34d9873079e4cb3d0c62926a225136824677073f" +uuid = "55437552-ac27-4d47-9aa3-63184e8fd398" +version = "1.0.0" + +[[deps.ClimaSeaIce]] +deps = ["Adapt", "JLD2", "KernelAbstractions", "Oceananigans", "RootSolvers", "Roots", "SeawaterPolynomials"] +git-tree-sha1 = "bc7928604980aaec717a1cf9029bac8026a9276c" +uuid = "6ba0ff68-24e6-4315-936c-2e99227c95a4" +version = "0.3.9" + +[[deps.CloseOpenIntervals]] +deps = ["Static", "StaticArrayInterface"] +git-tree-sha1 = "05ba0d07cd4fd8b7a39541e31a7b0254704ea581" +uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" +version = "0.1.13" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.8" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "67e11ee83a43eb71ddc950302c53bf33f0690dfe" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.12.1" + + [deps.ColorTypes.extensions] + StyledStringsExt = "StyledStrings" + + [deps.ColorTypes.weakdeps] + StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b" + +[[deps.ColorVectorSpace]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "Requires", "Statistics", "TensorCore"] +git-tree-sha1 = "8b3b6f87ce8f65a2b4f857528fd8d70086cd72b1" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.11.0" +weakdeps = ["SpecialFunctions"] + + [deps.ColorVectorSpace.extensions] + SpecialFunctionsExt = "SpecialFunctions" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "37ea44092930b1811e666c3bc38065d7d87fcc74" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.13.1" + +[[deps.CommonDataModel]] +deps = ["CFTime", "DataStructures", "Dates", "DiskArrays", "Preferences", "Printf", "Statistics"] +git-tree-sha1 = "675149c3c06350dabb9a807ca3dd473de8173703" +uuid = "1fbeeb36-5f17-413c-809b-666fb144f157" +version = "0.4.1" + +[[deps.CommonSolve]] +git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" +uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +version = "0.2.4" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools"] +git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.1" + +[[deps.CommonWorldInvalidations]] +git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0" +uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8" +version = "1.0.0" + +[[deps.Compat]] +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.18.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.1.1+0" + +[[deps.CompositionsBase]] +git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" +uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b" +version = "0.1.2" +weakdeps = ["InverseFunctions"] + + [deps.CompositionsBase.extensions] + CompositionsBaseInverseFunctionsExt = "InverseFunctions" + +[[deps.ConcurrentUtilities]] +deps = ["Serialization", "Sockets"] +git-tree-sha1 = "d9d26935a0bcffc87d2613ce14c527c99fc543fd" +uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" +version = "2.5.0" + +[[deps.ConstructionBase]] +git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.6.0" + + [deps.ConstructionBase.extensions] + ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseLinearAlgebraExt = "LinearAlgebra" + ConstructionBaseStaticArraysExt = "StaticArrays" + + [deps.ConstructionBase.weakdeps] + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.CpuId]] +deps = ["Markdown"] +git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406" +uuid = "adafc99b-e345-5852-983c-f28acb93d879" +version = "0.3.1" + +[[deps.Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + +[[deps.CubedSphere]] +deps = ["TaylorSeries"] +git-tree-sha1 = "afe9e8c11bf816a6fee878ddfc661e0bd138b747" +uuid = "7445602f-e544-4518-8976-18f8e8ae6cdb" +version = "0.3.2" + +[[deps.CubicSplines]] +deps = ["Random", "Test"] +git-tree-sha1 = "4875023d456ea37c581f406b8b1bc35bea95ae67" +uuid = "9c784101-8907-5a6d-9be6-98f00873c89b" +version = "0.2.1" + +[[deps.DataAPI]] +git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.16.0" + +[[deps.DataDeps]] +deps = ["HTTP", "Libdl", "Reexport", "SHA", "Scratch", "p7zip_jll"] +git-tree-sha1 = "8ae085b71c462c2cb1cfedcb10c3c877ec6cf03f" +uuid = "124859b0-ceae-595e-8997-d05f6a7a8dfe" +version = "0.7.13" + +[[deps.DataFrames]] +deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] +git-tree-sha1 = "c967271c27a95160e30432e011b58f42cd7501b5" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "1.8.0" + +[[deps.DataStructures]] +deps = ["OrderedCollections"] +git-tree-sha1 = "6c72198e6a101cccdd4c9731d3985e904ba26037" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.19.1" + +[[deps.DataValueInterfaces]] +git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" +uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" +version = "1.0.0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.15.1" + +[[deps.DiskArrays]] +deps = ["ConstructionBase", "LRUCache", "Mmap", "OffsetArrays"] +git-tree-sha1 = "cb1635fa91dfccffdaa09f99d457b581575e2a67" +uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3" +version = "0.4.17" + +[[deps.Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.12" + + [deps.Distances.extensions] + DistancesChainRulesCoreExt = "ChainRulesCore" + DistancesSparseArraysExt = "SparseArrays" + + [deps.Distances.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.DocStringExtensions]] +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.5" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.ExceptionUnwrapping]] +deps = ["Test"] +git-tree-sha1 = "d36f682e590a83d63d1c7dbd287573764682d12a" +uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4" +version = "0.1.11" + +[[deps.ExprTools]] +git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.10" + +[[deps.ExpressionExplorer]] +git-tree-sha1 = "4a8c0a9eebf807ac42f0f6de758e60a20be25ffb" +uuid = "21656369-7473-754a-2065-74616d696c43" +version = "1.1.3" + +[[deps.FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "97f08406df914023af55ade2f843c39e99c5d969" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.10.0" + +[[deps.FFTW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6d6219a004b8cf1e0b4dbe27a2860b8e04eba0be" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.11+0" + +[[deps.FileIO]] +deps = ["Pkg", "Requires", "UUIDs"] +git-tree-sha1 = "d60eb76f37d7e5a40cc2e7c36974d864b82dc802" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.17.1" +weakdeps = ["HTTP"] + + [deps.FileIO.extensions] + HTTPExt = "HTTP" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.5" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] +git-tree-sha1 = "dc41303865a16274ecb8450c220021ce1e0cf05f" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "1.2.1" +weakdeps = ["StaticArrays"] + + [deps.ForwardDiff.extensions] + ForwardDiffStaticArraysExt = "StaticArrays" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.GPUArrays]] +deps = ["Adapt", "GPUArraysCore", "KernelAbstractions", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "ScopedValues", "Serialization", "Statistics"] +git-tree-sha1 = "8ddb438e956891a63a5367d7fab61550fc720026" +uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +version = "11.2.6" +weakdeps = ["JLD2"] + + [deps.GPUArrays.extensions] + JLD2Ext = "JLD2" + +[[deps.GPUArraysCore]] +deps = ["Adapt"] +git-tree-sha1 = "83cf05ab16a73219e5f6bd1bdfa9848fa24ac627" +uuid = "46192b85-c4d5-4398-a991-12ede77f4527" +version = "0.2.0" + +[[deps.GPUCompiler]] +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "PrecompileTools", "Preferences", "Scratch", "Serialization", "TOML", "Tracy", "UUIDs"] +git-tree-sha1 = "9a8b92a457f55165923fcfe48997b7b93b712fca" +uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" +version = "1.7.2" + +[[deps.GPUToolbox]] +deps = ["LLVM"] +git-tree-sha1 = "5bfe837129bf49e2e049b4f1517546055cc16a93" +uuid = "096a3bc2-3ced-46d0-87f4-dd12716f4bfc" +version = "0.3.0" + +[[deps.Glob]] +git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496" +uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" +version = "1.3.1" + +[[deps.Graphs]] +deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "7a98c6502f4632dbe9fb1973a4244eaa3324e84d" +uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" +version = "1.13.1" + +[[deps.HDF5_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] +git-tree-sha1 = "e94f84da9af7ce9c6be049e9067e511e17ff89ec" +uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" +version = "1.14.6+0" + +[[deps.HTTP]] +deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] +git-tree-sha1 = "5e6fe50ae7f23d171f44e311c2960294aaa0beb5" +uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" +version = "1.10.19" + +[[deps.HashArrayMappedTries]] +git-tree-sha1 = "2eaa69a7cab70a52b9687c8bf950a5a93ec895ae" +uuid = "076d061b-32b6-4027-95e0-9a2c6f6d7e74" +version = "0.2.0" + +[[deps.HostCPUFeatures]] +deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] +git-tree-sha1 = "8e070b599339d622e9a081d17230d74a5c473293" +uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0" +version = "0.1.17" + +[[deps.Hwloc_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "XML2_jll", "Xorg_libpciaccess_jll"] +git-tree-sha1 = "3d468106a05408f9f7b6f161d9e7715159af247b" +uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" +version = "2.12.2+0" + +[[deps.IfElse]] +git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" +uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +version = "0.1.1" + +[[deps.ImageCore]] +deps = ["ColorVectorSpace", "Colors", "FixedPointNumbers", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "PrecompileTools", "Reexport"] +git-tree-sha1 = "8c193230235bbcee22c8066b0374f63b5683c2d3" +uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" +version = "0.10.5" + +[[deps.ImageMorphology]] +deps = ["DataStructures", "ImageCore", "LinearAlgebra", "LoopVectorization", "OffsetArrays", "Requires", "TiledIteration"] +git-tree-sha1 = "895205d762ae24a01689f8cc7ad584b55f1fd005" +uuid = "787d08f9-d448-5407-9aad-5290dd7ab264" +version = "0.4.7" + +[[deps.Inflate]] +git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.5" + +[[deps.InlineStrings]] +git-tree-sha1 = "8f3d257792a522b4601c24a577954b0a8cd7334d" +uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" +version = "1.4.5" + + [deps.InlineStrings.extensions] + ArrowTypesExt = "ArrowTypes" + ParsersExt = "Parsers" + + [deps.InlineStrings.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" + Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" + +[[deps.IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] +git-tree-sha1 = "ec1debd61c300961f98064cfb21287613ad7f303" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2025.2.0+0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.InverseFunctions]] +git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.17" +weakdeps = ["Dates", "Test"] + + [deps.InverseFunctions.extensions] + InverseFunctionsDatesExt = "Dates" + InverseFunctionsTestExt = "Test" + +[[deps.InvertedIndices]] +git-tree-sha1 = "6da3c4316095de0f5ee2ebd875df8721e7e0bdbe" +uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +version = "1.3.1" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "e2222959fbc6c19554dc15174c81bf7bf3aa691c" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.4" + +[[deps.IteratorInterfaceExtensions]] +git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" +uuid = "82899510-4779-5014-852e-03e436cf321d" +version = "1.0.0" + +[[deps.JLD2]] +deps = ["ChunkCodecLibZlib", "ChunkCodecLibZstd", "FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "ScopedValues"] +git-tree-sha1 = "da2e9b4d1abbebdcca0aa68afa0aa272102baad7" +uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +version = "0.6.2" +weakdeps = ["UnPack"] + + [deps.JLD2.extensions] + UnPackExt = "UnPack" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.7.1" + +[[deps.JuliaNVTXCallbacks_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "af433a10f3942e882d3c671aacb203e006a5808f" +uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e" +version = "0.2.1+0" + +[[deps.KernelAbstractions]] +deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"] +git-tree-sha1 = "83c617e9e9b02306a7acab79e05ec10253db7c87" +uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +version = "0.9.38" + + [deps.KernelAbstractions.extensions] + EnzymeExt = "EnzymeCore" + LinearAlgebraExt = "LinearAlgebra" + SparseArraysExt = "SparseArrays" + + [deps.KernelAbstractions.weakdeps] + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.Krylov]] +deps = ["LinearAlgebra", "Printf", "SparseArrays"] +git-tree-sha1 = "d1fc961038207e43982851e57ee257adc37be5e8" +uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +version = "0.10.2" + +[[deps.KrylovPreconditioners]] +deps = ["Adapt", "Graphs", "KernelAbstractions", "LinearAlgebra", "Metis", "SparseArrays"] +git-tree-sha1 = "77e0d2f1a250af347261c9aa89f74b0cfc530a71" +uuid = "45d422c2-293f-44ce-8315-2cb988662dec" +version = "0.3.6" + + [deps.KrylovPreconditioners.extensions] + KrylovPreconditionersAMDGPUExt = "AMDGPU" + KrylovPreconditionersCUDAExt = "CUDA" + KrylovPreconditionersOneAPIExt = "oneAPI" + + [deps.KrylovPreconditioners.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" + +[[deps.LLVM]] +deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Unicode"] +git-tree-sha1 = "ce8614210409eaa54ed5968f4b50aa96da7ae543" +uuid = "929cbde3-209d-540e-8aea-75f648917ca0" +version = "9.4.4" +weakdeps = ["BFloat16s"] + + [deps.LLVM.extensions] + BFloat16sExt = "BFloat16s" + +[[deps.LLVMExtra_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] +git-tree-sha1 = "8e76807afb59ebb833e9b131ebf1a8c006510f33" +uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" +version = "0.0.38+0" + +[[deps.LLVMLoopInfo]] +git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea" +uuid = "8b046642-f1f6-4319-8d3c-209ddc03c586" +version = "1.0.0" + +[[deps.LRUCache]] +git-tree-sha1 = "5519b95a490ff5fe629c4a7aa3b3dfc9160498b3" +uuid = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +version = "1.6.2" +weakdeps = ["Serialization"] + + [deps.LRUCache.extensions] + SerializationExt = ["Serialization"] + +[[deps.LaTeXStrings]] +git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.4.0" + +[[deps.LayoutPointers]] +deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] +git-tree-sha1 = "a9eaadb366f5493a5654e843864c13d8b107548c" +uuid = "10f19ff3-798f-405d-979b-55457f8fc047" +version = "0.1.17" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.4" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "8.4.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.11.0+1" + +[[deps.LibTracyClient_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "d2bc4e1034b2d43076b50f0e34ea094c2cb0a717" +uuid = "ad6e5548-8b26-5c9f-8ef3-ef0ad883f3a5" +version = "0.9.1+6" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "be484f5c92fad0bd8acfef35fe017900b0b73809" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.18.0+0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.29" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.LoggingExtras]] +deps = ["Dates", "Logging"] +git-tree-sha1 = "f00544d95982ea270145636c181ceda21c4e2575" +uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" +version = "1.2.0" + +[[deps.LoopVectorization]] +deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] +git-tree-sha1 = "a9fc7883eb9b5f04f46efb9a540833d1fad974b3" +uuid = "bdcacae8-1622-11e9-2a5c-532679323890" +version = "0.12.173" + + [deps.LoopVectorization.extensions] + ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"] + ForwardDiffNNlibExt = ["ForwardDiff", "NNlib"] + SpecialFunctionsExt = "SpecialFunctions" + + [deps.LoopVectorization.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" + SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[[deps.Lz4_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "191686b1ac1ea9c89fc52e996ad15d1d241d1e33" +uuid = "5ced341a-0733-55b8-9ab6-a4889d929147" +version = "1.10.1+0" + +[[deps.METIS_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "2eefa8baa858871ae7770c98c3c2a7e46daba5b4" +uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" +version = "5.1.3+0" + +[[deps.MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] +git-tree-sha1 = "282cadc186e7b2ae0eeadbd7a4dffed4196ae2aa" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2025.2.0+0" + +[[deps.MPI]] +deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"] +git-tree-sha1 = "a61ecf714d71064b766d481ef43c094d4c6e3c52" +uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" +version = "0.20.23" + + [deps.MPI.extensions] + AMDGPUExt = "AMDGPU" + CUDAExt = "CUDA" + + [deps.MPI.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + +[[deps.MPICH_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "d72d0ecc3f76998aac04e446547259b9ae4c265f" +uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" +version = "4.3.1+0" + +[[deps.MPIPreferences]] +deps = ["Libdl", "Preferences"] +git-tree-sha1 = "c105fe467859e7f6e9a852cb15cb4301126fac07" +uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" +version = "0.1.11" + +[[deps.MPItrampoline_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "e214f2a20bdd64c04cd3e4ff62d3c9be7e969a59" +uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" +version = "5.5.4+0" + +[[deps.MacroTools]] +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.16" + +[[deps.ManualMemory]] +git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" +uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667" +version = "0.1.8" + +[[deps.MappedArrays]] +git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e" +uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" +version = "0.4.2" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MbedTLS]] +deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] +git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf" +uuid = "739be429-bea8-5141-9913-cc70e7f3736d" +version = "1.1.9" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.2+1" + +[[deps.MeshArrays]] +deps = ["CatViews", "Dates", "Distributed", "Glob", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SharedArrays", "SparseArrays", "Statistics", "Unitful"] +git-tree-sha1 = "fea0859c1406389b3e6657b5376b577ffffdd03a" +uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683" +version = "0.3.23" + + [deps.MeshArrays.extensions] + MeshArraysDataDepsExt = ["DataDeps"] + MeshArraysGeoJSONExt = ["GeoJSON"] + MeshArraysJLD2Ext = ["JLD2"] + MeshArraysMakieExt = ["Makie"] + MeshArraysProjExt = ["Proj"] + MeshArraysShapefileExt = ["Shapefile"] + + [deps.MeshArrays.weakdeps] + DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" + GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9" + JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" + Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + Proj = "c94c279d-25a6-4763-9509-64d165bea63e" + Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" + +[[deps.Metis]] +deps = ["CEnum", "LinearAlgebra", "METIS_jll", "SparseArrays"] +git-tree-sha1 = "54aca4fd53d39dcd2c3f1bef367b6921e8178628" +uuid = "2679e427-3c69-5b7f-982b-ece356f1e94b" +version = "1.5.0" + + [deps.Metis.extensions] + MetisGraphs = "Graphs" + MetisLightGraphs = "LightGraphs" + MetisSimpleWeightedGraphs = ["SimpleWeightedGraphs", "Graphs"] + + [deps.Metis.weakdeps] + Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" + LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" + SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" + +[[deps.MicrosoftMPI_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "bc95bf4149bf535c09602e3acdf950d9b4376227" +uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf" +version = "10.1.4+3" + +[[deps.Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.2.0" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MosaicViews]] +deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"] +git-tree-sha1 = "7b86a5d4d70a9f5cdf2dacb3cbe6d251d1a61dbe" +uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389" +version = "0.3.4" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2023.1.10" + +[[deps.MuladdMacro]] +git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" +uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +version = "0.2.4" + +[[deps.NCDatasets]] +deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "DiskArrays", "NetCDF_jll", "NetworkOptions", "Printf"] +git-tree-sha1 = "5d5de6613f231d17cecb13a7dcdbb74740134ba8" +uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +version = "0.14.9" + +[[deps.NVTX]] +deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"] +git-tree-sha1 = "6b573a3e66decc7fc747afd1edbf083ff78c813a" +uuid = "5da4648a-3479-48b8-97b9-01cb529c0a1f" +version = "1.0.1" + +[[deps.NVTX_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "af2232f69447494514c25742ba1503ec7e9877fe" +uuid = "e98f9f5b-d649-5603-91fd-7774390e6439" +version = "3.2.2+0" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.1.3" + +[[deps.NearestNeighbors]] +deps = ["Distances", "StaticArrays"] +git-tree-sha1 = "ca7e18198a166a1f3eb92a3650d53d94ed8ca8a1" +uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +version = "0.4.22" + +[[deps.NetCDF_jll]] +deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML", "XML2_jll", "Zlib_jll", "Zstd_jll", "libaec_jll", "libzip_jll"] +git-tree-sha1 = "d574803b6055116af212434460adf654ce98e345" +uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3" +version = "401.900.300+0" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.Oceananigans]] +deps = ["Adapt", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "GPUArrays", "GPUArraysCore", "Glob", "InteractiveUtils", "JLD2", "KernelAbstractions", "Krylov", "KrylovPreconditioners", "LinearAlgebra", "Logging", "MPI", "MuladdMacro", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "ReactantCore", "Rotations", "SeawaterPolynomials", "SparseArrays", "StaticArrays", "Statistics", "StructArrays"] +git-tree-sha1 = "4780b18e5e966c85098ed576174cdf1cc94bf10b" +uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" +version = "0.100.5" + + [deps.Oceananigans.extensions] + OceananigansAMDGPUExt = "AMDGPU" + OceananigansCUDAExt = "CUDA" + OceananigansEnzymeExt = "Enzyme" + OceananigansMakieExt = ["MakieCore", "Makie"] + OceananigansMetalExt = "Metal" + OceananigansNCDatasetsExt = "NCDatasets" + OceananigansOneAPIExt = "oneAPI" + OceananigansReactantExt = ["Reactant", "KernelAbstractions", "ConstructionBase"] + OceananigansXESMFExt = ["XESMF"] + + [deps.Oceananigans.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b" + Metal = "dde4c033-4e86-420c-a63e-0dd931031962" + NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" + Reactant = "3c362404-f566-11ee-1572-e11a4b42c853" + XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" + oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" + +[[deps.OffsetArrays]] +git-tree-sha1 = "117432e406b5c023f665fa73dc26e79ec3630151" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.17.0" +weakdeps = ["Adapt"] + + [deps.OffsetArrays.extensions] + OffsetArraysAdaptExt = "Adapt" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.23+4" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.5+0" + +[[deps.OpenMPI_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"] +git-tree-sha1 = "ec764453819f802fc1e144bfe750c454181bd66d" +uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" +version = "5.0.8+0" + +[[deps.OpenSSL]] +deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] +git-tree-sha1 = "f1a7e086c677df53e064e0fdd2c9d0b0833e3f6e" +uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c" +version = "1.5.0" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f19301ae653233bc88b1810ae908194f07f8db9d" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "3.5.4+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.6+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.1" + +[[deps.PaddedViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f" +uuid = "5432bcbf-9aad-5242-b902-cca2824c8663" +version = "0.5.12" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.10.0" + +[[deps.PkgVersion]] +deps = ["Pkg"] +git-tree-sha1 = "f9501cc0430a26bc3d156ae1b5b0c1b47af4d6da" +uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" +version = "0.3.3" + +[[deps.PolyesterWeave]] +deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] +git-tree-sha1 = "645bed98cd47f72f67316fd42fc47dee771aefcd" +uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" +version = "0.2.2" + +[[deps.PooledArrays]] +deps = ["DataAPI", "Future"] +git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3" +uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" +version = "1.4.3" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "0f27480397253da18fe2c12a4ba4eb9eb208bf3d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.5.0" + +[[deps.PrettyTables]] +deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "1101cd475833706e4d0e7b122218257178f48f34" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.4.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.Quaternions]] +deps = ["LinearAlgebra", "Random", "RealDot"] +git-tree-sha1 = "994cc27cdacca10e68feb291673ec3a76aa2fae9" +uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" +version = "0.7.6" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.Random123]] +deps = ["Random", "RandomNumbers"] +git-tree-sha1 = "dbe5fd0b334694e905cb9fda73cd8554333c46e2" +uuid = "74087812-796a-5b5d-8853-05524746bad3" +version = "1.7.1" + +[[deps.RandomNumbers]] +deps = ["Random"] +git-tree-sha1 = "c6ec94d2aaba1ab2ff983052cf6a606ca5985902" +uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" +version = "1.6.0" + +[[deps.ReactantCore]] +deps = ["ExpressionExplorer", "MacroTools"] +git-tree-sha1 = "f3e31b90afcd152578a6c389eae46dd38b9a4f38" +uuid = "a3311ec8-5e00-46d5-b541-4f83e724a433" +version = "0.1.16" + +[[deps.RealDot]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9" +uuid = "c1ae055f-0cd5-4b69-90a6-9a35b1a98df9" +version = "0.1.0" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.1" + +[[deps.RootSolvers]] +deps = ["ForwardDiff", "Printf"] +git-tree-sha1 = "769388dbf7656e70f6ee250f35bb6cbca8f43203" +uuid = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" +version = "0.4.6" + +[[deps.Roots]] +deps = ["Accessors", "CommonSolve", "Printf"] +git-tree-sha1 = "8a433b1ede5e9be9a7ba5b1cc6698daa8d718f1d" +uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +version = "2.2.10" + + [deps.Roots.extensions] + RootsChainRulesCoreExt = "ChainRulesCore" + RootsForwardDiffExt = "ForwardDiff" + RootsIntervalRootFindingExt = "IntervalRootFinding" + RootsSymPyExt = "SymPy" + RootsSymPyPythonCallExt = "SymPyPythonCall" + RootsUnitfulExt = "Unitful" + + [deps.Roots.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" + SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" + SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[[deps.Rotations]] +deps = ["LinearAlgebra", "Quaternions", "Random", "StaticArrays"] +git-tree-sha1 = "5680a9276685d392c87407df00d57c9924d9f11e" +uuid = "6038ab10-8711-5258-84ad-4b1120ba62dc" +version = "1.7.1" + + [deps.Rotations.extensions] + RotationsRecipesBaseExt = "RecipesBase" + + [deps.Rotations.weakdeps] + RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SIMDTypes]] +git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" +uuid = "94e857df-77ce-4151-89e5-788b33177be4" +version = "0.1.0" + +[[deps.SLEEFPirates]] +deps = ["IfElse", "Static", "VectorizationBase"] +git-tree-sha1 = "456f610ca2fbd1c14f5fcf31c6bfadc55e7d66e0" +uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" +version = "0.6.43" + +[[deps.SciMLPublic]] +git-tree-sha1 = "ed647f161e8b3f2973f24979ec074e8d084f1bee" +uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b" +version = "1.0.0" + +[[deps.ScopedValues]] +deps = ["HashArrayMappedTries", "Logging"] +git-tree-sha1 = "c3b2323466378a2ba15bea4b2f73b081e022f473" +uuid = "7e506255-f358-4e82-b7e4-beb19740aa63" +version = "1.5.0" + +[[deps.Scratch]] +deps = ["Dates"] +git-tree-sha1 = "9b81b8393e50b7d4e6d0a9f14e192294d3b7c109" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.3.0" + +[[deps.SeawaterPolynomials]] +git-tree-sha1 = "e2671e9abe2a2faa51dcecd9d911522931c16012" +uuid = "d496a93d-167e-4197-9f49-d3af4ff8fe40" +version = "0.3.10" + +[[deps.SentinelArrays]] +deps = ["Dates", "Random"] +git-tree-sha1 = "712fb0231ee6f9120e005ccd56297abbc053e7e0" +uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" +version = "1.4.8" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.SimpleBufferStream]] +git-tree-sha1 = "f305871d2f381d21527c770d4788c06c097c9bc1" +uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" +version = "1.2.0" + +[[deps.SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "be8eeac05ec97d379347584fa9fe2f5f76795bcb" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.5" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "64d974c2e6fdf07f8155b5b2ca2ffa9069b608d9" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.2.2" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "f2685b435df2613e25fc10ad8c26dddb8640f547" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.6.1" + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + + [deps.SpecialFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.StackViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "be1cf4eb0ac528d96f5115b4ed80c26a8d8ae621" +uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" +version = "0.1.2" + +[[deps.Static]] +deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"] +git-tree-sha1 = "1e44e7b1dbb5249876d84c32466f8988a6b41bbb" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "1.3.0" + +[[deps.StaticArrayInterface]] +deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"] +git-tree-sha1 = "96381d50f1ce85f2663584c8e886a6ca97e60554" +uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" +version = "1.8.0" +weakdeps = ["OffsetArrays", "StaticArrays"] + + [deps.StaticArrayInterface.extensions] + StaticArrayInterfaceOffsetArraysExt = "OffsetArrays" + StaticArrayInterfaceStaticArraysExt = "StaticArrays" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] +git-tree-sha1 = "b8693004b385c842357406e3af647701fe783f98" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.9.15" + + [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" + StaticArraysStatisticsExt = "Statistics" + + [deps.StaticArrays.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.3" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.10.0" + +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "9d72a13a3f4dd3795a195ac5a44d7d6ff5f552ff" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.7.1" + +[[deps.StringManipulation]] +deps = ["PrecompileTools"] +git-tree-sha1 = "725421ae8e530ec29bcbdddbe91ff8053421d023" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.4.1" + +[[deps.StructArrays]] +deps = ["ConstructionBase", "DataAPI", "Tables"] +git-tree-sha1 = "8ad2e38cbb812e29348719cc63580ec1dfeb9de4" +uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +version = "0.7.1" +weakdeps = ["Adapt", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "SparseArrays", "StaticArrays"] + + [deps.StructArrays.extensions] + StructArraysAdaptExt = "Adapt" + StructArraysGPUArraysCoreExt = ["GPUArraysCore", "KernelAbstractions"] + StructArraysLinearAlgebraExt = "LinearAlgebra" + StructArraysSparseArraysExt = "SparseArrays" + StructArraysStaticArraysExt = "StaticArrays" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.2.1+1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.1" + +[[deps.Tables]] +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] +git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "1.12.1" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.TaylorSeries]] +deps = ["LinearAlgebra", "Markdown", "Requires", "SparseArrays"] +git-tree-sha1 = "abc13c4d3cccd1703335ba03abdaf0dc8ecc97e2" +uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" +version = "0.17.3" + + [deps.TaylorSeries.extensions] + TaylorSeriesIAExt = "IntervalArithmetic" + TaylorSeriesSAExt = "StaticArrays" + + [deps.TaylorSeries.weakdeps] + IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.TensorCore]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6" +uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50" +version = "0.1.1" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.Thermodynamics]] +deps = ["DocStringExtensions", "ForwardDiff", "KernelAbstractions", "Random", "RootSolvers"] +git-tree-sha1 = "bae39665e42c36e2e6bf5b814631915ec7d90232" +uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" +version = "0.15.0" + + [deps.Thermodynamics.extensions] + CreateParametersExt = "ClimaParams" + + [deps.Thermodynamics.weakdeps] + ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" + +[[deps.ThreadingUtilities]] +deps = ["ManualMemory"] +git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85" +uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" +version = "0.5.5" + +[[deps.TiledIteration]] +deps = ["OffsetArrays", "StaticArrayInterface"] +git-tree-sha1 = "1176cc31e867217b06928e2f140c90bd1bc88283" +uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" +version = "0.5.0" + +[[deps.Tracy]] +deps = ["ExprTools", "LibTracyClient_jll", "Libdl"] +git-tree-sha1 = "73e3ff50fd3990874c59fef0f35d10644a1487bc" +uuid = "e689c965-62c8-4b79-b2c5-8359227902fd" +version = "0.1.6" + + [deps.Tracy.extensions] + TracyProfilerExt = "TracyProfiler_jll" + + [deps.Tracy.weakdeps] + TracyProfiler_jll = "0c351ed6-8a68-550e-8b79-de6f926da83c" + +[[deps.TranscodingStreams]] +git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.11.3" + +[[deps.URIs]] +git-tree-sha1 = "bef26fb046d031353ef97a82e3fdb6afe7f21b1a" +uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" +version = "1.6.1" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.UnPack]] +git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b" +uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" +version = "1.0.2" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.Unitful]] +deps = ["Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "cec2df8cf14e0844a8c4d770d12347fda5931d72" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.25.0" + + [deps.Unitful.extensions] + ConstructionBaseUnitfulExt = "ConstructionBase" + ForwardDiffExt = "ForwardDiff" + InverseFunctionsUnitfulExt = "InverseFunctions" + LatexifyExt = ["Latexify", "LaTeXStrings"] + PrintfExt = "Printf" + + [deps.Unitful.weakdeps] + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" + Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" + Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.UnsafeAtomics]] +git-tree-sha1 = "b13c4edda90890e5b04ba24e20a310fbe6f249ff" +uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f" +version = "0.3.0" +weakdeps = ["LLVM"] + + [deps.UnsafeAtomics.extensions] + UnsafeAtomicsLLVM = ["LLVM"] + +[[deps.VectorizationBase]] +deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"] +git-tree-sha1 = "d1d9a935a26c475ebffd54e9c7ad11627c43ea85" +uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" +version = "0.21.72" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] +git-tree-sha1 = "80d3930c6347cfce7ccf96bd3bafdf079d9c0390" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.13.9+0" + +[[deps.XZ_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "fee71455b0aaa3440dfdd54a9a36ccef829be7d4" +uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800" +version = "5.8.1+0" + +[[deps.Xorg_libpciaccess_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] +git-tree-sha1 = "4909eb8f1cbf6bd4b1c30dd18b2ead9019ef2fad" +uuid = "a65dc6b1-eb27-53a1-bb3e-dea574b5389e" +version = "0.18.1+0" + +[[deps.ZipFile]] +deps = ["Libdl", "Printf", "Zlib_jll"] +git-tree-sha1 = "f492b7fe1698e623024e873244f10d89c95c340a" +uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" +version = "0.10.1" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+1" + +[[deps.Zstd_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "446b23e73536f84e8037f5dce465e92275f6a308" +uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" +version = "1.5.7+1" + +[[deps.demumble_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6498e3581023f8e530f34760d18f75a69e3a4ea8" +uuid = "1e29f10c-031c-5a83-9565-69cddfc27673" +version = "1.3.0+0" + +[[deps.libaec_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1aa23f01927b2dac46db77a56b31088feee0a491" +uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0" +version = "1.1.4+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.11.0+0" + +[[deps.libzip_jll]] +deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "OpenSSL_jll", "XZ_jll", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "86addc139bca85fdf9e7741e10977c45785727b7" +uuid = "337d8026-41b4-5cde-a456-74a10e5b31d1" +version = "1.11.3+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.52.0+1" + +[[deps.oneTBB_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "d5a767a3bb77135a99e433afe0eb14cd7f6914c3" +uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" +version = "2022.0.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+2" diff --git a/Project.toml b/Project.toml index 8a8396e72..c3f2a2e1f 100644 --- a/Project.toml +++ b/Project.toml @@ -56,6 +56,7 @@ ImageMorphology = "0.4" JLD2 = "0.4, 0.5, 0.6" KernelAbstractions = "0.9" MPI = "0.20" +MPIPreferences = "0.1.11" MeshArrays = "0.3" NCDatasets = "0.12, 0.13, 0.14" Oceananigans = "0.100.3" @@ -78,4 +79,4 @@ MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Coverage", "Test", "MPIPreferences", "CUDA_Runtime_jll", "Reactant", "PythonCall", "CondaPkg"] +test = ["Coverage", "Test", "MPIPreferences", "CUDA_Runtime_jll"] diff --git a/docs/make.jl b/docs/make.jl index 99afa9e46..1c48474ad 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,6 +5,22 @@ using DocumenterCitations, Literate +# temporary to enforce Oceananigans to use a specific branch +using Pkg +Pkg.add(PackageSpec( + name = "Oceananigans", + url = "https://github.com/CliMA/Oceananigans.jl", + rev = "ss/omip-branch-2" +)) + +Pkg.add(PackageSpec( + name = "ClimaSeaIce", + url = "https://github.com/CliMA/ClimaSeaIce.jl", + rev = "ss/omip-branch-2" +)) +Pkg.resolve() +@show Pkg.status() + ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" bib_filepath = joinpath(dirname(@__FILE__), "climaocean.bib") diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md index 047ccce2e..7d1a89e98 100644 --- a/docs/src/interface_fluxes.md +++ b/docs/src/interface_fluxes.md @@ -76,638 +76,3 @@ In both cases, computing turbulent fluxes requires: 5. Possibly, additional "bulk" properties of the surface media and radiation fluxes in order to compute an equilibrium "skin" surface temperature that differs from the bulk temperature below the surface. - -!!! note - In general, the surface specific humidity is typically related to the saturation specific humidity - at the the surface temperature ``T_s``, according to the Clausius-Claperyon relation. - For example, for ocean surfaces, the surface specific humidity is computed - according to via [Raoult's law](https://en.wikipedia.org/wiki/Raoult%27s_law) as - - ```math - q^\dagger(\rho, S, T) = x_{H_2O}(S) \frac{p_v^\dagger}{\rho R_v T} - ``` - - where ``x_{H_2O}(S)`` is the mole fraction of pure water in seawater with salinity ``S``, - and ``p_v^\dagger`` is the saturation vapor pressure, - - ```math - p_v^\dagger(T) = p_{tr} \left ( \frac{T}{T_{tr}} \right )^{\Delta c_p / R_v} \exp \left [ \frac{ℒ_{v0} - Δc_p T₀}{R_v} \left (\frac{1}{T_{tr}} - \frac{1}{T} \right ) \right ] - \quad \text{where} \quad - Δc_p = c_{p \ell} - c_{pv} \, . - ``` - - Many flux solvers (and the OMIP protocol) use a constant ``x_{H_2O} = 0.98``, which is equivalent to assuming - that the surface salinity is ``S \approx 35 \, \mathrm{g \, kg^{-1}}``, along with a reference seawater salinity composition. - Other surface specific humidity models may be used that take into account, for example, the microscopic structure - of snow, or the presence of a "dry skin" that buffers saturated soil from the atmosphere in a land model. - - Default values for the atmosphere thermodynamic parameters used to compute the saturation vapor pressure - and atmospheric equation of state are - - ```@example interface_fluxes - using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: AtmosphereThermodynamicsParameters - AtmosphereThermodynamicsParameters() - ``` - -### Coefficient-based fluxes - -Turbulent fluxes may be computed by prescribing "transfer coefficients" that relate differences -between the near-surface atmosphere and the ocean surface to fluxes, - -```math -\overline{\bm{u}' w'} ≈ C_D \, Δ \bm{u} \, U \\ -\overline{w' \theta'} ≈ C_\theta \, Δ \theta \, U \\ -\overline{w' q'} ≈ C_q \, Δ q \, U -``` - -The variable ``U`` is a characteristic velocity scale, which is most simply formulated as ``U = | Δ \bm{u}|``. -However, some parameterizations use formulations for ``U`` that -produce non-vanishing heat and moisture fluxes in zero-mean-wind conditions. -Usually these parameterizations are formulated as models for "gustiness" associated with atmospheric convection; -but more generally a common thread is that ``U`` may include contributions from unresolved turbulent motions -in addition to the relative mean velocity, ``Δ \bm{u}``. - -The variable ``C_D`` is often called the drag coefficient, while ``C_\theta`` and ``C_q`` are the heat transfer -coefficient and vapor flux coefficient. -The simplest method for computing fluxes is merely to prescribe ``C_D``, ``C_\theta``, and ``C_q`` -as constants -- typically with a magnitude around ``5 × 10^{-4}``--``2 × 10^{-3}``. -A comprehensive example is given below, but we note briefly here that -`ClimaOcean` supports the computation of turbulent fluxes with constant coefficients via - -```@example interface_fluxes -using ClimaOcean - -coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3, - heat_transfer_coefficient=2e-3, - vapor_flux_coefficient=1e-3) -``` - -### Similarity theory for neutral boundary layers - -The standard method for computing fluxes in realistic Earth system modeling contexts -uses a model for the structure of the near-surface atmosphere based on Monin--Obukhov similarity theory. -Similarity theory is essentially a dimensional argument and begins with the definition of "characteristic scales" -which are related to momentum, heat, and vapor fluxes through - -```math -u_\star^2 ≡ | \overline{\bm{u}' w'} | \\ -u_\star \theta_\star ≡ \overline{w' \theta'} \\ -u_\star q_\star ≡ \overline{w' q'} -``` - -where ``u_\star``, often called the "friction velocity", is the characteristic scale for velocity, -``\theta_\star`` is the characteristic scale for temperature, and ``q_\star`` is the characteristic scale -for water vapor. - -To introduce similarity theory, we first consider momentum fluxes in "neutral" conditions, -or with zero buoyancy flux. -We further simplify the situation by considering unidirectional flow with ``\bm{u} = u \, \bm{\hat x}``. -(To generalize our results to any flow direction, we simply rotate fluxes into the direction of the -relative velocity ``Δ \bm{u}``.) -The fundamental supposition of similarity theory is that the vertical shear depends only on -height above the boundary, such that by dimensional analysis, - -```math -\partial_z u \sim \frac{u_\star}{z} \, , -\qquad \text{and thus} \qquad -\partial_z u = \frac{u_\star}{\kappa z} \, , -``` - -where the second expression forms an equality by introducing the "Von Karman constant" ``\kappa``, -which is placed in the denominator by convention. -We can then integrate this expression from an inner scale ``z=\ell`` up to ``z=h`` to obtain - -```math -u_a(h) - u_a(\ell_u) = \frac{u_⋆}{\kappa} \log \left ( \frac{h}{\ell_u} \right ) -``` - -The inner length scale ``\ell_u``, which is called the "momentum roughness length", -can be interpreted as producing a characteristic upper value for the boundary layer shear, ``u_⋆ / \ell_u`` -in the region where similarity theory must be matched with the inner boundary layer (such as a viscous sublayer) -below. -Note that we take the inner velocity scale ``u_a(\ell_u)`` as being equal to the velocity of the surface, -so ``u_a(\ell_u) = u_s``. - -!!! note - We currently assume that the input to the surface flux computation is the - atmospheric velocity at ``z=h``. However, in coupled modeling context we are typically - instead given the atmosphere velocity _averaged_ over the height of the first layer, - or ``⟨u_a⟩_h = \frac{1}{h} \int_0^h \, u_a \, \mathrm{d} z``. - Formulating the flux computation in terms of ``⟨u_a⟩_h`` rather than ``u_a(h)`` - (e.g. [nishizawa2018surface](@citet)) - is a minor modification to the algorithm and an important avenue for future work. - -The roughness length in general depends on the physical nature of the surface. -For smooth, no-slip walls, experiments (cite) found agreement with a viscous sublayer model - -```math -\ell_ν = \mathbb{C}_\nu \frac{\nu}{u_\star} \, , -``` - -where ``\nu`` is the kinematic viscosity of the fluid (the air in our case) and ``\mathbb{C}_\nu`` is a free -parameter which was found to be around ``0.11``. -For air-water interfaces that develop a wind-forced spectrum of surface gravity waves, the alternative scaling - -```math -\ell_g = \mathbb{C}_g \frac{u_\star^2}{g} \, , -``` - -where ``g`` is gravitational acceleration, has been repeatedly (and perhaps shockingly due to its simplicity) confirmed by field campaigns. -The free parameter ``\mathbb{C}_g`` is often called - the "Charnock parameter" and takes typical values -between ``0`` and ``0.03`` [edson2013exchange](@citep). - -```@example -using ClimaOcean -using CairoMakie -set_theme!(Theme(fontsize=14, linewidth=4)) - -charnock_length = MomentumRoughnessLength(wave_formulation = 0.02, - smooth_wall_parameter = 0, - maximum_roughness_length = Inf) - -smooth_wall_length = MomentumRoughnessLength(wave_formulation = 0, - smooth_wall_parameter = 0.11) - -default_roughness_length = MomentumRoughnessLength() -modified_default_length = MomentumRoughnessLength(wave_formulation = 0.011) - -u★ = 1e-2:5e-3:3e-1 -ℓg = charnock_length.(u★) -ℓν = smooth_wall_length.(u★) -ℓd = default_roughness_length.(u★) -ℓ2 = modified_default_length.(u★) - -fig = Figure(size=(800, 400)) -ax1 = Axis(fig[1, 1], xlabel="Friction velocity, u★ (m s⁻¹)", ylabel="Momentum roughness length ℓᵤ (m)") -lines!(ax1, u★, ℓd, label="ClimaOcean default") -lines!(ax1, u★, ℓg, label="Charnock") -lines!(ax1, u★, ℓν, label="Smooth wall") -lines!(ax1, u★, ℓ2, color=:black, label="ClimaOcean default with ℂg = 0.011") - -ax2 = Axis(fig[1, 2], xlabel="Friction velocity, u★ (m s⁻¹)", ylabel="Momentum roughness length, ℓᵤ (m)") -u★ = 0.1:0.1:10 -ℓd = default_roughness_length.(u★) -ℓ2 = modified_default_length.(u★) -lines!(ax2, u★, ℓd) -lines!(ax2, u★, ℓ2, color=:black) - -Legend(fig[0, 1:2], ax1, orientation=:horizontal) - -fig -``` - -!!! note - The roughness length ``\ell`` should not be interpreted as a physical length scale, - a fact made clear by its submillimeter magnitude under (albeit calm) air-sea flux conditions. - -## Computing fluxes and the effective similarity drag coefficient - -ClimaOcean's default roughness length for air-sea fluxes is a function of the -friction velocity ``u_\star``. -This formulation produces a nonlinear equation for ``u_\star``, in terms of ``Δ u = u_a(h) - u_o``, -which we emphasize by rearranging the similarity profile - -```math -u_\star = \frac{\kappa \, Δ u}{\log \left [ h / \ell_u(u_\star) \right ]} \, . -``` - -The above equation is solved for ``u_\star`` using fixed-point iteration initialized with a reasonable -guess for ``u_\star``. -Once ``u_\star`` is obtained, the similarity drag coefficient may then be computed via - -```math -C_D(h) ≡ \frac{u_\star^2}{|Δ u(h)|^2} = \frac{\kappa^2}{\left ( \log \left [ h / \ell_u \right ] \right )^2} \, -``` - -where we have used the simple bulk velocity scale ``U = Δ u``. -We have also indicated that, the effective similarity drag "coefficient" depends on the height ``z=h`` -at which the atmospheric velocity is computed to form the relative velocity ``Δ u = u_a(h) - u_o``. -Most observational campaigns use ``h = 10 \, \mathrm{m}`` and most drag coefficients reported in the -literature pertain to ``h=10 \, \mathrm{m}``. - -To compute fluxes with ClimaOcean, we build an `OceanSeaIceModel` with an atmosphere and ocean state -concocted such that we can evaluate fluxes over a range of relative atmosphere and oceanic conditions. -For this we use a ``200 × 200`` horizontal grid and start with atmospheric winds that vary from -the relatively calm ``u_a(10 \, \mathrm{m}) = 0.5 \, \mathrm{m \, s^{-1}}`` to a -blustery ``u_a(10 \, \mathrm{m}) = 40 \, \mathrm{m \, s^{-1}}``. -We also initialize the ocean at rest with surface temperature ``T_o = 20 \, \mathrm{{}^∘ C}`` and -surface salinity ``S_o = 35 \, \mathrm{g \, kg^{-1}}`` -- but the surface temperature and salinity won't matter until later. - -```@example interface_fluxes -using Oceananigans -using ClimaOcean - -# Atmosphere velocities -Nx = Ny = 200 -uₐ = range(0.5, stop=40, length=Nx) # winds at 10 m, m/s - -# Ocean state parameters -T₀ = 20 # Surface temperature, ᵒC -S₀ = 35 # Surface salinity - -x = y = (0, 1) -z = (-1, 0) -atmos_grid = RectilinearGrid(size=(Nx, Ny); x, y, topology=(Periodic, Periodic, Flat)) -ocean_grid = RectilinearGrid(size=(Nx, Ny, 1); x, y, z, topology=(Periodic, Periodic, Bounded)) - -# Build the atmosphere -atmosphere = PrescribedAtmosphere(atmos_grid, surface_layer_height=10) -interior(atmosphere.tracers.T) .= 273.15 + T₀ # K -interior(atmosphere.velocities.u, :, :, 1, 1) .= uₐ # m/s - -kw = (momentum_advection=nothing, tracer_advection=nothing, closure=nothing) -ocean = ocean_simulation(ocean_grid; kw...) -set!(ocean.model, T=T₀, S=S₀) -``` - -Next we build two models with different flux formulations -- the default "similarity model" -that uses similarity theory with "Charnock" gravity wave parameter ``\mathbb{C}_g = 0.02``, -and a "coefficient model" with a constant drag coefficient ``C_D = 2 × 10^{-3}``: - -```@example interface_fluxes -neutral_similarity_fluxes = SimilarityTheoryFluxes(stability_functions=nothing) -interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=neutral_similarity_fluxes) -default_model = OceanSeaIceModel(ocean; atmosphere, interfaces) - -momentum_roughness_length = MomentumRoughnessLength(wave_formulation=0.04) -neutral_similarity_fluxes = SimilarityTheoryFluxes(stability_functions=nothing; momentum_roughness_length) -interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=neutral_similarity_fluxes) -increased_roughness_model = OceanSeaIceModel(ocean; atmosphere, interfaces) - -coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3) -interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=coefficient_fluxes) -coefficient_model = OceanSeaIceModel(ocean; atmosphere, interfaces) -``` - -Note that `OceanSeaIceModel` computes fluxes upon instantiation, so after constructing -the two models we are ready to analyze the results. -We first verify that the similarity model friction velocity has been computed successfully, - -```@example interface_fluxes -u★ = default_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity -u★ = interior(u★, :, 1, 1) -extrema(u★) -``` - -and it seems that we've obtained a range of friction velocities, which is expected -given that our atmospheric winds varied from ``0.5`` to ``40 \, \mathrm{m \, s^{-1}}``. -Computing the drag coefficient for the similarity model is as easy as - -```@example interface_fluxes -Cᴰ_default = @. (u★ / uₐ)^2 -extrema(Cᴰ_default) -``` - -We'll also re-compute the drag coefficient for the coefficient model -(which we specified as constant), which verifies that the coefficient was correctly -specified: - -```@example interface_fluxes -u★_coeff = coefficient_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity -u★_coeff = interior(u★_coeff, :, 1, 1) -Cᴰ_coeff = @. (u★_coeff / uₐ)^2 -extrema(Cᴰ_coeff) -``` - -We'll compare the computed fluxes and drag coefficients from our two models with -a polynomial expression due to [large2009global](@citet), and -an expression reported by [edson2013exchange](@citet) that was developed at ECMWF, - -```@example interface_fluxes -# From Large and Yeager (2009), equation 10 -c₁ = 0.0027 -c₂ = 0.000142 -c₃ = 0.0000764 -u★_LY = @. sqrt(c₁ * uₐ + c₂ * uₐ^2 + c₃ * uₐ^3) -Cᴰ_LY = @. (u★_LY / uₐ)^2 - -# From Edson et al. (2013), equation 20 -c₁ = 1.03e-3 -c₂ = 4e-5 -p₁ = 1.48 -p₂ = 0.21 -Cᴰ_EC = @. (c₁ + c₂ * uₐ^p₁) / uₐ^p₂ -u★_EC = @. sqrt(Cᴰ_EC) * uₐ -extrema(u★_EC) -``` - -Finally, we plot the results to compare the estimated friction velocity and effective -drag coefficient from the polynomials expressions with the two `OceanSeaIceModel`s: - -```@example interface_fluxes -using CairoMakie -set_theme!(Theme(fontsize=14, linewidth=4)) - -# Extract u★ and compute Cᴰ for increased roughness model -u★_rough = increased_roughness_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity -u★_rough = interior(u★_rough, :, 1, 1) -Cᴰ_rough = @. (u★_rough / uₐ)^2 - -fig = Figure(size=(800, 400)) -axu = Axis(fig[1:2, 1], xlabel="uₐ (m s⁻¹) at 10 m", ylabel="u★ (m s⁻¹)") -lines!(axu, uₐ, u★, label="ClimaOcean default") -lines!(axu, uₐ, u★_rough, label="Increased roughness model") -lines!(axu, uₐ, u★_LY, label="Large and Yeager (2009) polynomial fit") -lines!(axu, uₐ, u★_EC, label="ECMWF polynomial fit (Edson et al. 2013)") - -axd = Axis(fig[1:2, 2], xlabel="uₐ (m s⁻¹) at 10 m", ylabel="1000 × Cᴰ") -lines!(axd, uₐ, 1000 .* Cᴰ_default, label="ClimaOcean default") -lines!(axd, uₐ, 1000 .* Cᴰ_rough, label="Increased roughness model") -lines!(axd, uₐ, 1000 .* Cᴰ_LY, label="Large and Yeager (2009) polynomial fit") -lines!(axd, uₐ, 1000 .* Cᴰ_EC, label="ECMWF polynomial fit (Edson et al. 2013)") - -Legend(fig[3, 1:2], axd, nbanks = 2) - -fig -``` - -## Non-neutral boundary layers and stability functions - -The relationship between the relative air-sea state and turbulent fluxes -is modified by the presence of buoyancy fluxes -- "destabilizing" fluxes, which stimulate convection, -tend to increase turbulent exchange, while stabilizing fluxes suppress turbulence and turbulent exchange. -Monin--Obhukhov stability theory provides a scaling-argument-based framework -for modeling the effect of buoyancy fluxes on turbulent exchange. - -### Buoyancy flux and stability of the near-surface atmosphere - -Our next objective is to characterize the atmospheric statbility in terms of the buoyancy flux, ``\overline{w' b'}``, -which requires a bit of thermodynamics background to define the buoyancy perturbation, ``b'``. - -#### Buoyancy for a non-condensing mixture of dry air and water vapor - -The atmosphere is generally a mix of dry air, water vapor, non-vapor forms of water such as liquid droplets, -ice particles, rain, snow, hail, sleet, graupel, and more, and trace gases. -In the definition of buoyancy that follows, we neglect both the mass and volume of non-vapor water, -so that the specific humidity may be written - -```math -q \approx \frac{\rho_v}{\rho_v + \rho_d} \approx \frac{\rho_v}{\rho} \, , -``` - -where ``\rho_v`` is the density of water vapor, ``\rho_d`` is the density of dry air, and ``\rho \approx \rho_v + \rho_d`` -is the total density neglecting the mass of hypothetical condensed water species. - -!!! note - We endeavor to provide more information about the impact of this approximation. - Also, note that atmospheric data products like JRA55 do not explicitly provide - the mass ratio of condensed water, so the approximation is required in at least - some situations (such as simulations following the protocol of the Ocean Model - Intercomparison Project, OMIP). - On the other hand, generalizing the buoyancy formula that follow below to account - for the mass of condensed water is straightforward. - -The ideal gas law for a mixture of dry air and water vapor is then - -```math -p = \rho R_m(q) T \, -\qquad \text{where} \qquad -R_m(q) ≈ R_d \left (1 - q \right ) + R_v q = R_d \left ( 1 - \mathscr{M} q \right ) \, , -``` - -where ``\mathscr{M} = R_v/R_d - 1`` and ``R_m(q)`` is the effective mixture gas "constant" which varies with specific humidity ``q``, -and the ``\approx`` indicates that its form neglects the mass of condensed species. - -The buoyant perturbation experienced by air parcels advected by subgrid turbulent motions is then - -```math -b' ≡ - g \frac{\rho - \bar{\rho}}{\rho} = g \frac{\alpha - \bar{\alpha}}{\bar{\alpha}} -\qquad \text{where} \qquad -α = \frac{1}{\rho} = \frac{R_m T}{p} \, . -``` - -We neglect the effect of pressure perturbations to compute the buoyancy flux, so that ``p = \bar{p}`` and - -```math -\alpha - \bar{\alpha} = \frac{R_d}{p} \left [ T' - \mathscr{M} \left ( q' \bar{T} + \bar{q} T' + q' T' - \overline{q' T'} \right ) \right ] \, . -``` - -#### Buoyancy flux and the characteristic buoyancy scale - -In a computation whose details are reserved for an appendix, and which neglects ``\overline{q'T'}`` and the triple correlation ``\overline{w' q' T'}``, -we find that the buoyancy flux is approximately - -```math -\overline{w' b'} \approx g \frac{\overline{w'T'} - \mathscr{M} \left ( \overline{w' q'} \bar{T} + \bar{q} \overline{w' T'} \right )}{\bar{T} \left ( 1 - \mathscr{M} \bar q \right )} \, . -``` - -The characteristic buoyancy scale ``b_\star``, defined via ``u_\star b_\star \equiv \overline{w'b'}|_0``, is defined analogously to the temperature and vapor scales ``u_\star \theta_\star \equiv \overline{w' T'}`` and ``u_\star q_\star \equiv \overline{w' q'}``. -We therefore find - -```math -b_⋆ ≡ g \frac{\theta_\star - \mathscr{M} \left ( q_\star T_s + q_s \theta_\star \right ) }{ T_s \left (1 + \mathscr{M} q_s \right )} \, . -``` - -##### Stability of the near-surface atmosphere - -We use the ratio between the buoyancy flux and shear production at ``z=h`` -- which oceanographers often call -the "flux Richardson number", ``Ri_f`` -- to diagnose the stability of the atmosphere, - -```math -Ri_f(z) ≡ ζ(z) \equiv - \frac{\overline{w' b'}}{\partial_z \bar{\bm{u}} \, ⋅ \, \overline{\bm{u}' w'}} = - \frac{\kappa \, z}{u_\star^2} b_⋆ = \frac{z}{L_\star} -\qquad \text{where} \qquad -L_\star ≡ - \frac{u_\star^2}{\kappa b_\star} \, , -``` - -``\zeta`` is called the "stability parameter" and ``L_\star`` is called the "Monin--Obhukhov length scale". - -### The Monin--Obhukhov "stability functions" - -The fundamental premise of Monin--Obhkhov stability theory is that shear within a similarity layer affected by buoyancy fluxes may written - -```math -\frac{\kappa \, z}{u_\star} \partial_z \bar{u} = \tilde{\psi}_u(\zeta) \, , -``` - -where ``\tilde{\psi}_u(\zeta)`` is called the "stability function" (aka "dimensionless shear", and often denoted ``\phi``). -Comparing the Monin--Obukhov scaling to the neutral case expounded above, we find that ``\tilde{\psi}(0) = 1`` in neutral conditions with ``\zeta=0``. -In the presence of destabilizing fluxes, when ``ζ < 0``, observations show that ``\tilde{\psi}_u(\zeta) < 1`` (e.g. Businger 1971) -- representing an enhancement of turbulent exchange between the surface and atmosphere. -Conversely, ``\tilde{\psi}_u > 1`` when ``ζ > 0``, representing a suppression of turbulent fluxes by stable density stratification, or alternatively, an increase in the shear required to sustain a given friction velocity ``u_\star``. - -Monin and Obhukov's dimensional argument is also extended to potential temperature, so that for example - -```math -\frac{κ \, z}{\theta_\star} \partial_z \bar{\theta} = \tilde{\psi}_\theta (\zeta) \, . -``` - -Within the context of Monin--Obukhov stabilty theory, it can be shown that the neutral value ``\tilde{\psi}_\theta(0)`` is equal to the neutral turbulent Prandtl number, - -```math -Pr(\zeta=0) \equiv \frac{\tilde{\psi}_\theta(0)}{\tilde{\psi}_u(0)} = \tilde{\psi}_\theta(0) \, , -``` - -and observations suggest that ``\tilde{\psi}_θ(0) ≈ 0.7``. -Otherwise, the interpretation of variations in ``\tilde{\psi}_\theta`` (increased by stability, decreased by instability)is similar as for momentum. -We typically use the same "scalar" stability function to scale the vertical profiles of both temperature and water vapor, but neverthless ClimaOcean retains the possibility of an independent ``\tilde{\psi}_q``. - -### The Monin--Obhukhov self-similar vertical profiles - -To determine the implications of Monin--Obukhov similarity theory on the vertical profiles -of ``u``, ``\theta``, and ``q``, and therefore the implications for computing fluxes based on -the given differences ``Δ\bm{u}``, ``Δ \theta``, and ``Δ q``, we introduce "auxiliary stability functions" ``\psi_u(\zeta)``, which have derivatives ``\psi_u'(\zeta)`` and are related to ``\tilde{\psi}_u`` via - -```math -\tilde{ψ}_u(ζ) \equiv 1 - ζ ψ_u'(ζ) \, . -``` - -Inserting this transformation into the Monin--Obukhov scaling argument and rearranging terms yields - -```math -\partial_z u = \frac{u_\star}{\kappa \, z} + \frac{b_\star}{u_⋆} ψ' \left ( \frac{z}{L_⋆} \right ) \, , -``` - -which when integrated from ``z=\ell_u`` to ``z=h``, as for the neutral case, then produces - -```math -u_a(h) - u_a(\ell_u) = Δ u = \frac{u_\star}{\kappa} - \left [ \log \left (\frac{h}{\ell_u} \right ) - ψ_u \left ( \frac{h}{L_\star} \right ) + ψ_u \left (\frac{\ell_u}{L_\star} \right ) \right ] \, . -``` - -The term ``\psi_u(\ell_u / L_\star)`` is often neglected because ``\ell_u / L_\star`` is miniscule and because by definition, ``\psi_u(0) = 0``. -Similar formulas hold for temperature and water vapor, - -```math -Δ \theta = \frac{\theta_\star}{\kappa} \left [ \log \left (\frac{h}{\ell_\theta} \right ) - ψ_\theta \left ( \frac{h}{L_\star} \right ) + ψ_\theta \left (\frac{\ell_\theta}{L_\star} \right ) \right ] \, , \\[2ex] -Δ q = \frac{q_\star}{\kappa} \left [ \log \left (\frac{h}{\ell_q} \right ) - ψ_q \left ( \frac{h}{L_\star} \right ) + ψ_q \left (\frac{\ell_q}{L_\star} \right ) \right ] \, . -``` - -Let's plot some stability functions: - -```@example interface_fluxes -using ClimaOcean.OceanSeaIceModels.InterfaceComputations: - EdsonMomentumStabilityFunction, # Edson et al. 2013 - EdsonScalarStabilityFunction, # Edson et al. 2013 - ShebaMomentumStabilityFunction, # Grachev et al. 2007 - ShebaScalarStabilityFunction, # Grachev et al. 2007 - PaulsonMomentumStabilityFunction, # Paulson 1970 - PaulsonScalarStabilityFunction # Paulson 1970 - -edson_momentum = EdsonMomentumStabilityFunction() -edson_scalar = EdsonScalarStabilityFunction() -sheba_momentum = ShebaMomentumStabilityFunction() -sheba_scalar = ShebaScalarStabilityFunction() -paulson_momentum = PaulsonMomentumStabilityFunction() -paulson_scalar = PaulsonScalarStabilityFunction() - -ζstep = 0.01 -ζ = -4:ζstep:4 -ζ⁺ = first(ζ[ζ .≥ 0]):ζstep:last(ζ) -ζ⁻ = first(ζ):ζstep:last(ζ[ζ .≤ 0]) - -fig = Figure(size=(800, 400)) - -axm = Axis(fig[1, 1], xlabel="Stability parameter ζ", ylabel="Momentum auxiliary stability function ψₘ") -axs = Axis(fig[1, 2], xlabel="Stability parameter ζ", ylabel="Scalar auxiliary stability function ψₛ") - -lines!(axm, ζ, edson_momentum.(ζ), label="Edson et al. (2013)", alpha=0.7) -lines!(axm, ζ⁺, sheba_momentum.(ζ⁺), label="Grachev et al. (2007)", alpha=0.7) -lines!(axm, ζ⁻, paulson_momentum.(ζ⁻), label="Paulson (1970)", alpha=0.7) -axislegend(axm, position=:lb) - -lines!(axs, ζ, edson_scalar.(ζ), label="Edson et al. (2013)", alpha=0.7) -lines!(axs, ζ⁺, sheba_scalar.(ζ⁺), label="Grachev et al. (2007)", alpha=0.7) -lines!(axs, ζ⁻, paulson_scalar.(ζ⁻), label="Paulson (1970)", alpha=0.7) - -for ax in (axm, axs) - ylims!(ax, -14, 4) -end - -fig -``` - -#### Computing fluxes given atmopshere, surface, and bulk interior states - -We compute surface fluxes by solving the nonlinear set of equations for ``u_\star``, ``\theta_\star``. -We use fixed point iteration of the following three-variable system, - -```math -u_⋆^{n+1} = \, Δ u \, \, Ξ_u \left (h, \ell_u^n, L_⋆^n \right ) \\[2ex] -θ_⋆^{n+1} = \, Δ θ \, \, Ξ_θ \left (h, \ell_θ^n, L_⋆^n \right ) \\[2ex] -q_⋆^{n+1} = \, Δ q \, \, Ξ_q \left (h, \ell_q^n, L_⋆^n \right ) -``` - -where, for example, - -```math -\Xi_u \left ( h, \ell_u, L_⋆ \right ) ≡ \frac{κ}{\log \left ( \frac{h}{\ell_u} \right ) - \psi_u \left ( \frac{h}{L_\star} \right ) + \psi_u \left ( \frac{\ell_u}{L_\star} \right )} \, , -``` - -The above equations indicate how ``\ell_u``, ``\ell_\theta``, ``\ell_q``, and ``L_⋆ = - u_\star^2 / κ b_\star`` are all functions of ``u_\star, \theta_\star, q_\star``;\ -estimating the right-hand side requires using values at the previous iterate ``n``. -Note that if a skin temperature model is used, then we obtain a four-variable system, - -```math -u_⋆^{n+1} = \, Δ u \, \, Ξ_u \left (h, \ell_u^n, L_⋆^n \right ) \\[2ex] -θ_⋆^{n+1} = \, Δ θ^n \, \, Ξ_θ \left (h, \ell_θ^n, L_⋆^n \right ) \\[2ex] -q_⋆^{n+1} = \, Δ q^n \, \, Ξ_q \left (h, \ell_q^n, L_⋆^n \right ) \\[2ex] -T_s^{n+1} = F_T \left (θ_⋆, q_⋆, I_{sw}, I_{lw}, \cdots \right ) -``` - -where ``F_T`` denotes an estimate of the surface temperature that, in general, requires all incoming heat fluxes -including shortwave and longwave radiation ``I_{sw}`` and ``I_{lw}``. -In the skin temperature case, the air-surface temperature difference ``Δ \theta`` and the saturation specific humidity -that enters into the air-surface specific humidity difference ``Δ q`` also change each iterate. - -```@example interface_fluxes -using ClimaOcean.OceanSeaIceModels.InterfaceComputations: surface_specific_humidity - -ρₐ = 1.2 # guess -Tₒ = 273.15 + 20 # in Kelvin -Sₒ = 35 -interfaces = default_model.interfaces -ℂₐ = interfaces.atmosphere_properties -q_formulation = interfaces.atmosphere_ocean_interface.properties.specific_humidity_formulation -qₛ = surface_specific_humidity(q_formulation, ℂₐ, ρₐ, Tₒ, Sₒ) -@show qₛ -``` - -We then set the atmospheric state: - -```@example interface_fluxes -interior(atmosphere.pressure) .= 101352 -interior(atmosphere.tracers.q) .= qₛ - -Tₐ = 273.15 .+ range(-40, stop=40, length=Ny) -Tₐ = reshape(Tₐ, 1, Ny) -interior(atmosphere.tracers.T) .= Tₐ - -Oceananigans.TimeSteppers.update_state!(default_model) -u★ = default_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity -θ★ = default_model.interfaces.atmosphere_ocean_interface.fluxes.temperature_scale - -fig = Figure(size=(800, 600)) -axu = Axis(fig[2, 1], xlabel="Wind speed uₐ (m s⁻¹)", ylabel="Air-sea temperature difference (K)") -axθ = Axis(fig[2, 2], xlabel="Wind speed uₐ (m s⁻¹)", ylabel="Air-sea temperature difference (K)") -axC = Axis(fig[3, 1:2], xlabel="Wind speed uₐ (m s⁻¹)", ylabel="Cᴰ / neutral Cᴰ") - -ΔT = Tₐ .- Tₒ -ΔT = dropdims(ΔT, dims=1) - -hmu = heatmap!(axu, uₐ, ΔT, u★, colormap=:speed) -hmθ = heatmap!(axθ, uₐ, ΔT, θ★, colormap=:balance) - -Colorbar(fig[1, 1], hmu, label="u★ (m s⁻¹)", vertical=false) -Colorbar(fig[1, 2], hmθ, label="θ★ (K)", vertical=false) - -Cᴰ = [(u★[i, j] / uₐ[i])^2 for i in 1:Nx, j in 1:Ny] - -for j in (1, 20, 50, 100, 150, 200) - lines!(axC, uₐ, Cᴰ[:, j] ./ Cᴰ_default, label="ΔT = $(round(ΔT[j], digits=1)) K", alpha=0.8) -end - -axislegend(axC, orientation=:horizontal, nbanks=2, label="navid") - -xlims!(axC, 0, 10) -ylims!(axC, 0, 4) - -fig -``` - -The coefficient-based formula then takes the form - -```math -u_\star = \sqrt{C_D | Δ \bm{u} | \, U} \\ -\theta_\star = \frac{C_θ}{\sqrt{C_D}} \, Δ θ \, \sqrt{\frac{U}{|Δ \bm{u} |}} \\ -q_\star = \frac{C_q}{\sqrt{C_D}} \, Δ q \, \sqrt{\frac{U}{| Δ \bm{u} |}} \\ -``` diff --git a/experiments/omip_prototype/download_data.jl b/experiments/omip_prototype/download_data.jl new file mode 100644 index 000000000..2881bcd71 --- /dev/null +++ b/experiments/omip_prototype/download_data.jl @@ -0,0 +1,8 @@ +using ClimaOcean +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling: download_dataset + +atmosphere = JRA55PrescribedAtmosphere(; dir="forcing_data/", + dataset=JRA55MultipleYears(), + backend=JRA55NetCDFBackend(10), + include_rivers_and_icebergs=true) \ No newline at end of file diff --git a/experiments/omip_prototype/half_degree_simulation.jl b/experiments/omip_prototype/half_degree_simulation.jl new file mode 100644 index 000000000..45ac70037 --- /dev/null +++ b/experiments/omip_prototype/half_degree_simulation.jl @@ -0,0 +1,450 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency +using ArgParse + +import Oceananigans.OutputWriters: checkpointer_address + +function parse_commandline() + s = ArgParseSettings() + + @add_arg_table! s begin + "--kappa_skew" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--kappa_symmetric" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + end + return parse_args(s) +end + +args = parse_commandline() +κ_skew = args["kappa_skew"] +κ_symmetric = args["kappa_symmetric"] + +arch = GPU() + +Nx = 720 # longitudinal direction +Ny = 360 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0) + +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=40) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +@info "Built grid $(grid)" + +##### +##### A Prognostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity + +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8, fixed_Δt=65minutes) + +gm_closure = Oceananigans.TurbulenceClosures.EddyAdvectiveClosure(; κ_skew) +redi_closure = Oceananigans.TurbulenceClosures.IsopycnalDiffusivity(; κ_symmetric) + +# obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +obl_closure = RiBasedVerticalDiffusivity() +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-4), gm_closure, redi_closure) + +prefix = "halfdegree" +if obl_closure isa RiBasedVerticalDiffusivity + prefix *= "_RiBased" +else + prefix *= "_CATKE" +end + +prefix *= "_$(κ_skew)_$(κ_symmetric)" +prefix *= "_newgm_multiyearjra55" + +dir = joinpath(homedir(), "forcing_data_half_degree") +mkpath(dir) + +start_date = DateTime(1958, 1, 1) +end_date = start_date + Year(40) +simulation_period = Dates.value(Second(end_date - start_date)) +yearly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(1):end_date))))) +decadal_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(10):end_date))))) + +@info "Settting up salinity restoring..." +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset=EN4Monthly(), dir, start_date, end_date) +FS = DatasetRestoring(Smetadata, grid; rate = 1/18days, mask, time_indices_in_memory = 3) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +@info "Built ocean model $(ocean)" + +set!(ocean.model, T=Metadatum(:temperature; dataset=EN4Monthly(), date=start_date, dir), + S=Metadatum(:salinity; dataset=EN4Monthly(), date=start_date, dir)) +@info "Initialized T and S" + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) +@info "Built sea ice model $(sea_ice)" + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir)) + +@info "Initialized sea ice fields" + +##### +##### A Prescribed Atmosphere model +##### +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(6) + +@info "Setting up presctibed atmosphere $(dataset)" +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +@info "Built atmosphere model $(atmosphere)" + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built coupled model $(omip)" + +omip = Simulation(omip, Δt=20minutes, stop_time=60days) +@info "Built simulation $(omip)" + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1825days), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:time_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(1825days, window=1825days), + filename = "$(FILE_DIR)/ocean_complete_fields_10year_average", + overwrite_existing = true) + +sea_ice.output_writers[:time_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(1825days, window=1825days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_10year_average", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(100)) + +run!(omip) + +omip.Δt = 60minutes +omip.stop_time = simulation_period + +run!(omip) + +#%% +@info "Plotting results..." +using CairoMakie + +uo = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "u"; backend = OnDisk()) +vo = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "v"; backend = OnDisk()) +To = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "T"; backend = OnDisk()) + +# and sea ice fields with "i": +ui = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "u"; backend = OnDisk()) +vi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "v"; backend = OnDisk()) +hi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "h"; backend = OnDisk()) +ℵi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "ℵ"; backend = OnDisk()) +Ti = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "T"; backend = OnDisk()) + +times = uo.times +Nt = length(times) +n = Observable(Nt) + +# We create a land mask and use it to fill land points with `NaN`s. +land = interior(To.grid.immersed_boundary.bottom_height) .≥ 0 + +Toₙ = @lift begin + Tₙ = interior(To[$n]) + Tₙ[land] .= NaN + view(Tₙ, :, :, 1) +end + +heₙ = @lift begin + hₙ = interior(hi[$n]) + ℵₙ = interior(ℵi[$n]) + hₙ[land] .= NaN + view(hₙ, :, :, 1) .* view(ℵₙ, :, :, 1) +end + +# We compute the surface speeds for the ocean and the sea ice. +uoₙ = Field{Face, Center, Nothing}(uo.grid) +voₙ = Field{Center, Face, Nothing}(vo.grid) + +uiₙ = Field{Face, Center, Nothing}(ui.grid) +viₙ = Field{Center, Face, Nothing}(vi.grid) + +so = Field(sqrt(uoₙ^2 + voₙ^2)) +si = Field(sqrt(uiₙ^2 + viₙ^2)) + +soₙ = @lift begin + parent(uoₙ) .= parent(uo[$n]) + parent(voₙ) .= parent(vo[$n]) + compute!(so) + soₙ = interior(so) + soₙ[land] .= NaN + view(soₙ, :, :, 1) +end + +siₙ = @lift begin + parent(uiₙ) .= parent(ui[$n]) + parent(viₙ) .= parent(vi[$n]) + compute!(si) + siₙ = interior(si) + hₙ = interior(hi[$n]) + ℵₙ = interior(ℵi[$n]) + he = hₙ .* ℵₙ + siₙ[he .< 1e-7] .= 0 + siₙ[land] .= NaN + view(siₙ, :, :, 1) +end + +# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent +# eddy kinetic energy from the CATKE vertical mixing parameterization as well as the +# sea ice speed and the effective sea ice thickness. +fig = Figure(size=(1500, 1000)) + +title = @lift string("Global 1ᵒ ocean simulation after ", prettytime(times[$n] - times[1])) + +axso = Axis(fig[1, 1]) +axsi = Axis(fig[1, 3]) +axTo = Axis(fig[2, 1]) +axhi = Axis(fig[2, 3]) + +hmo = heatmap!(axso, soₙ, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) +hmi = heatmap!(axsi, siₙ, colorrange = (0, 0.5), colormap = :greys, nan_color=:lightgray) +Colorbar(fig[1, 2], hmo, label = "Ocean Surface speed (m s⁻¹)") +Colorbar(fig[1, 4], hmi, label = "Sea ice speed (m s⁻¹)") + +hmo = heatmap!(axTo, Toₙ, colorrange = (-1, 32), colormap = :magma, nan_color=:lightgray) +hmi = heatmap!(axhi, heₙ, colorrange = (0, 4), colormap = :blues, nan_color=:lightgray) +Colorbar(fig[2, 2], hmo, label = "Surface Temperature (ᵒC)") +Colorbar(fig[2, 4], hmi, label = "Effective ice thickness (m)") + +for ax in (axso, axsi, axTo, axhi) + hidedecorations!(ax) +end + +Label(fig[0, :], title) + +save("$(FILE_DIR)/global_snapshot.png", fig) +nothing #hide + +# ![](global_snapshot.png) + +# And now a movie: + +CairoMakie.record(fig, "$(FILE_DIR)/global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn + n[] = nn +end +nothing + +#%% +uoac = FieldTimeSeries("$(FILE_DIR)/ocean_complete_fields_10year_average.jld2", "u"; backend = OnDisk()) +voac = FieldTimeSeries("$(FILE_DIR)/ocean_complete_fields_10year_average.jld2", "v"; backend = OnDisk()) +woac = FieldTimeSeries("$(FILE_DIR)/ocean_complete_fields_10year_average.jld2", "w"; backend = OnDisk()) +Toac = FieldTimeSeries("$(FILE_DIR)/ocean_complete_fields_10year_average.jld2", "T"; backend = OnDisk()) +Soac = FieldTimeSeries("$(FILE_DIR)/ocean_complete_fields_10year_average.jld2", "S"; backend = OnDisk()) + +times = uoac.times +Nt = length(times) +n = Observable(Nt) + +# We create a land mask and use it to fill land points with `NaN`s. +land = interior(Toac.grid.immersed_boundary.bottom_height) .≥ 0 + +Toacₙ = @lift begin + Tₙ = interior(Toac[$n]) + Tₙ[land[:, :, 1], :] .= NaN + view(Tₙ, :, :, :) +end + +Soacₙ = @lift begin + Sₙ = interior(Soac[$n]) + Sₙ[land[:, :, 1], :] .= NaN + view(Sₙ, :, :, :) +end + +surface = Toac.grid.underlying_grid.Nz +middle = findlast(x -> x <= -500, Toac.grid.underlying_grid.z.cᵃᵃᶜ) +bottom = findlast(x -> x <= -2000, Toac.grid.underlying_grid.z.cᵃᵃᶜ) + +T_surfaceₙ = @lift $Toacₙ[:, :, surface] +T_middleₙ = @lift $Toacₙ[:, :, middle] +T_bottomₙ = @lift $Toacₙ[:, :, bottom] + +S_surfaceₙ = @lift $Soacₙ[:, :, surface] +S_middleₙ = @lift $Soacₙ[:, :, middle] +S_bottomₙ = @lift $Soacₙ[:, :, bottom] + +Tlim_surface = extrema(interior(Toac[Nt], :, :, surface)) +Tlim_middle = extrema(interior(Toac[Nt], :, :, middle)) +Tlim_bottom = extrema(interior(Toac[Nt], :, :, bottom)) + +Slim_surface = extrema(interior(Soac[Nt], :, :, surface)[interior(Soac[Nt], :, :, surface) .!= 0]) +Slim_middle = extrema(interior(Soac[Nt], :, :, middle)[interior(Soac[Nt], :, :, middle) .!= 0]) +Slim_bottom = extrema(interior(Soac[Nt], :, :, bottom)[interior(Soac[Nt], :, :, bottom) .!= 0]) + +# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent +# eddy kinetic energy from the CATKE vertical mixing parameterization as well as the +# sea ice speed and the effective sea ice thickness. +fig = Figure(size=(2400, 1000)) + +title = @lift string("10 year average starting from year ", round((times[$n] - times[1]) / 365 / 24 / 60^2)) + +axTs = Axis(fig[1, 1], title = "z = $(round(Toac.grid.underlying_grid.z.cᵃᵃᶜ[surface])) m") +axTm = Axis(fig[1, 3], title = "z = $(round(Toac.grid.underlying_grid.z.cᵃᵃᶜ[middle])) m") +axTb = Axis(fig[1, 5], title = "z = $(round(Toac.grid.underlying_grid.z.cᵃᵃᶜ[bottom])) m") + +axSs = Axis(fig[2, 1]) +axSm = Axis(fig[2, 3]) +axSb = Axis(fig[2, 5]) + +T_colorscheme = :turbo +S_colorscheme = :turbo + +hmTs = heatmap!(axTs, T_surfaceₙ, colormap = T_colorscheme, nan_color=:lightgray, colorrange = Tlim_surface) +hmTm = heatmap!(axTm, T_middleₙ, colormap = T_colorscheme, nan_color=:lightgray, colorrange = Tlim_middle) +hmTb = heatmap!(axTb, T_bottomₙ, colormap = T_colorscheme, nan_color=:lightgray, colorrange = Tlim_bottom) +Colorbar(fig[1, 2], hmTs) +Colorbar(fig[1, 4], hmTm) +Colorbar(fig[1, 6], hmTb, label = "Temperature (ᵒC)") + +hmSs = heatmap!(axSs, S_surfaceₙ, colormap = S_colorscheme, nan_color=:lightgray, colorrange = Slim_surface) +hmSm = heatmap!(axSm, S_middleₙ, colormap = S_colorscheme, nan_color=:lightgray, colorrange = Slim_middle) +hmSb = heatmap!(axSb, S_bottomₙ, colormap = S_colorscheme, nan_color=:lightgray, colorrange = Slim_bottom) +Colorbar(fig[2, 2], hmSs) +Colorbar(fig[2, 4], hmSm) +Colorbar(fig[2, 6], hmSb, label = "Salinity (psu)") + +for ax in (axTs, axTm, axTb, axSs, axSm, axSb) + hidedecorations!(ax) +end + +Label(fig[0, :], title) + +save("$(FILE_DIR)/10year_time_average_level_snapshot.png", fig) +nothing #hide + +# ![](global_snapshot.png) + +# And now a movie: + +CairoMakie.record(fig, "$(FILE_DIR)/10year_time_average_level.mp4", 1:Nt, framerate = 1) do nn + n[] = nn +end +nothing +#%% \ No newline at end of file diff --git a/experiments/omip_prototype/minimal_script_MWE.jl b/experiments/omip_prototype/minimal_script_MWE.jl new file mode 100644 index 000000000..c6c7f2bf9 --- /dev/null +++ b/experiments/omip_prototype/minimal_script_MWE.jl @@ -0,0 +1,172 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using Oceananigans.TurbulenceClosures: DiffusiveFormulation, AdvectiveFormulation +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency + +import Oceananigans.OutputWriters: checkpointer_address + +CUDA.versioninfo() + +ngpus = Int(length(CUDA.devices())) +@info "$ngpus GPUs" +# arch = GPU() +arch = Distributed(GPU(), partition=Partition(1, ngpus), synchronized_communication=true) +# arch = Distributed(GPU(), partition=Partition(1, ngpus)) +# arch = Distributed(CPU(), partition=Partition(1, 4), synchronized_communication=true) +@info "Architecture $(arch)" + +Nx = 2880 # longitudinal direction +Ny = 1440 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0) + +const z_surf = z_faces(Nz) +halo_size = 7 + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (halo_size, halo_size, halo_size)) + +@info "halo size: $(halo_size)" + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=10) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) +# grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=false) + +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; substeps = 200) + +using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity + +obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +# obl_closure = RiBasedVerticalDiffusivity() +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-4)) + +dir = joinpath(homedir(), "forcing_data_minimal") +mkpath(dir) + +dataset = EN4Monthly() +start_date = DateTime(1993, 1, 1) +end_date = start_date + Month(3) + +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset, dir) + +ocean = ocean_simulation(grid; Δt=5minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + closure) + +dataset = EN4Monthly() + +set!(ocean.model, T=Metadatum(:temperature; dataset, date=start_date, dir), + S=Metadatum(:salinity; dataset, date=start_date, dir)) + +@info ocean.model.clock + +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir, date=start_date), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir, date=start_date)) + +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(10) + +@info "dataset: $dataset" + +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) + +@info "atmosphere: $atmosphere" +radiation = Radiation() + +@info "Setting up Ocean-SeaIce-Atmosphere model..." +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +# omip = OceanSeaIceModel(ocean; atmosphere, radiation) +@info "sea ice model: $(omip.sea_ice)" +@info omip + +@info "Setting up OMIP simulation..." +omip = Simulation(omip, Δt=5minutes, stop_time=60days) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(interior(sea_ice.model.ice_thickness)) + ℵmax = maximum(interior(sea_ice.model.ice_concentration)) + Tmax = maximum(interior(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)) + Tmin = minimum(interior(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)) + # Tmax = maximum(ocean.model.tracers.T) + # Tmin = minimum(ocean.model.tracers.T) + Smax = maximum(interior(ocean.model.tracers.S)) + Smin = minimum(interior(ocean.model.tracers.S)) + umax = maximum(interior(ocean.model.velocities.u)) + vmax = maximum(interior(ocean.model.velocities.v)) + wmax = maximum(interior(ocean.model.velocities.w)) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + # msg0 = @sprintf("local rank: %d, ", arch.local_rank) + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("extrema(S): (%.2f, %.2f) psu, ", Smax, Smin) + msg6 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg7 = @sprintf("wall time: %s \n", prettytime(step_time)) + + # @info msg0 * msg1 + @info msg1 * msg2 * msg4 * msg5 * msg6 * msg7 + # @info msg0 * msg1 * msg2 * msg4 * msg5 * msg6 * msg7 + # @info msg1 * msg2 * msg4 * msg5 * msg6 * msg7 + # @info msg1 * msg4 * msg5 * msg6 * msg7 + + CUDA.memory_status() + wall_time[] = time_ns() + return nothing +end + +function gpus_gc(sim) + @info "Starting GPU garbage collection for iteration $(iteration(sim))..." + @info "before GC:" + CUDA.memory_status() + wall_time[] = time_ns() + + GC.gc() + CUDA.reclaim() + @info "after GC:" + CUDA.memory_status() + step_time = 1e-9 * (time_ns() - wall_time[]) + @info "GPU garbage collection took $(prettytime(step_time))" + wall_time[] = time_ns() + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(10)) +# add_callback!(omip, gpus_gc, IterationInterval(10)) + +@info "Starting OMIP simulation..." +run!(omip) \ No newline at end of file diff --git a/experiments/omip_prototype/one_degree_omip.jl b/experiments/omip_prototype/one_degree_omip.jl new file mode 100644 index 000000000..1744af824 --- /dev/null +++ b/experiments/omip_prototype/one_degree_omip.jl @@ -0,0 +1,155 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA + +import Oceananigans.OutputWriters: checkpointer_address + +arch = GPU() + +Nx = 360 # longitudinal direction +Ny = 180 # meridional direction +Nz = 60 + +z_faces = ExponentialCoordinate(Nz, -6000, 0) +# z_faces = MutableVerticalDiscretization(z_faces) + +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=75) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=5) + +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7, fixed_Δt=45minutes) + +eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3) +catke_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +closure = (catke_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-4), eddy_closure) + +dataset = EN4Monthly() +date = DateTime(1958, 1, 1) +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset) + +FS = DatasetRestoring(Smetadata, grid; rate = 1/18days, mask, time_indices_in_memory = 10) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +dataset = EN4Monthly() +date = DateTime(1958, 1, 1) + +set!(ocean.model, T=Metadatum(:temperature; dataset, date), + S=Metadatum(:salinity; dataset, date)) + +@info ocean.model.clock + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly()), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly())) + +##### +##### A Prescribed Atmosphere model +##### + +dir = "./forcing_data" +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +atmosphere = JRA55PrescribedAtmosphere(arch; dir, dataset, backend, include_rivers_and_icebergs=true) +radiation = Radiation() + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +omip = Simulation(omip, Δt=20minutes, stop_time=60days) + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = IterationInterval(1000), + prefix = "ocean_checkpoint_onedegree", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = IterationInterval(1000), + prefix = "sea_ice_checkpoint_onedegree", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(1)) + +run!(omip) + +omip.Δt = 30minutes +omip.stop_time = 58 * 365days + +run!(omip) diff --git a/experiments/omip_prototype/oneeighth_degree_simulation.jl b/experiments/omip_prototype/oneeighth_degree_simulation.jl new file mode 100644 index 000000000..8db222c29 --- /dev/null +++ b/experiments/omip_prototype/oneeighth_degree_simulation.jl @@ -0,0 +1,348 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using Oceananigans.Architectures: on_architecture +using ClimaOcean.OceanSimulations +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaOcean.DataWrangling: NearestNeighborInpainting +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using PythonCall +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency + +import Oceananigans.OutputWriters: checkpointer_address + +CUDA.versioninfo() + +ngpus = Int(length(CUDA.devices())) +@info "$ngpus GPUs" +# arch = GPU() +arch = Distributed(GPU(), partition=Partition(1, ngpus), synchronized_communication=true) +# arch = Distributed(GPU(), partition=Partition(1, ngpus)) +# arch = Distributed(CPU(), partition=Partition(1, 4), synchronized_communication=true) +@info "Architecture $(arch)" + +Nx = 2880 # longitudinal direction +Ny = 1440 # meridional direction +Nz = 100 + +Δt = 10minutes + +z_faces = ExponentialDiscretization(Nz, -6000, 0) + +const z_surf = z_faces(Nz) + +@info "Building grid..." +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +@info "Regridding bathymetry..." +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=10) +fitted_bottom = GridFittedBottom(bottom_height) + +@info "Building immersed boundary grid..." +grid = ImmersedBoundaryGrid(grid, fitted_bottom; active_cells_map=true) +# grid = ImmersedBoundaryGrid(grid, fitted_bottom) +@info grid +@info "Created ImmersedBoundaryGrid" + +##### +##### A Propgnostic Ocean model +##### +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; substeps = 200) +@info "Free surface", free_surface + +obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() # CATKE +# obl_closure = RiBasedVerticalDiffusivity() +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-4)) + +if obl_closure isa RiBasedVerticalDiffusivity + prefix = "RiBased" +else + prefix = "CATKE" +end + +prefix *= "oneeighth_degree_$(Δt)" + +dir = joinpath(homedir(), "oneeighth_degree_forcing_data") +mkpath(dir) + +@info "Building ocean component..." +ocean = ocean_simulation(grid; Δt=5minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + closure) + +start_date = DateTime(1993, 1, 1) +end_date = DateTime(2003, 4, 1) +simulation_period = Dates.value(Second(end_date - start_date)) +monthly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Month(1):end_date))))) + +dataset = EN4Monthly() + +@info "Setting initial conditions..." +set!(ocean.model, T=Metadatum(:temperature; dataset, date=start_date, dir), + S=Metadatum(:salinity; dataset, date=start_date, dir)) + +@info ocean.model.clock + +##### +##### A Prognostic Sea-ice model +##### + +@info "Building sea-ice component..." +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) + +@info "Setting sea-ice initial conditions..." +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir, date=start_date), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir, date=start_date)) + +##### +##### A Prescribed Atmosphere model +##### + +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(20) + +@info "Building atmospheric forcing..." +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +##### +##### An ocean-sea ice coupled model +##### + +@info "Building coupled ocean-sea ice model..." +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +omip = Simulation(omip, Δt=5minutes, stop_time=60days) + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +@info "Setting up output writers..." +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(monthly_times), + prefix = "$(FILE_DIR)/ocean_checkpoint_oneeighthdegree", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(monthly_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint_oneeighthdegree", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) + +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(30days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(30days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = SpecifiedTimes(monthly_times), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +sea_ice.output_writers[:full] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = SpecifiedTimes(monthly_times), + filename = "$(FILE_DIR)/sea_ice_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:monthly_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(30days, window=30days), + filename = "$(FILE_DIR)/ocean_complete_fields_monthly_average", + overwrite_existing = true) + +sea_ice.output_writers[:monthly_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(30days, window=30days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_monthly_average", + overwrite_existing = true) + +ocean.output_writers[:yearly_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(365days, window=365days), + filename = "$(FILE_DIR)/ocean_complete_fields_yearly_average", + overwrite_existing = true) + +sea_ice.output_writers[:yearly_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(365days, window=365days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_yearly_average", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("local rank: %d, ", arch.local_rank) + msg2 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg3 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 + + CUDA.memory_status() + + total_mem = Sys.total_memory() + free_mem = Sys.free_memory() + used_mem = total_mem - free_mem + @info msg1 * "Memory usage: $(round(100 * used_mem / total_mem, digits=1))%, $(round(used_mem / 1e9, digits=1)) / $(round(total_mem / 1e9, digits=1)) GB" + wall_time[] = time_ns() + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(400)) + +@info "Starting simulation..." +run!(omip) + +@info "Initialization complete. Running the rest..." + +omip.Δt = Δt +omip.stop_time = simulation_period + +run!(omip) + +# #%% +# @info "Plotting results..." +# using CairoMakie + +# uo = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "u"; backend = OnDisk()) +# vo = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "v"; backend = OnDisk()) +# To = FieldTimeSeries("$(FILE_DIR)/ocean_surface_fields.jld2", "T"; backend = OnDisk()) + +# # and sea ice fields with "i": +# ui = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "u"; backend = OnDisk()) +# vi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "v"; backend = OnDisk()) +# hi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "h"; backend = OnDisk()) +# ℵi = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "ℵ"; backend = OnDisk()) +# Ti = FieldTimeSeries("$(FILE_DIR)/sea_ice_surface_fields.jld2", "T"; backend = OnDisk()) + +# times = uo.times +# Nt = length(times) +# n = Observable(Nt) + +# # We create a land mask and use it to fill land points with `NaN`s. +# land = interior(To.grid.immersed_boundary.bottom_height) .≥ 0 + +# Toₙ = @lift begin +# Tₙ = interior(To[$n]) +# Tₙ[land] .= NaN +# view(Tₙ, :, :, 1) +# end + +# heₙ = @lift begin +# hₙ = interior(hi[$n]) +# ℵₙ = interior(ℵi[$n]) +# hₙ[land] .= NaN +# view(hₙ, :, :, 1) .* view(ℵₙ, :, :, 1) +# end + +# # We compute the surface speeds for the ocean and the sea ice. +# uoₙ = Field{Face, Center, Nothing}(uo.grid) +# voₙ = Field{Center, Face, Nothing}(vo.grid) + +# uiₙ = Field{Face, Center, Nothing}(ui.grid) +# viₙ = Field{Center, Face, Nothing}(vi.grid) + +# so = Field(sqrt(uoₙ^2 + voₙ^2)) +# si = Field(sqrt(uiₙ^2 + viₙ^2)) + +# soₙ = @lift begin +# parent(uoₙ) .= parent(uo[$n]) +# parent(voₙ) .= parent(vo[$n]) +# compute!(so) +# soₙ = interior(so) +# soₙ[land] .= NaN +# view(soₙ, :, :, 1) +# end + +# siₙ = @lift begin +# parent(uiₙ) .= parent(ui[$n]) +# parent(viₙ) .= parent(vi[$n]) +# compute!(si) +# siₙ = interior(si) +# hₙ = interior(hi[$n]) +# ℵₙ = interior(ℵi[$n]) +# he = hₙ .* ℵₙ +# siₙ[he .< 1e-7] .= 0 +# siₙ[land] .= NaN +# view(siₙ, :, :, 1) +# end + +# # Finally, we plot a snapshot of the surface speed, temperature, and the turbulent +# # eddy kinetic energy from the CATKE vertical mixing parameterization as well as the +# # sea ice speed and the effective sea ice thickness. +# fig = Figure(size=(2000, 1000)) + +# title = @lift string("Global 1/8ᵒ ocean simulation after ", prettytime(times[$n] - times[1])) + +# axso = Axis(fig[1, 1]) +# axsi = Axis(fig[1, 3]) +# axTo = Axis(fig[2, 1]) +# axhi = Axis(fig[2, 3]) + +# hmo = heatmap!(axso, soₙ, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) +# hmi = heatmap!(axsi, siₙ, colorrange = (0, 0.5), colormap = :greys, nan_color=:lightgray) +# Colorbar(fig[1, 2], hmo, label = "Ocean Surface speed (m s⁻¹)") +# Colorbar(fig[1, 4], hmi, label = "Sea ice speed (m s⁻¹)") + +# hmo = heatmap!(axTo, Toₙ, colorrange = (-1, 32), colormap = :magma, nan_color=:lightgray) +# hmi = heatmap!(axhi, heₙ, colorrange = (0, 4), colormap = :blues, nan_color=:lightgray) +# Colorbar(fig[2, 2], hmo, label = "Surface Temperature (ᵒC)") +# Colorbar(fig[2, 4], hmi, label = "Effective ice thickness (m)") + +# for ax in (axso, axsi, axTo, axhi) +# hidedecorations!(ax) +# end + +# Label(fig[0, :], title) + +# save("$(FILE_DIR)/global_snapshot.png", fig) +# nothing #hide + +# # ![](global_snapshot.png) + +# # And now a movie: + +# CairoMakie.record(fig, "$(FILE_DIR)/oneeighth_degree_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn +# n[] = nn +# end +# nothing +# #%% \ No newline at end of file diff --git a/experiments/omip_prototype/oneeighth_degree_simulation_minimal.jl b/experiments/omip_prototype/oneeighth_degree_simulation_minimal.jl new file mode 100644 index 000000000..818dbf2a9 --- /dev/null +++ b/experiments/omip_prototype/oneeighth_degree_simulation_minimal.jl @@ -0,0 +1,166 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.Architectures: on_architecture +using Oceananigans.DistributedComputations: Equal +using Oceananigans.BoundaryConditions: fill_halo_regions! +using ClimaOcean.OceanSimulations +using ClimaOcean.DataWrangling +using ClimaOcean.DataWrangling: NearestNeighborInpainting +using Printf +using Dates +using CUDA +using PythonCall +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency + +import Oceananigans.OutputWriters: checkpointer_address + +@info "GLORYS: Running one-eighth degree simulation" +# arch = GPU() +# arch = Distributed(GPU(), partition=Partition(1, 4), synchronized_communication=true) +arch = Distributed(CPU(), partition=Partition(1, 4), synchronized_communication=true) + +@info "Architecture $(arch)" + +Nx = 2880 # longitudinal direction +Ny = 1440 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0) + +const z_surf = z_faces(Nz) + +@info "Building grid..." +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +@info "Regridding bathymetry..." +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=10) + +fitted_bottom = GridFittedBottom(bottom_height) + +@info "Building immersed boundary grid..." +grid = ImmersedBoundaryGrid(grid, fitted_bottom; active_cells_map=true) + +@info grid +@info "Created ImmersedBoundaryGrid" + +##### +##### A Prognostic Ocean model +##### +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; substeps = 70) +@info "Free surface", free_surface + +obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() # CATKE +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-4)) + +glorys_dir = joinpath(homedir(), "GLORYS_data") +mkpath(glorys_dir) + +glorys_dataset = GLORYSMonthly() + +@info "Building ocean component..." +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + closure) + +start_date = DateTime(1993, 1, 1) +end_date = start_date + Month(3) +simulation_period = Dates.value(Second(end_date - start_date)) + +inpainting = NearestNeighborInpainting(50) +@info "Setting initial conditions..." + +Tᵢ = Metadatum(:temperature; dataset=glorys_dataset, date=start_date, dir=glorys_dir) +Sᵢ = Metadatum(:salinity; dataset=glorys_dataset, date=start_date, dir=glorys_dir) + +set!(ocean.model.tracers.T, Tᵢ; inpainting) +set!(ocean.model.tracers.S, Sᵢ; inpainting) + +fill_halo_regions!(ocean.model.tracers.T) +fill_halo_regions!(ocean.model.tracers.S) + +##### +##### A Prognostic Sea-ice model +##### + +@info "Building sea-ice component..." +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) + +@info "Setting sea-ice initial conditions..." +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=glorys_dataset, dir=glorys_dir, date=start_date), + ℵ=Metadatum(:sea_ice_concentration; dataset=glorys_dataset, dir=glorys_dir, date=start_date)) + +##### +##### A Prescribed Atmosphere model +##### + +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +jra55_dataset = MultiYearJRA55() +jra55_backend = JRA55NetCDFBackend(3) + +@info "Building atmospheric forcing..." +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset=jra55_dataset, backend=jra55_backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +##### +##### An ocean-sea ice coupled model +##### + +@info "Building coupled ocean-sea ice model..." +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built ocean-sea ice model" +omip = Simulation(omip, Δt=5minutes, stop_time=60days) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +add_callback!(omip, progress, IterationInterval(1)) + +@info "Starting simulation..." +run!(omip) + +# @info "Initialization complete. Running the rest..." +# omip.Δt = 10minutes +# omip.stop_time = simulation_period + +# run!(omip) \ No newline at end of file diff --git a/experiments/omip_prototype/sixth_degree_omip.jl b/experiments/omip_prototype/sixth_degree_omip.jl new file mode 100644 index 000000000..c1339d82c --- /dev/null +++ b/experiments/omip_prototype/sixth_degree_omip.jl @@ -0,0 +1,153 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA + +import Oceananigans.OutputWriters: checkpointer_address + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) + +arch = GPU() + +Nx = 2160 # longitudinal direction +Ny = 1080 # meridional direction +Nz = 60 + +r_faces = ClimaOcean.ExponentialCoordinate(Nz, -6000) +z_faces = MutableVerticalDiscretization(r_faces) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=15) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7, fixed_Δt=10minutes) + +mixing_length = CATKEMixingLength(Cᵇ=0.01) +turbulent_kinetic_energy_equation = CATKEEquation(Cᵂϵ=1.0) + +catke_closure = CATKEVerticalDiffusivity(; mixing_length, turbulent_kinetic_energy_equation) +closure = (catke_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-5)) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + free_surface, + closure) + +dataset = ECCO4Monthly() + +set!(ocean.model, T=Metadatum(:temperature; dataset), + S=Metadatum(:salinity; dataset)) + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +sea_ice = sea_ice_simulation(grid; advection=WENO(order=7)) + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset), + ℵ=Metadatum(:sea_ice_concentration; dataset)) + +##### +##### A Prescribed Atmosphere model +##### + +dir = "./forcing_data" +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +atmosphere = JRA55PrescribedAtmosphere(arch; dir, dataset, backend, include_rivers_and_icebergs=true) +radiation = Radiation() + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +omip = Simulation(omip, Δt=20, stop_time=60days) + +# Figure out the outputs.... + +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = IterationInterval(10000), + prefix = "ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = IterationInterval(10000), + prefix = "sea_ice_checkpoint", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(50)) + +run!(omip) + +omip.Δt = 6minutes +omip.stop_time = 58 * 365days + +run!(omip) diff --git a/experiments/omip_prototype/sixth_degree_omip_nodynamics.jl b/experiments/omip_prototype/sixth_degree_omip_nodynamics.jl new file mode 100644 index 000000000..6e6aa7a56 --- /dev/null +++ b/experiments/omip_prototype/sixth_degree_omip_nodynamics.jl @@ -0,0 +1,233 @@ +using Pkg +Pkg.update() +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Grids: AbstractGrid +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using Oceananigans.Fields: ConstantField +using Oceananigans.Operators + +import Oceananigans.OutputWriters: checkpointer_address +import Oceananigans.OutputWriters: saveproperty! + +saveproperty!(file, address, p::ConstantField) = file[address] = p + +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) + +arch = GPU() +r_faces = ClimaOcean.exponential_z_faces(; Nz=60, depth=6200) +z_faces = MutableVerticalDiscretization(r_faces) + +Nx = 2160 # longitudinal direction +Ny = 1080 # meridional direction +Nz = length(r_faces) - 1 + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=10) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; cfl=0.7, fixed_Δt=20minutes) + +# closure = (catke_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-5)) +# closure = (ClimaOcean.OceanSimulations.default_ocean_closure(), VerticalScalarDiffusivity(κ=1e-5, ν=1e-5)) + +closure = RiBasedVerticalDiffusivity() + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + closure) + +dataset = ECCO4Monthly() + +set!(ocean.model, T=Metadatum(:temperature; dataset), + S=Metadatum(:salinity; dataset)) + +##### +##### A Prognostic Sea-ice model +##### + +sea_ice = sea_ice_simulation(grid, ocean; dynamics = nothing) #advection=WENO(order=7)) + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset), + ℵ=Metadatum(:sea_ice_concentration; dataset)) + +# using JLD2 +# file = jldopen("sea_ice_checkpoint_iteration10000.jld2") +# parent(sea_ice.model.ice_thickness) .= CuArray(file["SeaIceModel/h/data"]) +# parent(sea_ice.model.ice_concentration) .= CuArray(file["SeaIceModel/ℵ/data"]) +# parent(sea_ice.model.velocities.u) .= CuArray(file["SeaIceModel/u/data"]) +# parent(sea_ice.model.velocities.v) .= CuArray(file["SeaIceModel/v/data"]) + +##### +##### A Prescribed Atmosphere model +##### + +dir = "./forcing_data" +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(40) + +atmosphere = JRA55PrescribedAtmosphere(arch; dir, dataset, backend, include_rivers_and_icebergs=true) +radiation = Radiation(sea_ice_albedo=0.7) + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +omip = Simulation(omip, Δt=10, stop_time=Inf) + +synch!(sea_ice.model, ocean.model) +synch!(omip.model, ocean.model) + +# Figure out the outputs.... +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = IterationInterval(10000), + prefix = "ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = IterationInterval(10000), + prefix = "sea_ice_checkpoint", + overwrite_existing = true) + +uo = ocean.model.velocities.u +vo = ocean.model.velocities.v +wo = ocean.model.velocities.w +T, S = ocean.model.tracers +η = ocean.model.free_surface.η + +ui = sea_ice.model.velocities.u +vi = sea_ice.model.velocities.v +h = sea_ice.model.ice_thickness +ℵ = sea_ice.model.ice_concentration + +omip.output_writers[:ocean_free_surface] = JLD2Writer(ocean.model, (; η), + schedule = TimeInterval(1days), + filename = "ocean_free_surface", + overwrite_existing = true) + +omip.output_writers[:ocean_surface_fields] = JLD2Writer(ocean.model, (; uo, vo, wo, T, S), + schedule = TimeInterval(1days), + filename = "ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +omip.output_writers[:sea_ice_surface_fields] = JLD2Writer(sea_ice.model, (; ui, vi, h, ℵ), + schedule = TimeInterval(1days), + filename = "sea_ice_surface_fields", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(abs, ocean.model.velocities.u) + vmax = maximum(abs, ocean.model.velocities.v) + uimax = maximum(abs, sea_ice.model.velocities.u) + vimax = maximum(abs, sea_ice.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(50)) +wizard = TimeStepWizard(cfl=0.7, max_change=1.1, max_Δt=20minutes) + +function sea_ice_cell_advection_timescale(grid, velocities) + u, v = velocities + τ = KernelFunctionOperation{Center, Center, Center}(cell_advection_timescaleᶜᶜ, grid, u, v) + return minimum(τ) +end + +@inline _inverse_timescale(i, j, k, Δ⁻¹, U, topo) = @inbounds abs(U[i, j, k]) * Δ⁻¹ +@inline _inverse_timescale(i, j, k, Δ⁻¹, U, topo::Flat) = 0 + +@inline function cell_advection_timescaleᶜᶜ(i, j, k, grid::AbstractGrid{FT, TX, TY}, u, v) where {FT, TX, TY} + Δx⁻¹ = Δx⁻¹ᶠᶜᶜ(i, j, k, grid) + Δy⁻¹ = Δy⁻¹ᶜᶠᶜ(i, j, k, grid) + + inverse_timescale_x = _inverse_timescale(i, j, k, Δx⁻¹, u, TX()) + inverse_timescale_y = _inverse_timescale(i, j, k, Δy⁻¹, v, TY()) + + inverse_timescale = inverse_timescale_x + inverse_timescale_y + + return 1 / inverse_timescale +end + +function add_wizard!(sim) + wizard(sim.model.ocean) + sea_ice = sim.model.sea_ice + Δti = 0.15 * sea_ice_cell_advection_timescale(sea_ice.model.grid, sea_ice.model.velocities) + @info "Wizard says: ocean Δt: $(ocean.Δt), sea ice Δt: $(Δti)" + iter = sea_ice.model.clock.iteration + Δtw = min(ocean.Δt, Δti) + + if iter < 5000 + sim.Δt = 60 + else + sim.Δt = Δtw + end + + @info "Final Δt is $(sim.Δt)" +end + +omip.callbacks[:wizard] = Callback(add_wizard!, IterationInterval(10)) + +run!(omip) diff --git a/half_degree_omip.jl b/half_degree_omip.jl new file mode 100644 index 000000000..6a8723594 --- /dev/null +++ b/half_degree_omip.jl @@ -0,0 +1,267 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using JLD2 +using ArgParse + +import Oceananigans.OutputWriters: checkpointer_address + +using Libdl +ucx_libs = filter(lib -> occursin("ucx", lowercase(lib)), Libdl.dllist()) +if isempty(ucx_libs) + @info "✓ No UCX - safe to run!" +else + @warn "✗ UCX libraries detected! This can cause issues with MPI+CUDA. Detected libs:\n$(join(ucx_libs, "\n"))" +end + +function parse_commandline() + s = ArgParseSettings() + + @add_arg_table! s begin + "--kappa_skew" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--kappa_symmetric" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + end + return parse_args(s) +end + +args = parse_commandline() +κ_skew = args["kappa_skew"] +κ_symmetric = args["kappa_symmetric"] + +CUDA.versioninfo() + +@info "Using κ_skew = $(κ_skew) m²/s and κ_symmetric = $(κ_symmetric) m²/s" + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) +# restart_iteration = "446000" +restart_iteration = nothing + +arch = GPU() + +Nx = 720 # longitudinal direction +Ny = 360 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0; scale=1800) +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=55) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization, AdvectiveFormulation, IsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=7) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8, fixed_Δt=50minutes) + +using Oceananigans.Operators: Δx, Δy + +@inline Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz) = 2 * (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2)) +@inline geometric_νhb(i, j, k, grid, lx, ly, lz, clock, fields, λ) = Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz)^2 / λ + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, κ_symmetric, skew_flux_formulation=AdvectiveFormulation()) +# obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +obl_closure = RiBasedVerticalDiffusivity() +visc_closure = HorizontalScalarBiharmonicDiffusivity(ν=geometric_νhb, discrete_form=true, parameters=25days) + +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=3e-4), visc_closure, eddy_closure) + +prefix = "halfdegree" +if obl_closure isa RiBasedVerticalDiffusivity + prefix *= "_RiBased" +else + prefix *= "_CATKE" +end + +prefix *= "_$(κ_skew)_$(κ_symmetric)" +prefix *= "_advectiveGM_multiyearjra55" + +dir = joinpath(homedir(), "forcing_data_half_degree") +mkpath(dir) + +start_date = DateTime(1958, 1, 1) +end_date = start_date + Year(20) +simulation_period = Dates.value(Second(end_date - start_date)) +yearly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(1):end_date))))) +decadal_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(10):end_date))))) + +@info "Settting up salinity restoring..." +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset=EN4Monthly(), dir, start_date, end_date) +FS = DatasetRestoring(Smetadata, grid; rate = 1/30days, mask, time_indices_in_memory = 10) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +@info "Built ocean model $(ocean)" + +set!(ocean.model, T=Metadatum(:temperature; dataset=EN4Monthly(), date=start_date, dir), + S=Metadatum(:salinity; dataset=EN4Monthly(), date=start_date, dir)) +@info "Initialized T and S" + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) +@info "Built sea ice model $(sea_ice)" + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir)) + +@info "Initialized sea ice fields" + +##### +##### A Prescribed Atmosphere model +##### +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +@info "Setting up presctibed atmosphere $(dataset)" +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +@info "Built atmosphere model $(atmosphere)" + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built coupled model $(omip)" + +omip = Simulation(omip, Δt=30minutes, stop_time=simulation_period) +@info "Built simulation $(omip)" + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1825days), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:time_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(1825days, window=1825days), + filename = "$(FILE_DIR)/ocean_complete_fields_10year_average", + overwrite_existing = true) + +sea_ice.output_writers[:time_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(1825days, window=1825days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_10year_average", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(100)) + +run!(omip) \ No newline at end of file diff --git a/half_degree_omip_calibrationsamples.jl b/half_degree_omip_calibrationsamples.jl new file mode 100644 index 000000000..1e61d09d2 --- /dev/null +++ b/half_degree_omip_calibrationsamples.jl @@ -0,0 +1,291 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using JLD2 +using ArgParse +using Oceananigans.OutputWriters: AveragedSpecifiedTimes + +import Oceananigans.OutputWriters: checkpointer_address + +using Libdl +ucx_libs = filter(lib -> occursin("ucx", lowercase(lib)), Libdl.dllist()) +if isempty(ucx_libs) + @info "✓ No UCX - safe to run!" +else + @warn "✗ UCX libraries detected! This can cause issues with MPI+CUDA. Detected libs:\n$(join(ucx_libs, "\n"))" +end + +function parse_commandline() + s = ArgParseSettings() + + @add_arg_table! s begin + "--kappa_skew" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--kappa_symmetric" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--start_year" + help = "Start year of the simulation" + arg_type = Int + default = 1962 + "--simulation_length" + help = "Length of the simulation in years" + arg_type = Int + default = 20 + "--sampling_length" + help = "Length of the sampling window in years" + arg_type = Int + default = 10 + end + return parse_args(s) +end + +args = parse_commandline() +κ_skew = args["kappa_skew"] +κ_symmetric = args["kappa_symmetric"] +start_year = args["start_year"] +simulation_length = args["simulation_length"] +sampling_length = args["sampling_length"] + +@info "Using κ_skew = $(κ_skew) m²/s and κ_symmetric = $(κ_symmetric) m²/s, starting in year $(start_year) for a length of $(simulation_length) years with a $(sampling_length)-year sample." + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) +# restart_iteration = "446000" +restart_iteration = nothing + +arch = GPU() + +Nx = 720 # longitudinal direction +Ny = 360 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0; scale=1800) +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=55) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization, AdvectiveFormulation, IsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=7) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8, fixed_Δt=50minutes) + +using Oceananigans.Operators: Δx, Δy + +@inline Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz) = 2 * (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2)) +@inline geometric_νhb(i, j, k, grid, lx, ly, lz, clock, fields, λ) = Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz)^2 / λ + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, κ_symmetric, skew_flux_formulation=AdvectiveFormulation()) +# obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +obl_closure = RiBasedVerticalDiffusivity() +visc_closure = HorizontalScalarBiharmonicDiffusivity(ν=geometric_νhb, discrete_form=true, parameters=25days) + +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=3e-4), visc_closure, eddy_closure) + +prefix = "halfdegree" +if obl_closure isa RiBasedVerticalDiffusivity + prefix *= "_RiBased" +else + prefix *= "_CATKE" +end + +prefix *= "_$(κ_skew)_$(κ_symmetric)" +prefix *= "_$(start_year)" +prefix *= "_$(simulation_length)year_$(sampling_length)yearsample" +prefix *= "_advectiveGM_multiyearjra55_calibrationsamples" + +dir = joinpath(homedir(), "forcing_data_half_degree") +mkpath(dir) + +start_date = DateTime(start_year, 1, 1) +end_date = start_date + Year(simulation_length) +simulation_period = Dates.value(Second(end_date - start_date)) +yearly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(1):end_date))))) +decadal_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(10):end_date))))) +# sampling_endtimes = decadal_times[3:end] +sampling_start_date = end_date - Year(sampling_length) +sampling_window = Dates.value(Second(end_date - sampling_start_date)) + +@info "Settting up salinity restoring..." +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset=EN4Monthly(), dir, start_date, end_date) +FS = DatasetRestoring(Smetadata, grid; rate = 1/30days, mask, time_indices_in_memory = 10) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +@info "Built ocean model $(ocean)" + +set!(ocean.model, T=Metadatum(:temperature; dataset=EN4Monthly(), date=start_date, dir), + S=Metadatum(:salinity; dataset=EN4Monthly(), date=start_date, dir)) +@info "Initialized T and S" + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) +@info "Built sea ice model $(sea_ice)" + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir)) + +@info "Initialized sea ice fields" + +##### +##### A Prescribed Atmosphere model +##### +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +@info "Setting up presctibed atmosphere $(dataset)" +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +@info "Built atmosphere model $(atmosphere)" + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built coupled model $(omip)" + +omip = Simulation(omip, Δt=30minutes, stop_time=simulation_period) +@info "Built simulation $(omip)" + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1825days), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:time_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/ocean_complete_fields_10year_average", + overwrite_existing = true) + +sea_ice.output_writers[:time_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_10year_average", + overwrite_existing = true) + +ocean.output_writers[:sample_decadal_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(simulation_period, window=sampling_window), + filename = "$(FILE_DIR)/ocean_complete_fields_$(sampling_length)year_average_calibrationsample", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(100)) + +run!(omip) \ No newline at end of file diff --git a/half_degree_omip_calibrationsamples_eccoinitialization.jl b/half_degree_omip_calibrationsamples_eccoinitialization.jl new file mode 100644 index 000000000..02c96a10a --- /dev/null +++ b/half_degree_omip_calibrationsamples_eccoinitialization.jl @@ -0,0 +1,290 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using JLD2 +using ArgParse +using Oceananigans.OutputWriters: AveragedSpecifiedTimes + +import Oceananigans.OutputWriters: checkpointer_address + +using Libdl +ucx_libs = filter(lib -> occursin("ucx", lowercase(lib)), Libdl.dllist()) +if isempty(ucx_libs) + @info "✓ No UCX - safe to run!" +else + @warn "✗ UCX libraries detected! This can cause issues with MPI+CUDA. Detected libs:\n$(join(ucx_libs, "\n"))" +end + +function parse_commandline() + s = ArgParseSettings() + + @add_arg_table! s begin + "--kappa_skew" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--kappa_symmetric" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 5e2 + "--start_year" + help = "Start year of the simulation" + arg_type = Int + default = 1962 + "--simulation_length" + help = "Length of the simulation in years" + arg_type = Int + default = 20 + "--sampling_length" + help = "Length of the sampling window in years" + arg_type = Int + default = 10 + end + return parse_args(s) +end + +args = parse_commandline() +κ_skew = args["kappa_skew"] +κ_symmetric = args["kappa_symmetric"] +start_year = args["start_year"] +simulation_length = args["simulation_length"] +sampling_length = args["sampling_length"] + +@info "Using κ_skew = $(κ_skew) m²/s and κ_symmetric = $(κ_symmetric) m²/s, starting in year $(start_year) for a length of $(simulation_length) years with a $(sampling_length)-year sample." + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) +# restart_iteration = "446000" +restart_iteration = nothing + +arch = GPU() + +Nx = 720 # longitudinal direction +Ny = 360 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0; scale=1800) +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=55) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization, AdvectiveFormulation, IsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=7) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8, fixed_Δt=50minutes) + +using Oceananigans.Operators: Δx, Δy + +@inline Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz) = 2 * (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2)) +@inline geometric_νhb(i, j, k, grid, lx, ly, lz, clock, fields, λ) = Δ²ᵃᵃᵃ(i, j, k, grid, lx, ly, lz)^2 / λ + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, κ_symmetric, skew_flux_formulation=AdvectiveFormulation()) +obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +visc_closure = HorizontalScalarBiharmonicDiffusivity(ν=geometric_νhb, discrete_form=true, parameters=25days) + +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=3e-4), visc_closure, eddy_closure) + +prefix = "halfdegree_eccoinitial" +if obl_closure isa RiBasedVerticalDiffusivity + prefix *= "_RiBased" +else + prefix *= "_CATKE" +end + +prefix *= "_$(κ_skew)_$(κ_symmetric)" +prefix *= "_$(start_year)" +prefix *= "_$(simulation_length)year_$(sampling_length)yearsample" +prefix *= "_advectiveGM_multiyearjra55_calibrationsamples" + +dir = joinpath(homedir(), "ECCO_data") +mkpath(dir) + +start_date = DateTime(start_year, 1, 1) +end_date = start_date + Year(simulation_length) +simulation_period = Dates.value(Second(end_date - start_date)) +yearly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(1):end_date))))) +decadal_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(10):end_date))))) +# sampling_endtimes = decadal_times[3:end] +sampling_start_date = end_date - Year(sampling_length) +sampling_window = Dates.value(Second(end_date - sampling_start_date)) + +@info "Settting up salinity restoring..." +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset=ECCO4Monthly(), dir, start_date, end_date) +FS = DatasetRestoring(Smetadata, grid; rate = 1/30days, mask, time_indices_in_memory = 10) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +@info "Built ocean model $(ocean)" + +set!(ocean.model, T=Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, dir), + S=Metadatum(:salinity; dataset=ECCO4Monthly(), date=start_date, dir)) +@info "Initialized T and S" + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) +@info "Built sea ice model $(sea_ice)" + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), date=start_date, dir), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), date=start_date, dir)) + +@info "Initialized sea ice fields" + +##### +##### A Prescribed Atmosphere model +##### +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +@info "Setting up presctibed atmosphere $(dataset)" +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +@info "Built atmosphere model $(atmosphere)" + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built coupled model $(omip)" + +omip = Simulation(omip, Δt=30minutes, stop_time=simulation_period) +@info "Built simulation $(omip)" + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1825days), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:time_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/ocean_complete_fields_10year_average", + overwrite_existing = true) + +sea_ice.output_writers[:time_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_10year_average", + overwrite_existing = true) + +ocean.output_writers[:sample_decadal_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(simulation_period, window=sampling_window), + filename = "$(FILE_DIR)/ocean_complete_fields_$(sampling_length)year_average_calibrationsample", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(100)) + +run!(omip) \ No newline at end of file diff --git a/one_degree_omip_calibrationsamples.jl b/one_degree_omip_calibrationsamples.jl new file mode 100644 index 000000000..d6722022e --- /dev/null +++ b/one_degree_omip_calibrationsamples.jl @@ -0,0 +1,288 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using Oceananigans.BuoyancyFormulations: buoyancy, buoyancy_frequency +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA +using JLD2 +using ArgParse +using Oceananigans.OutputWriters: AveragedSpecifiedTimes + +import Oceananigans.OutputWriters: checkpointer_address + +using Libdl +ucx_libs = filter(lib -> occursin("ucx", lowercase(lib)), Libdl.dllist()) +if isempty(ucx_libs) + @info "✓ No UCX - safe to run!" +else + @warn "✗ UCX libraries detected! This can cause issues with MPI+CUDA. Detected libs:\n$(join(ucx_libs, "\n"))" +end + +function parse_commandline() + s = ArgParseSettings() + + @add_arg_table! s begin + "--kappa_skew" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 1e3 + "--kappa_symmetric" + help = "Isopycnal skew diffusivity (m^2/s)" + arg_type = Float64 + default = 1e3 + "--start_year" + help = "Start year of the simulation" + arg_type = Int + default = 1992 + "--simulation_length" + help = "Length of the simulation in years" + arg_type = Int + default = 25 + "--sampling_length" + help = "Length of the sampling window in years" + arg_type = Int + default = 10 + end + return parse_args(s) +end + +args = parse_commandline() +κ_skew = args["kappa_skew"] +κ_symmetric = args["kappa_symmetric"] +start_year = args["start_year"] +simulation_length = args["simulation_length"] +sampling_length = args["sampling_length"] + +@info "1-degree omip" +@info "Using κ_skew = $(κ_skew) m²/s and κ_symmetric = $(κ_symmetric) m²/s, starting in year $(start_year) for a length of $(simulation_length) years with a $(sampling_length)-year sample." + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) +# restart_iteration = "446000" +restart_iteration = nothing + +arch = GPU() + +Nx = 360 # longitudinal direction +Ny = 180 # meridional direction +Nz = 100 + +z_faces = ExponentialDiscretization(Nz, -6000, 0; scale=1800) +const z_surf = z_faces(Nz) + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1, interpolation_passes=75) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization, AdvectiveFormulation, IsopycnalSkewSymmetricDiffusivity +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=5) +free_surface = SplitExplicitFreeSurface(grid; cfl=0.8, fixed_Δt=50minutes) + +using Oceananigans.Operators: Δx, Δy + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, κ_symmetric, skew_flux_formulation=AdvectiveFormulation()) +# obl_closure = ClimaOcean.OceanSimulations.default_ocean_closure() +obl_closure = RiBasedVerticalDiffusivity() + +closure = (obl_closure, VerticalScalarDiffusivity(κ=1e-5, ν=3e-4), eddy_closure) + +prefix = "onedegree" +if obl_closure isa RiBasedVerticalDiffusivity + prefix *= "_RiBased" +else + prefix *= "_CATKE" +end + +prefix *= "_$(κ_skew)_$(κ_symmetric)" +prefix *= "_$(start_year)" +prefix *= "_$(simulation_length)year_$(sampling_length)yearsample" +prefix *= "_advectiveGM_multiyearjra55_calibrationsamples" + +dir = joinpath(homedir(), "ECCO_data") +mkpath(dir) + +start_date = DateTime(start_year, 1, 1) +end_date = start_date + Year(simulation_length) +simulation_period = Dates.value(Second(end_date - start_date)) +yearly_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(1):end_date))))) +decadal_times = cumsum(vcat([0.], Dates.value.(Second.(diff(start_date:Year(10):end_date))))) +# sampling_endtimes = decadal_times[3:end] +sampling_start_date = end_date - Year(sampling_length) +sampling_window = Dates.value(Second(end_date - sampling_start_date)) + +@info "Settting up salinity restoring..." +@inline mask(x, y, z, t) = z ≥ z_surf - 1 +Smetadata = Metadata(:salinity; dataset=ECCO4Monthly(), dir, start_date, end_date) +FS = DatasetRestoring(Smetadata, grid; rate = 1/18days, mask, time_indices_in_memory = 10) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + timestepper = :SplitRungeKutta3, + free_surface, + forcing = (; S = FS), + closure) + +@info "Built ocean model $(ocean)" + +set!(ocean.model, T=Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, dir), + S=Metadatum(:salinity; dataset=ECCO4Monthly(), date=start_date, dir)) +@info "Initialized T and S" + +##### +##### A Prognostic Sea-ice model +##### + +# Default sea-ice dynamics and salinity coupling are included in the defaults +# sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) +sea_ice = sea_ice_simulation(grid, ocean; dynamics=nothing) +@info "Built sea ice model $(sea_ice)" + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset=ECCO4Monthly(), dir), + ℵ=Metadatum(:sea_ice_concentration; dataset=ECCO4Monthly(), dir)) + +@info "Initialized sea ice fields" + +##### +##### A Prescribed Atmosphere model +##### +jra55_dir = joinpath(homedir(), "JRA55_data") +mkpath(jra55_dir) +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(100) + +@info "Setting up presctibed atmosphere $(dataset)" +atmosphere = JRA55PrescribedAtmosphere(arch; dir=jra55_dir, dataset, backend, include_rivers_and_icebergs=true, start_date, end_date) +radiation = Radiation() + +@info "Built atmosphere model $(atmosphere)" + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +@info "Built coupled model $(omip)" + +omip = Simulation(omip, Δt=40minutes, stop_time=simulation_period) +@info "Built simulation $(omip)" + +# Figure out the outputs.... +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +FILE_DIR = "./Data/$(prefix)/" +mkpath(FILE_DIR) + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = SpecifiedTimes(decadal_times), + prefix = "$(FILE_DIR)/sea_ice_checkpoint", + overwrite_existing = true) + +u, v, w = ocean.model.velocities +T, S = ocean.model.tracers +b = Field(buoyancy(ocean.model)) +N² = Field(buoyancy_frequency(ocean.model)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities, (; b, N²)) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/ocean_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs; + schedule = TimeInterval(180days), + filename = "$(FILE_DIR)/sea_ice_surface_fields", + overwrite_existing = true) + +ocean.output_writers[:full] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1825days), + filename = "$(FILE_DIR)/ocean_complete_fields", + overwrite_existing = true) + +ocean.output_writers[:time_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/ocean_complete_fields_10year_average", + overwrite_existing = true) + +sea_ice.output_writers[:time_average] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = AveragedTimeInterval(3650days, window=3650days), + filename = "$(FILE_DIR)/sea_ice_complete_fields_10year_average", + overwrite_existing = true) + +ocean.output_writers[:sample_decadal_average] = JLD2Writer(ocean.model, ocean_outputs; + schedule = AveragedTimeInterval(simulation_period, window=sampling_window), + filename = "$(FILE_DIR)/ocean_complete_fields_$(sampling_length)year_average_calibrationsample", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(100)) + +run!(omip) \ No newline at end of file diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index a647ebe37..e54ac2594 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -101,7 +101,7 @@ function _regrid_bathymetry(target_grid, metadata; arch = architecture(target_grid) - bathymetry_native_grid = native_grid(metadata, arch; halo = (10, 10, 1)) + bathymetry_native_grid = native_grid(metadata, arch; halo = (20, 20, 1)) FT = eltype(target_grid) filepath = metadata_path(metadata) @@ -195,7 +195,8 @@ end # Fix active cells to be at least `-minimum_depth`. active = z < 0 # it's a wet cell - z = ifelse(active, min(z, -minimum_depth), z) + above_minimum_depth = z > -minimum_depth + z = ifelse(active, ifelse(above_minimum_depth, zero(z), z), z) @inbounds target_z[i, j, 1] = z end diff --git a/src/DataWrangling/Copernicus/Copernicus.jl b/src/DataWrangling/Copernicus/Copernicus.jl index 2ee7f0116..8c7ac0dce 100644 --- a/src/DataWrangling/Copernicus/Copernicus.jl +++ b/src/DataWrangling/Copernicus/Copernicus.jl @@ -22,7 +22,8 @@ import ClimaOcean.DataWrangling: metadata_filename, inpainted_metadata_path, reversed_vertical_axis, - available_variables + available_variables, + is_three_dimensional using Scratch @@ -129,20 +130,31 @@ end inpainted_metadata_path(metadata::CopernicusMetadata) = joinpath(metadata.dir, inpainted_metadata_filename(metadata)) -location(::CopernicusMetadata) = (Center, Center, Center) +location(metadata::CopernicusMetadata) = is_three_dimensional(metadata) ? (Center, Center, Center) : (Center, Center, Nothing) longitude_interfaces(::CopernicusMetadata) = (0, 360) latitude_interfaces(::CopernicusMetadata) = (-80, 90) +is_three_dimensional(metadata::CopernicusMetadata) = + metadata.name == :temperature || + metadata.name == :salinity || + metadata.name == :u_velocity || + metadata.name == :v_velocity + function z_interfaces(metadata::CopernicusMetadata) - path = metadata_path(metadata) - ds = Dataset(path) - zc = - reverse(ds["depth"][:]) - close(ds) - dz = zc[2] - zc[1] - zf = zc[1:end-1] .+ zc[2:end] - push!(zf, 0) - pushfirst!(zf, zf[1] - dz) - return zf + # Check if this is a 2D surface variable + if !is_three_dimensional(metadata) + return (0, 1) + else + path = metadata_path(metadata) + ds = Dataset(path) + zc = - reverse(ds["depth"][:]) + close(ds) + dz = zc[2] - zc[1] + zf = zc[1:end-1] .+ zc[2:end] + push!(zf, 0) + pushfirst!(zf, zf[1] - dz) + return zf + end end -end # module Copernicus +end # module Copernicus \ No newline at end of file diff --git a/src/DataWrangling/DataWrangling.jl b/src/DataWrangling/DataWrangling.jl index e5da9795f..ca0d7ced2 100644 --- a/src/DataWrangling/DataWrangling.jl +++ b/src/DataWrangling/DataWrangling.jl @@ -206,7 +206,7 @@ function default_inpainting(metadata) if metadata.name in (:temperature, :salinity) return NearestNeighborInpainting(Inf) elseif metadata.name in (:sea_ice_thickness, :sea_ice_concentration) - return nothing + return ValueInpainting(0) else return NearestNeighborInpainting(5) end diff --git a/src/DataWrangling/JRA55/JRA55_field_time_series.jl b/src/DataWrangling/JRA55/JRA55_field_time_series.jl index 2def6129f..3a5ca820f 100644 --- a/src/DataWrangling/JRA55/JRA55_field_time_series.jl +++ b/src/DataWrangling/JRA55/JRA55_field_time_series.jl @@ -189,7 +189,6 @@ function set!(fts::JRA55NetCDFFTSMultipleYears, backend=fts.backend) LX, LY, LZ = location(fts) i₁, i₂, j₁, j₂, TX = compute_bounding_indices(nothing, nothing, fts.grid, LX, LY, λc, φc) - if issorted(nn) data = ds[name][i₁:i₂, j₁:j₂, nn] else @@ -458,7 +457,7 @@ function JRA55FieldTimeSeries(metadata::JRA55Metadata, architecture=CPU(), FT=Fl ds = Dataset(filepath) data = ds[shortname][i₁:i₂, j₁:j₂, time_indices_in_memory] close(ds) - + copyto!(interior(fts, :, :, 1, :), data) fill_halo_regions!(fts) diff --git a/src/DataWrangling/inpainting.jl b/src/DataWrangling/inpainting.jl index dde0f1f91..93414de91 100644 --- a/src/DataWrangling/inpainting.jl +++ b/src/DataWrangling/inpainting.jl @@ -11,6 +11,30 @@ struct NearestNeighborInpainting{M} maxiter :: M end +""" + ValueInpainting{V} + +A structure representing a simple value inpainting algorithm, where all missing values are +substituted with a specified constant value. + +# Fields +- `value :: V`: The constant value to use for inpainting missing data. +- `maxiter`: For cache compatibility (not used in the algorithm itself). Defaults to 0. + +# Example +```julia +inpainting = ValueInpainting(0.0) # Fill missing values with 0 +inpaint_mask!(field, mask; inpainting) +``` +""" +struct ValueInpainting{V, M} + value :: V + maxiter :: M +end + +# Convenience constructor that sets maxiter to a default value +ValueInpainting(value) = ValueInpainting(value, 0) + propagate_horizontally!(field, ::Nothing, substituting_field=deepcopy(field); kw...) = field function propagating(field, mask, iter, inpainting::NearestNeighborInpainting) @@ -18,6 +42,9 @@ function propagating(field, mask, iter, inpainting::NearestNeighborInpainting) return nans > 0 && iter < inpainting.maxiter end +# ValueInpainting doesn't need iteration +propagating(field, mask, iter, inpainting::ValueInpainting) = false + """ propagate_horizontally!(inpainting, field, mask [, substituting_field=deepcopy(field)]) @@ -57,6 +84,17 @@ function propagate_horizontally!(inpainting::NearestNeighborInpainting, field, m return field end +function propagate_horizontally!(inpainting::ValueInpainting, field, mask, + substituting_field=deepcopy(field)) + grid = field.grid + arch = architecture(grid) + + launch!(arch, grid, size(field), _fill_with_number!, field, mask, inpainting.value) + fill_halo_regions!(field) + + return field +end + # Maybe we can remove this propagate field in lieu of a diffusion, # Still we'll need to do this a couple of steps on the original grid @kernel function _propagate_field!(substituting_field, ::NearestNeighborInpainting, field) @@ -102,6 +140,12 @@ end @inbounds field[i, j, k] *= !isnan(field[i, j, k]) end +@kernel function _fill_with_number!(field, mask, value) + i, j, k = @index(Global, NTuple) + FT_value = convert(eltype(field), value) + @inbounds field[i, j, k] = ifelse(mask[i, j, k], FT_value, field[i, j, k]) +end + """ inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf)) @@ -117,10 +161,24 @@ Arguments - `mask`: Boolean-valued `Field`, values where `mask[i, j, k] == true` are inpainted. -- `inpainting`: The inpainting algorithm to use. The only option is - `NearestNeighborInpainting(maxiter)`, where an average - of the valid surrounding values is used `maxiter` times. - Default: `NearestNeighborInpainting(Inf)`. +- `inpainting`: The inpainting algorithm to use. Options include: + * `NearestNeighborInpainting(maxiter)`: Uses an average + of the valid surrounding values, repeated `maxiter` times. + Default: `NearestNeighborInpainting(Inf)`. + * `ValueInpainting(value)`: Fills all missing values with + the specified constant `value`. + +# Examples +```julia +# Use nearest neighbor inpainting with default settings +inpaint_mask!(field, mask) + +# Fill missing values with zero +inpaint_mask!(field, mask; inpainting=ValueInpainting(0)) + +# Fill missing values with a specific temperature +inpaint_mask!(field, mask; inpainting=ValueInpainting(273.15)) +``` """ function inpaint_mask!(field, mask; inpainting=NearestNeighborInpainting(Inf)) diff --git a/src/DataWrangling/metadata_field.jl b/src/DataWrangling/metadata_field.jl index a57b6622c..1802585d3 100644 --- a/src/DataWrangling/metadata_field.jl +++ b/src/DataWrangling/metadata_field.jl @@ -211,6 +211,7 @@ function set_metadata_field!(field, data, metadatum) return nothing end +# TODO: sort out elegant way to handle NaNs for sea ice data in GLORYS @kernel function _set_2d_metadata_field!(field, data, mangling, temp_units, conc_units) i, j = @index(Global, NTuple) d = mangle(i, j, data, mangling) @@ -223,6 +224,8 @@ end elseif !isnothing(conc_units) d = convert_concentration(d, conc_units) end + + d = convert_temperature(d, temp_units) @inbounds field[i, j, 1] = d end @@ -241,6 +244,7 @@ end elseif !isnothing(conc_units) d = convert_concentration(d, conc_units) end + @inbounds field[i, j, k] = d end diff --git a/src/InitialConditions/InitialConditions.jl b/src/InitialConditions/InitialConditions.jl index 6a815f056..a529e6046 100644 --- a/src/InitialConditions/InitialConditions.jl +++ b/src/InitialConditions/InitialConditions.jl @@ -18,7 +18,7 @@ using JLD2 # Implementation of 3-dimensional regridding # TODO: move all the following to Oceananigans! -using Oceananigans.Fields: regrid!, interpolate! +using Oceananigans.Fields: interpolate! using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z, @@ -31,44 +31,6 @@ construct_grid(::Type{<:RectilinearGrid}, arch, size, extent, topology) = construct_grid(::Type{<:LatitudeLongitudeGrid}, arch, size, extent, topology) = LatitudeLongitudeGrid(arch; size, longitude = extent[1], latitude = extent[2], z = extent[3], topology) -# Regrid a field in three dimensions -function three_dimensional_regrid!(a, b) - target_grid = a.grid isa ImmersedBoundaryGrid ? a.grid.underlying_grid : a.grid - source_grid = b.grid isa ImmersedBoundaryGrid ? b.grid.underlying_grid : b.grid - - topo = topology(target_grid) - arch = architecture(target_grid) - arch = child_architecture(arch) - - target_y = yt = cpu_face_constructor_y(target_grid) - target_z = zt = cpu_face_constructor_z(target_grid) - - target_size = Nt = size(target_grid) - - source_x = xs = cpu_face_constructor_x(source_grid) - source_y = ys = cpu_face_constructor_y(source_grid) - - source_size = Ns = size(source_grid) - - # Start by regridding in z - @debug "Regridding in z" - zgrid = construct_grid(typeof(target_grid), arch, (Ns[1], Ns[2], Nt[3]), (xs, ys, zt), topo) - field_z = Field(location(b), zgrid) - regrid!(field_z, zgrid, source_grid, b) - - # regrid in y - @debug "Regridding in y" - ygrid = construct_grid(typeof(target_grid), arch, (Ns[1], Nt[2], Nt[3]), (xs, yt, zt), topo) - field_y = Field(location(b), ygrid); - regrid!(field_y, ygrid, zgrid, field_z); - - # Finally regrid in x - @debug "Regridding in x" - regrid!(a, target_grid, ygrid, field_y) - - return a -end - include("diffuse_tracers.jl") end # module diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 1532a3be7..dff50271b 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -261,8 +261,8 @@ end @inbounds begin Ts = surface_temperature[i, j, kᴺ] Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) - ℵi = ice_concentration[i, j, 1] - + ℵi = ice_concentration[i, j, kᴺ] + Qs = downwelling_radiation.Qs[i, j, 1] Qℓ = downwelling_radiation.Qℓ[i, j, 1] Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux diff --git a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index d12ae69b3..9f0b3de73 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -14,6 +14,7 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model) v = ZeroField(), h = sea_ice.model.ice_thickness, ℵ = sea_ice.model.ice_concentration, + hc = sea_ice.model.ice_consolidation_thickness, Tₒ = ocean.model.tracers.T, Sₒ = ocean.model.tracers.S) @@ -97,6 +98,7 @@ end uᵢ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u) vᵢ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v) hᵢ = interior_state.h[i, j, 1] + hc = interior_state.hc[i, j, 1] ℵᵢ = interior_state.ℵ[i, j, 1] Tₛ = interface_temperature[i, j, 1] Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ) @@ -117,7 +119,7 @@ end h_bℓ = atmosphere_state.h_bℓ) downwelling_radiation = (; Qs, Qℓ) - local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ) + local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, hc=hc) # Estimate initial interface state u★ = convert(FT, 1e-4) diff --git a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl index 0ad3243ee..fd6143295 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/interface_states.jl @@ -269,9 +269,10 @@ end # where Ωc (the sensible heat transfer coefficient) is given by Ωc = Qc / (Tₐ - Tˢ) # ⟹ Tₛ = (Tᵢ * k - (Qv + Qu + Qd + Ωc * Tₐ) * h / (k - Ωc * h) @inline function flux_balance_temperature(st::SkinTemperature{<:ClimaSeaIce.ConductiveFlux}, Ψₛ, ℙₛ, Qc, Qv, Qu, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) - F = st.internal_flux - k = F.conductivity - h = Ψᵢ.h + F = st.internal_flux + k = F.conductivity + h = Ψᵢ.h + hc = Ψᵢ.hc # Critical thickness for ice consolidation # Bottom temperature at the melting temperature Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) @@ -303,6 +304,9 @@ end Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) Tₛ⁺ = min(Tₛ⁺, Tₘ) + # If the ice is not consolidated, use the bottom temperature + Tₛ⁺ = ifelse(h ≥ hc, Tₛ⁺, Tᵢ) + return Tₛ⁺ end diff --git a/src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl b/src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl index 62845744a..d7ca790a2 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl @@ -20,23 +20,25 @@ Base.show(io::IO, ::ScalarRoughnessLength{FT}) where FT = print(io, "ScalarRough struct WindDependentWaveFormulation{FT} Umax :: FT + αmin :: FT ℂ₁ :: FT ℂ₂ :: FT end """ WindDependentWaveFormulation(FT = Oceananigans.defaults.FloatType; - Umax = 19, ℂ₁ = 0.0017, ℂ₂ = -0.005) + Umax = 19, αmin = 0.011, ℂ₁ = 0.0017, ℂ₂ = -0.005) -A gravity wave parameter based on the wind speed `ΔU` with the formula `ℂ₁ * max(ΔU, Umax) + ℂ₂`. +A gravity wave parameter based on the wind speed `ΔU` with the formula `max(αmin, ℂ₁ * min(ΔU, Umax) + ℂ₂`). """ -WindDependentWaveFormulation(FT=Oceananigans.defaults.FloatType; Umax = 19, ℂ₁ = 0.0017, ℂ₂ = -0.005) = +WindDependentWaveFormulation(FT=Oceananigans.defaults.FloatType; Umax = 19, αmin = 0.011, ℂ₁ = 0.0017, ℂ₂ = -0.005) = WindDependentWaveFormulation(convert(FT, Umax), + convert(FT, αmin), convert(FT, ℂ₁), convert(FT, ℂ₂)) gravity_wave_parameter(α::Number, args...) = α -gravity_wave_parameter(α::WindDependentWaveFormulation, ΔU) = α.ℂ₁ * max(ΔU, α.Umax) + α.ℂ₂ +gravity_wave_parameter(α::WindDependentWaveFormulation, ΔU) = max(α.αmin, α.ℂ₁ * min(ΔU, α.Umax) + α.ℂ₂) """ ScalarRoughnessLength(FT = Float64; @@ -88,7 +90,7 @@ function MomentumRoughnessLength(FT=Oceananigans.defaults.FloatType; gravitational_acceleration = default_gravitational_acceleration, maximum_roughness_length = 1, air_kinematic_viscosity = 1.5e-5, - wave_formulation = 0.02, + wave_formulation = WindDependentWaveFormulation(FT), smooth_wall_parameter = 0.11) if wave_formulation isa Number diff --git a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl index 1000486ac..a1c54426c 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl @@ -23,6 +23,7 @@ struct SimilarityTheoryFluxes{FT, UF, R, B, S} von_karman_constant :: FT # parameter turbulent_prandtl_number :: FT # parameter gustiness_parameter :: FT # bulk velocity parameter + minimum_velocity_scale :: FT # minimum velocity scale stability_functions :: UF # functions for turbulent fluxes roughness_lengths :: R # parameterization for turbulent fluxes similarity_form :: B # similarity profile relating atmosphere to interface state @@ -33,6 +34,7 @@ Adapt.adapt_structure(to, fluxes::SimilarityTheoryFluxes) = SimilarityTheoryFluxes(adapt(to, fluxes.von_karman_constant), adapt(to, fluxes.turbulent_prandtl_number), adapt(to, fluxes.gustiness_parameter), + adapt(to, fluxes.minimum_velocity_scale), adapt(to, fluxes.stability_functions), adapt(to, fluxes.roughness_lengths), adapt(to, fluxes.similarity_form), @@ -86,11 +88,12 @@ Keyword Arguments function SimilarityTheoryFluxes(FT::DataType = Oceananigans.defaults.FloatType; von_karman_constant = 0.4, turbulent_prandtl_number = 1, - gustiness_parameter = 1, + gustiness_parameter = 2, stability_functions = atmosphere_ocean_stability_functions(FT), momentum_roughness_length = MomentumRoughnessLength(FT), temperature_roughness_length = ScalarRoughnessLength(FT), water_vapor_roughness_length = ScalarRoughnessLength(FT), + minimum_velocity_scale = 0.5, similarity_form = LogarithmicSimilarityProfile(), solver_stop_criteria = nothing, solver_tolerance = 1e-8, @@ -113,6 +116,7 @@ function SimilarityTheoryFluxes(FT::DataType = Oceananigans.defaults.FloatType; return SimilarityTheoryFluxes(convert(FT, von_karman_constant), convert(FT, turbulent_prandtl_number), convert(FT, gustiness_parameter), + convert(FT, minimum_velocity_scale), stability_functions, roughness_lengths, similarity_form, @@ -202,7 +206,8 @@ function iterate_interface_fluxes(flux_formulation::SimilarityTheoryFluxes, approximate_interface_state) U = sqrt(Δu^2 + Δv^2 + Uᴳ^2) - + U = max(U, flux_formulation.minimum_velocity_scale) + # Compute roughness length scales ℓu₀ = roughness_length(ℓu, u★, U, 𝒬ₛ, ℂₐ) ℓq₀ = roughness_length(ℓq, ℓu₀, u★, U, 𝒬ₛ, ℂₐ) diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl index e7a3ca18c..15764b14c 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl @@ -5,6 +5,7 @@ using Oceananigans.Fields: Center using Oceananigans.Grids: grid_name using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series using Oceananigans.TimeSteppers: Clock, tick! +using Oceananigans.Simulations: TimeStepWizard using Oceananigans.Utils: prettysummary, Time using Adapt @@ -371,6 +372,8 @@ end return nothing end +(wizard::TimeStepWizard)(atmos::PrescribedAtmosphere) = Inf + @inline thermodynamics_parameters(atmos::Nothing) = nothing @inline thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters @inline surface_layer_height(atmos::PrescribedAtmosphere) = atmos.surface_layer_height diff --git a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl index 0c8d3cc47..937d11654 100644 --- a/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl @@ -7,6 +7,7 @@ using .InterfaceComputations: using ClimaSeaIce: SeaIceModel, SeaIceThermodynamics using Oceananigans.Grids: φnode +using Oceananigans.Simulations: TimeStepWizard using Printf @@ -20,10 +21,10 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_ # TODO after ice time-step: # - Adjust ocean heat flux if the ice completely melts? - time_step!(ocean, Δt) + !isnothing(ocean) && time_step!(ocean, Δt) # Time step the atmosphere - time_step!(atmosphere, Δt) + !isnothing(atmosphere) && time_step!(atmosphere, Δt) # TODO: # - Store fractional ice-free / ice-covered _time_ for more @@ -51,3 +52,27 @@ function update_state!(coupled_model::OceanSeaIceModel, callbacks=[]; compute_te return nothing end + +function (wizard::TimeStepWizard)(simulation::Simulation{<:OceanSeaIceModel}) + model = simulation.model + ocean_Δt = wizard(model.ocean) + + sea_ice_Δt = if isnothing(model.sea_ice) + Inf + else + wizard(model.sea_ice) + end + + atmosphere_Δt = if isnothing(model.atmosphere) + Inf + else + wizard(model.atmosphere) + end + + Δt = min(ocean_Δt, sea_ice_Δt, atmosphere_Δt) + + simulation.Δt = Δt + + return nothing +end + diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index c09803132..0434dc0ec 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -14,10 +14,11 @@ using Oceananigans.Operators using ClimaSeaIce using ClimaSeaIce: SeaIceModel, SlabSeaIceThermodynamics, PhaseTransitions, ConductiveFlux using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using ClimaSeaIce.Rheologies using ClimaSeaIce.SeaIceDynamics: SplitExplicitSolver, SemiImplicitStress, SeaIceMomentumEquation, StressBalanceFreeDrift using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology -using ClimaOcean.OceanSimulations: Default +using ClimaOcean.OceanSimulations: Default, u_immersed_bottom_drag, v_immersed_bottom_drag g_Earth = Oceananigans.defaults.gravitational_acceleration Ω_Earth = Oceananigans.defaults.planet_rotation_rate @@ -62,6 +63,14 @@ function sea_ice_simulation(grid, ocean=nothing; bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) + u_bc = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true, parameters=1e-1) + v_bc = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true, parameters=1e-1) + + # immersed_u_bc = ImmersedBoundaryConditions(top=nothing, bottom=nothing, west=nothing, east=nothing, south=u_bc, north=u_bc) + # immersed_v_bc = ImmersedBoundaryConditions(top=nothing, bottom=nothing, south=nothing, north=nothing, west=v_bc, east=v_bc) + # u_bcs = FieldBoundaryConditions(grid, (Face(), Center(), nothing); immersed = immersed_u_bc) + # v_bcs = FieldBoundaryConditions(grid, (Center(), Face(), nothing); immersed = immersed_v_bc) + # Build the sea ice model sea_ice_model = SeaIceModel(grid; ice_salinity, @@ -81,12 +90,17 @@ function sea_ice_simulation(grid, ocean=nothing; return sea_ice end +default_solver(::Nothing) = SplitExplicitSolver(120) +default_solver(ocean::Simulation) = default_solver(ocean.model.timestepper) +default_solver(::Oceananigans.TimeSteppers.QuasiAdamsBashforth2TimeStepper) = SplitExplicitSolver(120) +default_solver(::Oceananigans.TimeSteppers.SplitRungeKutta3TimeStepper) = SplitExplicitSolver(360) + function sea_ice_dynamics(grid, ocean=nothing; - sea_ice_ocean_drag_coefficient = 5.5e-3, + sea_ice_ocean_drag_coefficient = 2.5e-3, rheology = ElastoViscoPlasticRheology(), coriolis = nothing, free_drift = nothing, - solver = SplitExplicitSolver(120)) + solver = default_solver(ocean)) if isnothing(ocean) SSU = Oceananigans.Fields.ZeroField() diff --git a/test/test_jra55.jl b/test/test_jra55.jl index ec250e7ba..0007d0660 100644 --- a/test/test_jra55.jl +++ b/test/test_jra55.jl @@ -1,5 +1,6 @@ include("runtests_setup.jl") +using ClimaOcean.JRA55 using ClimaOcean.JRA55: download_JRA55_cache using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere