From 8ea4b6f216ddd6c7e152f46cfb576edfb2fae664 Mon Sep 17 00:00:00 2001 From: LenkaNovak Date: Wed, 26 Jul 2023 08:53:11 -0700 Subject: [PATCH] Add flux partitioning src + unit tests implemented in AMIP driver, conserving deps amip runs + custom sfc state new pipeline run clean up clean up revs rev2 rev3 clean format codecov code cov --- .buildkite/pipeline.yml | 7 + Project.toml | 4 + docs/src/fluxcalculator.md | 14 +- experiments/AMIP/modular/Manifest.toml | 2 +- experiments/AMIP/modular/cli_options.jl | 4 + .../components/atmosphere/climaatmos_init.jl | 137 +++++--- .../modular/components/land/bucket_init.jl | 3 +- .../modular/components/land/bucket_utils.jl | 34 +- .../components/ocean/slab_ocean_init.jl | 12 +- .../components/ocean/slab_seaice_init.jl | 12 +- .../AMIP/modular/coupler_driver_modular.jl | 68 +++- .../coupler_driver.jl | 2 +- .../ClimaCore/sea_breeze/Manifest.toml | 302 ++++++++--------- perf/Manifest.toml | 2 +- src/FieldExchanger.jl | 56 ++-- src/FluxCalculator.jl | 314 +++++++++++++++++- src/Interfacer.jl | 1 - test/Project.toml | 4 + test/field_exchanger_tests.jl | 14 +- test/flux_calculator_tests.jl | 299 ++++++++++++++++- 20 files changed, 1001 insertions(+), 290 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index ecdc0285c..0b07fb55c 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -139,6 +139,13 @@ steps: agents: slurm_mem: 20GB + - label: "Slabplanet: partitioned turbulent fluxes" + key: "slabplanet_partitioned_fluxes" + command: "julia --color=yes --project=experiments/AMIP/modular/ experiments/AMIP/modular/coupler_driver_modular.jl --run_name slabplanet_partitioned_fluxes --FLOAT_TYPE Float64 --coupled true --turb_flux_partition PartitionedStateFluxes --surface_setup PrescribedSurface --moist equil --vert_diff true --rad gray --energy_check true --mode_name slabplanet --t_end 10days --dt_save_to_sol 9days --dt_cpl 200 --dt 200secs --mono_surface true --h_elem 4 --precip_model 0M --anim true --job_id slabplanet_partitioned_fluxes" + artifact_paths: "experiments/AMIP/modular/output/slabplanet/slabplanet_partitioned_fluxes_artifacts/*" + agents: + slurm_mem: 20GB + # Test: non-monotonous remapping for land mask - label: "Slabplanet: non-monotonous surface remap" key: "slabplanet_non-monotonous" diff --git a/Project.toml b/Project.toml index 7e5c2ddf1..d57b61f24 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["CliMA Contributors "] version = "0.1.0" [deps] +CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" ClimaAtmos = "b2c96348-7fb7-4fe0-8da9-78d88439e717" ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d" ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" @@ -18,6 +19,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f" TempestRemap_jll = "8573a8c5-1df0-515e-a024-abad257ee284" @@ -26,6 +28,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] +CLIMAParameters = "0.7" ClimaAtmos = "0.15" ClimaComms = "0.5" ClimaCore = "0.10" @@ -39,6 +42,7 @@ OrdinaryDiffEq = "5, 6" Plots = "1" PrettyTables = "1, 2" SciMLBase = "1" +StaticArrays = "1" SurfaceFluxes = "0.6" TempestRemap_jll = "2" TerminalLoggers = "0.1" diff --git a/docs/src/fluxcalculator.md b/docs/src/fluxcalculator.md index 6dea0c856..d2f9016be 100644 --- a/docs/src/fluxcalculator.md +++ b/docs/src/fluxcalculator.md @@ -3,18 +3,18 @@ This modules contains abstract types and functions to calculate surface fluxes in the coupler, or to call flux calculating functions from the component models. Fluxes over a heterogeneous surface (e.g., from a gridpoint where atmospheric cell is overlying both land and ocean) can be handled in two different ways: -1. **Combined fluxes** (called with `CombinedAtmosGrid()`) +1. **Combined fluxes** (called with `CombinedStateFluxes()`) - these are calculated by averaging the surface properties for each gridpoint (e.g., land and ocean temperatures, albedos and roughness lengths are averaged, based on their respective area fractions), so the flux is calculated only once per gridpoint of the grid where we calculate fluxes. This is computationally faster, but it makes the fluxes on surface boundaries more diffuse. Currently, we use this method for calculating radiative fluxes in the atmosphere, and turbulent fluxes in the coupler (on the atmospheric grid). The fluxes are calculated in the atmospheric (host) model's cache, which can be setup to avoid allocating coupler fields. -2. **Partitioned fluxes** (called with `PartitionedComponentModelGrid()`) - - these are calculated separately for each surface type. It is then the fluxes (rather than the surface states) that are combined and passed to the atmospheric model as a boundary condition. This method ensures that each surface model receives fluxes that correspond to its state properties, resulting in a more accurate model evolution. (TODO: implementation in the next PR) +2. **Partitioned fluxes** (called with `PartitionedStateFluxes()`) + - these are calculated separately for each surface type. It is then the fluxes (rather than the surface states) that are combined and passed to the atmospheric model as a boundary condition. This method ensures that each surface model receives fluxes that correspond to its state properties, resulting in a more accurate model evolution. However, it is more computationally expensive, and requires more communication between the component models. ## FluxCalculator API ```@docs ClimaCoupler.FluxCalculator.TurbulentFluxPartition - ClimaCoupler.FluxCalculator.PartitionedComponentModelGrid - ClimaCoupler.FluxCalculator.CombinedAtmosGrid - ClimaCoupler.FluxCalculator.compute_combined_turbulent_fluxes! - ClimaCoupler.FluxCalculator.compute_atmos_turbulent_fluxes! + ClimaCoupler.FluxCalculator.PartitionedStateFluxes + ClimaCoupler.FluxCalculator.CombinedStateFluxes + ClimaCoupler.FluxCalculator.combined_turbulent_fluxes! + ClimaCoupler.FluxCalculator.atmos_turbulent_fluxes! ``` diff --git a/experiments/AMIP/modular/Manifest.toml b/experiments/AMIP/modular/Manifest.toml index 96ea15733..ada6a93d4 100644 --- a/experiments/AMIP/modular/Manifest.toml +++ b/experiments/AMIP/modular/Manifest.toml @@ -239,7 +239,7 @@ uuid = "d934ef94-cdd4-4710-83d6-720549644b70" version = "0.3.5" [[deps.ClimaCoupler]] -deps = ["ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] +deps = ["CLIMAParameters", "ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] path = "../../.." uuid = "4ade58fe-a8da-486c-bd89-46df092ec0c7" version = "0.1.0" diff --git a/experiments/AMIP/modular/cli_options.jl b/experiments/AMIP/modular/cli_options.jl index 0538d6731..4aa925e08 100644 --- a/experiments/AMIP/modular/cli_options.jl +++ b/experiments/AMIP/modular/cli_options.jl @@ -27,6 +27,10 @@ function argparse_settings() help = "Boolean flag indicating whether (1st order) monotone and conservative remapping is applied." arg_type = Bool default = false + "--turb_flux_partition" + help = "Method to partition turbulent fluxes. [`PartitionedStateFluxes`, `CombinedStateFluxes`]" + arg_type = String + default = "CombinedStateFluxes" "--albedo_from_file" help = "Access land surface albedo information from data file" arg_type = Bool diff --git a/experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl b/experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl index 97c82c9ca..423669261 100644 --- a/experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl +++ b/experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl @@ -1,10 +1,17 @@ # atmos_init: for ClimaAtmos pre-AMIP interface import ClimaAtmos using ClimaAtmos: RRTMGPI -import ClimaCoupler.FluxCalculator: compute_atmos_turbulent_fluxes!, calculate_surface_air_density +import ClimaAtmos: CT1, CT2, CT12, CT3, C3, C12, unit_basis_vector_data, ⊗ +import ClimaCoupler.FluxCalculator: + atmos_turbulent_fluxes!, + calculate_surface_air_density, + PartitionedStateFluxes, + extrapolate_ρ_to_sfc, + get_surface_params using ClimaCore: Fields.level, Geometry import ClimaCoupler.FieldExchanger: get_thermo_params import ClimaCoupler.Interfacer: get_field, update_field!, name, get_model_state_vector +using StaticArrays # the clima atmos `integrator` is now defined struct ClimaAtmosSimulation{P, Y, D, I} <: AtmosModelSimulation @@ -39,7 +46,7 @@ function atmos_init(::Type{FT}, parsed_args::Dict) where {FT} ClimaAtmosSimulation(integrator.p.params, Y, spaces, integrator) end -# required by the Interfacer +# extensions required by the Interfacer get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux}) = Fields.level(sim.integrator.p.ᶠradiation_flux, half) get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) = @@ -53,10 +60,28 @@ get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_moisture_flux}) = get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) = Spaces.level(sim.integrator.p.ᶜts, 1) +# extensions required by FluxCalculator (partitioned fluxes) +get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) = + Spaces.level(Fields.coordinate_field(sim.integrator.u.c).z, 1) +get_field(sim::ClimaAtmosSimulation, ::Val{:height_sfc}) = + Spaces.level(Fields.coordinate_field(sim.integrator.u.f).z, half) +function get_field(sim::ClimaAtmosSimulation, ::Val{:uv_int}) + uₕ_int = Geometry.UVVector.(Spaces.level(sim.integrator.u.c.uₕ, 1)) + return @. StaticArrays.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2) +end +get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) = TD.air_density.(thermo_params, sim.integrator.p.ᶜts) +get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) = TD.air_temperature.(thermo_params, sim.integrator.p.ᶜts) +get_field(sim::ClimaAtmosSimulation, ::Val{:cv_m}) = TD.cv_m.(thermo_params, sim.integrator.p.ᶜts) +get_field(sim::ClimaAtmosSimulation, ::Val{:gas_constant_air}) = + TD.gas_constant_air.(thermo_params, sim.integrator.p.ᶜts) + +get_surface_params(sim::ClimaAtmosSimulation) = CAP.surface_fluxes_params(sim.integrator.p.params) -function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_temperature}, field) - sim.integrator.p.radiation_model.surface_temperature .= RRTMGPI.field2array(field) +# extensions required by the Interfacer +function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_temperature}, csf) + sim.integrator.p.radiation_model.surface_temperature .= RRTMGPI.field2array(csf.T_S) end + function update_field!(sim::ClimaAtmosSimulation, ::Val{:albedo}, field) sim.integrator.p.radiation_model.diffuse_sw_surface_albedo .= reshape(RRTMGPI.field2array(field), 1, length(parent(field))) @@ -64,22 +89,45 @@ function update_field!(sim::ClimaAtmosSimulation, ::Val{:albedo}, field) reshape(RRTMGPI.field2array(field), 1, length(parent(field))) end -step!(sim::ClimaAtmosSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true) +# get_surface_params required by FluxCalculator (partitioned fluxes) +function update_field!(sim::ClimaAtmosSimulation, ::Val{:turbulent_fluxes}, fields) + (; F_turb_energy, F_turb_moisture, F_turb_ρτxz, F_turb_ρτyz) = fields + + Y = sim.integrator.u + surface_local_geometry = Fields.level(Fields.local_geometry_field(Y.f), Fields.half) + surface_normal = @. C3(unit_basis_vector_data(C3, surface_local_geometry)) + + # get template objects for the contravariant components of the momentum fluxes (required by Atmos boundary conditions) + vec_ct12_ct1 = @. CT12(CT2(unit_basis_vector_data(CT1, surface_local_geometry)), surface_local_geometry) + vec_ct12_ct2 = @. CT12(CT2(unit_basis_vector_data(CT2, surface_local_geometry)), surface_local_geometry) + + sim.integrator.p.sfc_conditions.ρ_flux_uₕ .= ( + surface_normal .⊗ + C12.( + swap_space!(ones(axes(vec_ct12_ct1)), F_turb_ρτxz) .* vec_ct12_ct1 .+ + swap_space!(ones(axes(vec_ct12_ct2)), F_turb_ρτyz) .* vec_ct12_ct2, + surface_local_geometry, + ) + ) -reinit!(sim::ClimaAtmosSimulation) = reinit!(sim.integrator) + parent(sim.integrator.p.sfc_conditions.ρ_flux_h_tot) .= parent(F_turb_energy) .* parent(surface_normal) # (shf + lhf) + parent(sim.integrator.p.sfc_conditions.ρ_flux_q_tot) .= parent(F_turb_moisture) .* parent(surface_normal) # (evap) -""" - update_sim!(atmos_sim::ClimaAtmosSimulation, csf) + # TODO: see if Atmos can rever to a simpler solution +end -Updates the surface fields for temperature and albedo. +# extensions required by FieldExchanger +step!(sim::ClimaAtmosSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true) + +reinit!(sim::ClimaAtmosSimulation) = reinit!(sim.integrator) -# Arguments -- `atmos_sim`: [ClimaAtmosSimulation] containing an atmospheric model simulation object. -- `csf`: [NamedTuple] containing coupler fields. -""" function update_sim!(atmos_sim::ClimaAtmosSimulation, csf, turbulent_fluxes) update_field!(atmos_sim, Val(:albedo), csf.albedo) - update_field!(atmos_sim, Val(:surface_temperature), csf.T_S) + update_field!(atmos_sim, Val(:surface_temperature), csf) + + if turbulent_fluxes isa PartitionedStateFluxes + update_field!(atmos_sim, Val(:turbulent_fluxes), csf) + end end @@ -112,35 +160,48 @@ function coupler_surface_setup( end """ - compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) + get_new_cache(atmos_sim::ClimaAtmosSimulation, csf) -Computes turbulent surface fluxes using ClimaAtmos's `update_surface_conditions!`. This -requires that we define a new temporary parameter Tuple, `new_p`, and save the new surface state -in it. We do not want `new_p` to live in the atmospheric model permanently, because that would also -trigger flux calculation during Atmos `step!`. We only want to trigger this once per coupling -timestep from ClimaCoupler. +Returns a new `p` with the updated surface conditions. """ +function get_new_cache(atmos_sim::ClimaAtmosSimulation, csf) -function compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) + p = atmos_sim.integrator.p - function modified_atmos_p(atmos_sim, csf_sfc) + csf_sfc = (; T = csf.T_S, z0m = csf.z0m_S, z0b = csf.z0b_S, beta = csf.beta, q_vap = csf.q_sfc) + modified_atmos_cache(atmos_sim, csf_sfc) +end - p = atmos_sim.integrator.p +""" + modified_atmos_cache(atmos_sim, csf_sfc) - coupler_sfc_setup = coupler_surface_setup(CoupledMoninObukhov(), p; csf_sfc = csf_sfc) +Returns a new `p` with the updated surface conditions. +""" +function modified_atmos_cache(atmos_sim, csf_sfc) - p_names = propertynames(p) - p_values = map(x -> x == :sfc_setup ? coupler_sfc_setup : getproperty(p, x), p_names) # TODO: use merge here + p = atmos_sim.integrator.p - (; zip(p_names, p_values)...) - end + coupler_sfc_setup = coupler_surface_setup(CoupledMoninObukhov(), p; csf_sfc = csf_sfc) - p = atmos_sim.integrator.p + p_names = propertynames(p) + p_values = map(x -> x == :sfc_setup ? coupler_sfc_setup : getproperty(p, x), p_names) # TODO: use merge here - csf_sfc = (; T = csf.T_S, z0m = csf.z0m_S, z0b = csf.z0b_S, beta = csf.beta, q_vap = csf.q_sfc) - new_p = modified_atmos_p(atmos_sim, csf_sfc) + (; zip(p_names, p_values)...) +end + +""" + atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) + +Computes turbulent surface fluxes using ClimaAtmos's `update_surface_conditions!`. This +requires that we define a new temporary parameter Tuple, `new_p`, and save the new surface state +in it. We do not want `new_p` to live in the atmospheric model permanently, because that would also +trigger flux calculation during Atmos `step!`. We only want to trigger this once per coupling +timestep from ClimaCoupler. +""" +function atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) + new_p = get_new_cache(atmos_sim, csf) ClimaAtmos.SurfaceConditions.update_surface_conditions!(atmos_sim.integrator.u, new_p, atmos_sim.integrator.t) - p.sfc_conditions .= new_p.sfc_conditions + atmos_sim.integrator.p.sfc_conditions .= new_p.sfc_conditions end """ @@ -161,18 +222,6 @@ function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::Fie extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), T_S)) end -""" - extrapolate_ρ_to_sfc(thermo_params, ts_int, T_sfc) - -Uses the ideal gas law and hydrostatic balance to extrapolate for surface density pointwise. -""" -function extrapolate_ρ_to_sfc(thermo_params, ts_in, T_sfc) - T_int = TD.air_temperature(thermo_params, ts_in) - Rm_int = TD.gas_constant_air(thermo_params, ts_in) - ρ_air = TD.air_density(thermo_params, ts_in) - ρ_air * (T_sfc / T_int)^(TD.cv_m(thermo_params, ts_in) / Rm_int) -end - """ get_model_state_vector(sim::ClimaAtmosSimulation) diff --git a/experiments/AMIP/modular/components/land/bucket_init.jl b/experiments/AMIP/modular/components/land/bucket_init.jl index d5fbbe2d4..6315a9efb 100644 --- a/experiments/AMIP/modular/components/land/bucket_init.jl +++ b/experiments/AMIP/modular/components/land/bucket_init.jl @@ -21,6 +21,7 @@ using ClimaLSM: import ClimaCoupler.Interfacer: LandModelSimulation, get_field, update_field!, name import ClimaCoupler.FieldExchanger: step!, reinit! +import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, surface_thermo_state """ BucketSimulation{M, Y, D, I} @@ -85,8 +86,6 @@ function ClimaLSM.Bucket.surface_fluxes( ) end - - """ net_radiation(radiation::CoupledRadiativeFluxes{FT}, model::BucketModel{FT}, diff --git a/experiments/AMIP/modular/components/land/bucket_utils.jl b/experiments/AMIP/modular/components/land/bucket_utils.jl index d5af95a56..0dfdce677 100644 --- a/experiments/AMIP/modular/components/land/bucket_utils.jl +++ b/experiments/AMIP/modular/components/land/bucket_utils.jl @@ -80,7 +80,7 @@ function make_lsm_domain( ) end -# required by Interfacer +# extensions required by Interfacer get_field(sim::BucketSimulation, ::Val{:surface_temperature}) = ClimaLSM.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t) get_field(sim::BucketSimulation, ::Val{:surface_humidity}) = @@ -92,11 +92,7 @@ get_field(sim::BucketSimulation, ::Val{:beta}) = get_field(sim::BucketSimulation, ::Val{:albedo}) = ClimaLSM.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p) get_field(sim::BucketSimulation, ::Val{:area_fraction}) = sim.area_fraction - - -# The surface air density is computed using the atmospheric state at the first level and making ideal gas -# and hydrostatic balance assumptions. The land model does not compute the surface air density so this is -# a reasonable stand-in. +get_field(sim::BucketSimulation, ::Val{:air_density}) = sim.integrator.p.bucket.ρ_sfc function update_field!(sim::BucketSimulation, ::Val{:turbulent_energy_flux}, field) parent(sim.integrator.p.bucket.turbulent_energy_flux) .= parent(field) @@ -123,6 +119,32 @@ step!(sim::BucketSimulation, t) = step!(sim.integrator, t - sim.integrator.t, tr reinit!(sim::BucketSimulation) = reinit!(sim.integrator) +# extensions required by FluxCalculator (partitioned fluxes) +function update_turbulent_fluxes_point!(sim::BucketSimulation, fields::NamedTuple, colidx::Fields.ColumnIndex) + (; F_turb_energy, F_turb_moisture) = fields + sim.integrator.p.bucket.turbulent_energy_flux[colidx] .= F_turb_energy + sim.integrator.p.bucket.evaporation[colidx] .= + F_turb_moisture ./ LSMP.ρ_cloud_liq(sim.model.parameters.earth_param_set) + return nothing +end + +# extension of FluxCalculator.surface_thermo_state, overriding the saturated-surface default +function surface_thermo_state( + sim::BucketSimulation, + thermo_params::TD.Parameters.ThermodynamicsParameters, + thermo_state_int, + colidx::ClimaCore.Fields.ColumnIndex, +) + + T_sfc = get_field(sim, Val(:surface_temperature), colidx) + # Note that the surface air density, ρ_sfc, is computed using the atmospheric state at the first level and making ideal gas + # and hydrostatic balance assumptions. The land model does not compute the surface air density so this is + # a reasonable stand-in. + ρ_sfc = get_field(sim, Val(:air_density), colidx) + q_sfc = get_field(sim, Val(:surface_humidity), colidx) # already calculated in rhs! (cache) + @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc) +end + """ get_model_state_vector(sim::BucketSimulation) diff --git a/experiments/AMIP/modular/components/ocean/slab_ocean_init.jl b/experiments/AMIP/modular/components/ocean/slab_ocean_init.jl index 0539847a9..4375de861 100644 --- a/experiments/AMIP/modular/components/ocean/slab_ocean_init.jl +++ b/experiments/AMIP/modular/components/ocean/slab_ocean_init.jl @@ -2,6 +2,7 @@ using ClimaCore import ClimaCoupler.Interfacer: OceanModelSimulation, get_field, update_field!, name import ClimaCoupler.FieldExchanger: step!, reinit! +import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point! """ SlabOceanSimulation{P, Y, D, I} @@ -17,6 +18,7 @@ struct SlabOceanSimulation{P, Y, D, I} <: OceanModelSimulation domain::D integrator::I end +name(::SlabOceanSimulation) = "SlabOceanSimulation" # ocean parameters struct OceanSlabParameters{FT <: AbstractFloat} @@ -103,7 +105,7 @@ Ensures that the space of the SST struct matches that of the land_fraction, and """ clean_sst(SST, _info) = (swap_space!(zeros(axes(_info.land_fraction)), SST) .+ float_type_bcf(_info)(273.15)) -# required by Interfacer +# extensions required by Interfacer get_field(sim::SlabOceanSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc get_field(sim::SlabOceanSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m @@ -111,6 +113,7 @@ get_field(sim::SlabOceanSimulation, ::Val{:roughness_buoyancy}) = sim.integrator get_field(sim::SlabOceanSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0) get_field(sim::SlabOceanSimulation, ::Val{:albedo}) = sim.integrator.p.params.α get_field(sim::SlabOceanSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction +get_field(sim::SlabOceanSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc function update_field!(sim::SlabOceanSimulation, ::Val{:area_fraction}, field::Fields.Field) sim.integrator.p.area_fraction .= field @@ -126,10 +129,17 @@ function update_field!(sim::SlabOceanSimulation, ::Val{:air_density}, field) parent(sim.integrator.p.ρ_sfc) .= parent(field) end +# extensions required by FieldExchanger step!(sim::SlabOceanSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true) reinit!(sim::SlabOceanSimulation) = reinit!(sim.integrator) +# extensions required by FluxCalculator (partitioned fluxes) +function update_turbulent_fluxes_point!(sim::SlabOceanSimulation, fields::NamedTuple, colidx::Fields.ColumnIndex) + (; F_turb_energy) = fields + @. sim.integrator.p.F_turb_energy[colidx] = F_turb_energy +end + """ get_model_state_vector(sim::SlabOceanSimulation) diff --git a/experiments/AMIP/modular/components/ocean/slab_seaice_init.jl b/experiments/AMIP/modular/components/ocean/slab_seaice_init.jl index 03aaae960..31b20a769 100644 --- a/experiments/AMIP/modular/components/ocean/slab_seaice_init.jl +++ b/experiments/AMIP/modular/components/ocean/slab_seaice_init.jl @@ -1,5 +1,6 @@ import ClimaCoupler.Interfacer: SeaIceModelSimulation, get_field, update_field!, name import ClimaCoupler.FieldExchanger: step!, reinit! +import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point! """ PrescribedIceSimulation{P, Y, D, I} @@ -26,6 +27,7 @@ struct PrescribedIceSimulation{P, Y, D, I} <: SeaIceModelSimulation domain::D integrator::I end +name(::PrescribedIceSimulation) = "PrescribedIceSimulation" # sea-ice parameters struct IceSlabParameters{FT <: AbstractFloat} @@ -116,7 +118,7 @@ clean_sic(SIC, _info) = swap_space!(zeros(axes(_info.land_fraction)), SIC) ./ fl get_ice_fraction(h_ice::FT, mono::Bool, threshold = 0.5) where {FT} = mono ? h_ice : Regridder.binary_mask(h_ice, threshold = FT(threshold)) -# required by Interfacer +# extensions required by Interfacer get_field(sim::PrescribedIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m @@ -124,6 +126,7 @@ get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) = sim.integr get_field(sim::PrescribedIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0) get_field(sim::PrescribedIceSimulation, ::Val{:albedo}) = sim.integrator.p.params.α get_field(sim::PrescribedIceSimulation, ::Val{:area_fraction}) = sim.integrator.p.area_fraction +get_field(sim::PrescribedIceSimulation, ::Val{:air_density}) = sim.integrator.p.ρ_sfc function update_field!(sim::PrescribedIceSimulation, ::Val{:area_fraction}, field::Fields.Field) sim.integrator.p.area_fraction .= field @@ -139,9 +142,16 @@ function update_field!(sim::PrescribedIceSimulation, ::Val{:air_density}, field) parent(sim.integrator.p.ρ_sfc) .= parent(field) end +# extensions required by FieldExchanger step!(sim::PrescribedIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true) reinit!(sim::PrescribedIceSimulation) = reinit!(sim.integrator) +# extensions required by FluxCalculator (partitioned fluxes) +function update_turbulent_fluxes_point!(sim::PrescribedIceSimulation, fields::NamedTuple, colidx::Fields.ColumnIndex) + (; F_turb_energy) = fields + @. sim.integrator.p.F_turb_energy[colidx] = F_turb_energy +end + """ get_model_state_vector(sim::PrescribedIceSimulation) diff --git a/experiments/AMIP/modular/coupler_driver_modular.jl b/experiments/AMIP/modular/coupler_driver_modular.jl index 1a1064744..7b5deb4be 100644 --- a/experiments/AMIP/modular/coupler_driver_modular.jl +++ b/experiments/AMIP/modular/coupler_driver_modular.jl @@ -83,6 +83,7 @@ if isinteractive() parsed_args["dt_cpl"] = 200 #hide parsed_args["dt"] = "200secs" #hide parsed_args["mono_surface"] = true #hide + parsed_args["turb_flux_partition"] = "PartitionedStateFluxes" #hide parsed_args["h_elem"] = 4 #hide # parsed_args["dt_save_restart"] = "5days" #hide parsed_args["precip_model"] = "0M" #hide @@ -131,7 +132,12 @@ import ClimaCoupler.Interfacer: get_field, update_field!, update_sim! -import ClimaCoupler.FluxCalculator: PartitionedComponentModelGrid, CombinedAtmosGrid, compute_combined_turbulent_fluxes! +import ClimaCoupler.FluxCalculator: + PartitionedStateFluxes, + CombinedStateFluxes, + combined_turbulent_fluxes!, + MoninObukhovScheme, + partitioned_turbulent_fluxes! import ClimaCoupler.FieldExchanger: import_atmos_fields!, import_combined_surface_fields!, @@ -182,7 +188,7 @@ We use a common `Space` for all global surfaces. This enables the MPI processes the atmospheric and surface components, so exchanges are parallelized. Note this is only possible when the atmosphere and surface are of the same horizontal resolution. =# -## init a 2D bounary space at the surface +## init a 2D boundary space at the surface boundary_space = atmos_sim.domain.face_space.horizontal_space # init land-sea fraction @@ -326,6 +332,8 @@ coupler_field_names = ( :beta, :F_turb_energy, :F_turb_moisture, + :F_turb_ρτxz, + :F_turb_ρτyz, :F_radiative, :P_liq, :P_snow, @@ -410,8 +418,14 @@ end #= ## Initialize Component Model Exchange =# - -turbulent_fluxes = CombinedAtmosGrid() +turbulent_fluxes = nothing +if parsed_args["turb_flux_partition"] == "PartitionedStateFluxes" + turbulent_fluxes = PartitionedStateFluxes() +elseif parsed_args["turb_flux_partition"] == "CombinedStateFluxes" + turbulent_fluxes = CombinedStateFluxes() +else + error("turb_flux_partition must be either PartitionedStateFluxes or CombinedStateFluxes") +end # 1) coupler combines surface states and calculates rho_sfc using surface and atmos variables update_surface_fractions!(cs) @@ -419,20 +433,35 @@ import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, tur import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) update_model_sims!(cs.model_sims, cs.fields, turbulent_fluxes) -# 2) each component model calculates its own vapor specific humidity (q_sfc) +# 2) each surface component model calculates its own vapor specific humidity (q_sfc) # TODO: the q_sfc calculation follows the design of the bucket q_sfc, but it would be neater to abstract this from step! step!(land_sim, Δt_cpl) step!(ocean_sim, Δt_cpl) step!(ice_sim, Δt_cpl) -# 3) coupler re-imports updated surface fields and calculates turbulent fluxes on the atmos grid, while updating atmos sfc_conditions -import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta -compute_combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) +# 3) coupler re-imports updated surface fields and calculates turbulent fluxes, while updating atmos sfc_conditions +if turbulent_fluxes isa CombinedStateFluxes + # calculate fluxes using combined surface states on the atmos grid + import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta, q_sfc + combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) # this updates the atmos thermo state, sfc_ts +elseif turbulent_fluxes isa PartitionedStateFluxes + # calculate turbulent fluxes in surface models and save the weighted average in coupler fields + partitioned_turbulent_fluxes!(cs.model_sims, cs.fields, cs.boundary_space, MoninObukhovScheme(), thermo_params) + + # update atmos sfc_conditions for surface temperature + # TODO: this is hard coded and needs to be simplified (need CA modification) + new_p = get_new_cache(atmos_sim, cs.fields) + ClimaAtmos.SurfaceConditions.update_surface_conditions!(atmos_sim.integrator.u, new_p, atmos_sim.integrator.t) # sets T_sfc (but SF calculation not necessary - CA) + atmos_sim.integrator.p.sfc_conditions .= new_p.sfc_conditions +end # 4) given the new sfc_conditions, atmos calls the radiative flux callback reinit!(atmos_sim) # sets a nonzero radiation flux -# 5) coupler re-imports updated atmos fluxes, and sends them to surface models +# 5) coupler re-imports updated atmos fluxes (radiative fluxes for both `turbulent_fluxes` types +# and also turbulent fluxes if `turbulent_fluxes isa CombinedStateFluxes`, +# and sends them to the surface component models. If `turbulent_fluxes isa PartitionedStateFluxes` +# atmos receives the turbulent fluxes from the coupler. import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) update_model_sims!(cs.model_sims, cs.fields, turbulent_fluxes) @@ -494,11 +523,24 @@ function solve_coupler!(cs) ## exchange combined fields and (if specified) calculate fluxes using combined states update_surface_fractions!(cs) - import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) - compute_combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) - import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # radiative and/or turbulent + import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta + if turbulent_fluxes isa CombinedStateFluxes + combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) # this updates the surface thermo state, sfc_ts, in ClimaAtmos (but also unnecessarily calculates fluxes) + elseif turbulent_fluxes isa PartitionedStateFluxes + # calculate turbulent fluxes in surfaces and save the weighted average in coupler fields + partitioned_turbulent_fluxes!(cs.model_sims, cs.fields, cs.boundary_space, MoninObukhovScheme(), thermo_params) + + # update atmos sfc_conditions for surface temperature - TODO: this needs to be simplified (need CA modification) + new_p = get_new_cache(atmos_sim, cs.fields) + ClimaAtmos.SurfaceConditions.update_surface_conditions!( + atmos_sim.integrator.u, + new_p, + atmos_sim.integrator.t, + ) # to set T_sfc (but SF calculation not necessary - CA modification) + atmos_sim.integrator.p.sfc_conditions .= new_p.sfc_conditions + end - # TODO: compute_and_send_partitioned_turbulent_fluxes!(cs) + import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # radiative and/or turbulent ## monthly callbacks if trigger_callback(cs, Monthly()) diff --git a/experiments/AMIP/moist_mpi_earth_dynamical_sea_ice/coupler_driver.jl b/experiments/AMIP/moist_mpi_earth_dynamical_sea_ice/coupler_driver.jl index 5fe8bc230..d6241cc24 100644 --- a/experiments/AMIP/moist_mpi_earth_dynamical_sea_ice/coupler_driver.jl +++ b/experiments/AMIP/moist_mpi_earth_dynamical_sea_ice/coupler_driver.jl @@ -35,7 +35,7 @@ include("mpi/mpi_init.jl") include("atmos/atmos_init.jl") atmos_sim = atmos_init(FT, Y, spaces, integrator, params = params); -# init a 2D bounary space at the surface, assuming the same instance (and MPI distribution if applicable) as the atmos domain above +# init a 2D boundary space at the surface, assuming the same instance (and MPI distribution if applicable) as the atmos domain above boundary_space = ClimaCore.Fields.level(atmos_sim.domain.face_space, half) # global surface grid # init land-sea mask diff --git a/experiments/ClimaCore/sea_breeze/Manifest.toml b/experiments/ClimaCore/sea_breeze/Manifest.toml index 1850f2d90..ec04e7619 100644 --- a/experiments/ClimaCore/sea_breeze/Manifest.toml +++ b/experiments/ClimaCore/sea_breeze/Manifest.toml @@ -5,9 +5,9 @@ manifest_format = "2.0" project_hash = "4aa4c869f9d982e7450ef4fd4b4c13c84420cf0d" [[deps.ADTypes]] -git-tree-sha1 = "e58c18d2312749847a74f5be80bb0fa53da102bd" +git-tree-sha1 = "f5c25e8a5b29b5e941b7408bc8cc79fea4d9ef9a" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "0.1.5" +version = "0.1.6" [[deps.AMD]] deps = ["Libdl", "LinearAlgebra", "SparseArrays", "Test"] @@ -17,9 +17,9 @@ version = "0.5.0" [[deps.AbstractFFTs]] deps = ["ChainRulesCore", "LinearAlgebra"] -git-tree-sha1 = "8bc0aaec0ca548eb6cf5f0d7d16351650c1ee956" +git-tree-sha1 = "cad4c758c0038eea30394b1b671526921ca85b21" uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" -version = "1.3.2" +version = "1.4.0" [[deps.AbstractTrees]] git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c" @@ -62,9 +62,9 @@ version = "0.1.29" [[deps.ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "9e5ad0f651c1be8a355e4d327e9cb3bc135d1b73" +git-tree-sha1 = "0adeb373cb2268a0f4e62d349ce9957b509d4d7e" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "1.0.8" +version = "1.1.0" [[deps.ArtifactWrappers]] deps = ["Downloads", "Pkg"] @@ -115,9 +115,9 @@ version = "0.1.5" [[deps.BlockArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"] -git-tree-sha1 = "ffb0191c5d1872a8c91a6181e94a365b7d08d5f2" +git-tree-sha1 = "174b4970af15a500a29e76151f5c53195784b9d4" uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -version = "0.16.35" +version = "0.16.36" [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -138,9 +138,9 @@ version = "0.1.2" [[deps.CLIMAParameters]] deps = ["DocStringExtensions", "TOML", "Test"] -git-tree-sha1 = "119ef280867dde6b8a3f351b36b1f5748db12faf" +git-tree-sha1 = "d44564edc0950638a868e26472d474f717575b3f" uuid = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" -version = "0.7.5" +version = "0.7.7" [[deps.CPUSummary]] deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] @@ -197,37 +197,37 @@ uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" version = "1.16.0" [[deps.ChangesOfVariables]] -deps = ["LinearAlgebra", "Test"] -git-tree-sha1 = "f84967c4497e0e1955f9a582c232b02847c5f589" +deps = ["InverseFunctions", "LinearAlgebra", "Test"] +git-tree-sha1 = "2fba81a302a7be671aefe194f0525ef231104e7f" uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" -version = "0.1.7" +version = "0.1.8" [[deps.ClimaAtmos]] deps = ["ArgParse", "ArtifactWrappers", "Artifacts", "AtmosphericProfilesLibrary", "CLIMAParameters", "CUDA", "ClimaComms", "ClimaCore", "ClimaTimeSteppers", "CloudMicrophysics", "Colors", "Dates", "Dierckx", "DiffEqBase", "DiffEqCallbacks", "Distributions", "DocStringExtensions", "FastGaussQuadrature", "ImageFiltering", "Insolation", "Interpolations", "IntervalSets", "JLD2", "LambertW", "LinearAlgebra", "Logging", "NCDatasets", "NVTX", "OrdinaryDiffEq", "Pkg", "Printf", "RRTMGP", "Random", "RootSolvers", "StaticArrays", "Statistics", "StatsBase", "SurfaceFluxes", "TerminalLoggers", "Test", "Thermodynamics", "YAML"] -git-tree-sha1 = "c6033e6377d1af82709871be44a82db0a2a520ce" +git-tree-sha1 = "6cf1b2c438239cc3bfd3b8cffe59fa50442f1996" uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717" -version = "0.15.1" +version = "0.15.3" [[deps.ClimaComms]] deps = ["CUDA", "MPI"] -git-tree-sha1 = "7f3f28296a39cf99cb1372a6f1fc47b0baa1c66d" +git-tree-sha1 = "6cbe1b835f91033e6326a92ac1a1f6a8c9b2215f" uuid = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d" -version = "0.4.2" +version = "0.5.3" [[deps.ClimaCore]] deps = ["Adapt", "BlockArrays", "CUDA", "ClimaComms", "CubedSphere", "DataStructures", "DiffEqBase", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "LinearAlgebra", "PkgVersion", "RecursiveArrayTools", "Requires", "RootSolvers", "SparseArrays", "Static", "StaticArrays", "Statistics", "UnPack"] -git-tree-sha1 = "fe69ca07e8f9db7034d8cd6fa816332da3463576" +git-tree-sha1 = "6d3ea05148235667c3e51bd0d8b03cd6416ce919" uuid = "d414da3d-4745-48bb-8d80-42e94e092884" -version = "0.10.41" +version = "0.10.44" [[deps.ClimaCoreTempestRemap]] -deps = ["ClimaComms", "ClimaCore", "Dates", "LinearAlgebra", "NCDatasets", "PkgVersion", "TempestRemap_jll", "Test"] -git-tree-sha1 = "40eec8cff925a004281f47dd5b32809db9e8d5e0" +deps = ["ClimaCore", "Dates", "LinearAlgebra", "NCDatasets", "PkgVersion", "TempestRemap_jll", "Test"] +git-tree-sha1 = "12713936cdc28a3442f61d92d782d5bd329ded7e" uuid = "d934ef94-cdd4-4710-83d6-720549644b70" -version = "0.3.9" +version = "0.3.5" [[deps.ClimaCoupler]] -deps = ["ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] +deps = ["CLIMAParameters", "ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] path = "../../.." uuid = "4ade58fe-a8da-486c-bd89-46df092ec0c7" version = "0.1.0" @@ -258,15 +258,15 @@ version = "0.10.3" [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83" +git-tree-sha1 = "02aa26a4cf76381be7f66e020a3eddeb27b0a092" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.1" +version = "0.7.2" [[deps.ColorSchemes]] deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] -git-tree-sha1 = "be6ab11021cd29f0344d5c4357b163af05a48cba" +git-tree-sha1 = "dd3000d954d483c1aad05fe1eb9e6a715c97013e" uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.21.0" +version = "3.22.0" [[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] @@ -275,10 +275,10 @@ uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" version = "0.11.4" [[deps.ColorVectorSpace]] -deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] -git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589" +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "Requires", "Statistics", "TensorCore"] +git-tree-sha1 = "a1f44953f2382ebb937d60dafbe2deea4bd23249" uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" -version = "0.9.10" +version = "0.10.0" [[deps.Colors]] deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] @@ -286,12 +286,6 @@ git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" version = "0.12.10" -[[deps.CommonDataModel]] -deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf"] -git-tree-sha1 = "2678b3fc170d582655a14d22867b031b6e43c2d4" -uuid = "1fbeeb36-5f17-413c-809b-666fb144f157" -version = "0.2.4" - [[deps.CommonSolve]] git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" @@ -305,9 +299,9 @@ version = "0.3.0" [[deps.Compat]] deps = ["Dates", "LinearAlgebra", "UUIDs"] -git-tree-sha1 = "4e88377ae7ebeaf29a047aa1ee40826e0b708a5d" +git-tree-sha1 = "5ce999a19f4ca23ea484e92a1774a61b8ca4cf8e" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.7.0" +version = "4.8.0" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -321,15 +315,15 @@ version = "0.3.2" [[deps.ConcurrentUtilities]] deps = ["Serialization", "Sockets"] -git-tree-sha1 = "96d823b94ba8d187a6d8f0826e731195a74b90e9" +git-tree-sha1 = "5372dbbf8f0bdb8c700db5367132925c0771ef7e" uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" -version = "2.2.0" +version = "2.2.1" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] -git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.2" +version = "1.5.3" [[deps.Contour]] git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" @@ -402,9 +396,9 @@ version = "0.1.0+0" [[deps.DiffEqBase]] deps = ["ArrayInterface", "ChainRulesCore", "DataStructures", "Distributions", "DocStringExtensions", "EnumX", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "Reexport", "Requires", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces", "ZygoteRules"] -git-tree-sha1 = "62c41421bd0facc43dfe4e9776135fe21fd1e1b9" +git-tree-sha1 = "c5692436e7f8279503466db216c74165d1b301e4" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.126.0" +version = "6.127.0" [[deps.DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "Markdown", "NLsolve", "Parameters", "RecipesBase", "RecursiveArrayTools", "SciMLBase", "StaticArraysCore"] @@ -426,9 +420,9 @@ version = "1.15.1" [[deps.Distances]] deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "49eba9ad9f7ead780bfb7ee319f962c811c6d3b2" +git-tree-sha1 = "b6def76ffad15143924a2199f72a5cd883a2e8a9" uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" -version = "0.10.8" +version = "0.10.9" [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] @@ -522,9 +516,9 @@ version = "3.3.10+0" [[deps.FastBroadcast]] deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"] -git-tree-sha1 = "d1248fceea0b26493fd33e8e9e8c553270da03bd" +git-tree-sha1 = "aa9925a229d45fe3018715238956766fa21804d1" uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" -version = "0.2.5" +version = "0.2.6" [[deps.FastClosures]] git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" @@ -539,9 +533,9 @@ version = "0.5.1" [[deps.FastLapackInterface]] deps = ["LinearAlgebra"] -git-tree-sha1 = "c1293a93193f0ae94be7cf338d33e162c39d8788" +git-tree-sha1 = "b12f05108e405dadcc2aff0008db7f831374e051" uuid = "29a986be-02c6-4525-aec4-84b980013641" -version = "1.2.9" +version = "2.0.0" [[deps.FileIO]] deps = ["Pkg", "Requires", "UUIDs"] @@ -554,9 +548,9 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" [[deps.FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] -git-tree-sha1 = "0b3b52afd0f87b0a3f5ada0466352d125c9db458" +git-tree-sha1 = "f372472e8672b1d993e93dada09e23139b509f9e" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "1.2.1" +version = "1.5.0" [[deps.FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays", "StaticArrays"] @@ -589,10 +583,10 @@ uuid = "f6369f11-7733-5829-9624-2563aa707210" version = "0.10.35" [[deps.FreeType2_jll]] -deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9" +deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"] +git-tree-sha1 = "d8db6a5a2fe1381c1ea4ef2cab7c69c2de7f9ea0" uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" -version = "2.10.4+0" +version = "2.13.1+0" [[deps.FriBidi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -635,21 +629,21 @@ version = "0.1.5" [[deps.GPUCompiler]] deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] -git-tree-sha1 = "d60b5fe7333b5fa41a0378ead6614f1ab51cf6d0" +git-tree-sha1 = "72b2e3c2ba583d1a7aa35129e56cf92e07c083e3" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" -version = "0.21.3" +version = "0.21.4" [[deps.GR]] deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"] -git-tree-sha1 = "8b8a2fd4536ece6e554168c21860b6820a8a83db" +git-tree-sha1 = "d73afa4a2bb9de56077242d98cf763074ab9a970" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.72.7" +version = "0.72.9" [[deps.GR_jll]] -deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "19fad9cd9ae44847fe842558a744748084a722d1" +deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt6Base_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "f61f768bf090d97c532d24b64e07b237e9bb7b6b" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.72.7+0" +version = "0.72.9+0" [[deps.GaussQuadrature]] deps = ["SpecialFunctions"] @@ -680,12 +674,6 @@ git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23" uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" version = "2.74.0+2" -[[deps.Graphics]] -deps = ["Colors", "LinearAlgebra", "NaNMath"] -git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd" -uuid = "a2bd30eb-e257-5431-a919-1863eab51364" -version = "1.1.2" - [[deps.Graphite2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" @@ -717,9 +705,9 @@ version = "1.12.2+2" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "7f5ef966a02a8fdf3df2ca03108a88447cb3c6f0" +git-tree-sha1 = "cb56ccdd481c0dd7f975ad2b3b62d9eda088f7e2" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.9.8" +version = "1.9.14" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] @@ -735,9 +723,9 @@ version = "0.1.15" [[deps.HypergeometricFunctions]] deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] -git-tree-sha1 = "0ec02c648befc2f94156eaef13b0f38106212f3f" +git-tree-sha1 = "83e95aaab9dc184a6dcd9c4c52aa0dc26cd14a1d" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.17" +version = "0.3.21" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -746,21 +734,21 @@ version = "0.1.1" [[deps.ImageBase]] deps = ["ImageCore", "Reexport"] -git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835" +git-tree-sha1 = "eb49b82c172811fd2c86759fa0553a2221feb909" uuid = "c817782e-172a-44cc-b673-b171935fbb9e" -version = "0.1.5" +version = "0.1.7" [[deps.ImageCore]] -deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"] -git-tree-sha1 = "acf614720ef026d38400b3817614c45882d75500" +deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "PrecompileTools", "Reexport"] +git-tree-sha1 = "fc5d1d3443a124fde6e92d0260cd9e064eba69f8" uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" -version = "0.9.4" +version = "0.10.1" [[deps.ImageFiltering]] -deps = ["CatIndices", "ComputationalResources", "DataStructures", "FFTViews", "FFTW", "ImageBase", "ImageCore", "LinearAlgebra", "OffsetArrays", "Reexport", "SnoopPrecompile", "SparseArrays", "StaticArrays", "Statistics", "TiledIteration"] -git-tree-sha1 = "d90867cbe037730a73c9a9499b3591eedbe387a0" +deps = ["CatIndices", "ComputationalResources", "DataStructures", "FFTViews", "FFTW", "ImageBase", "ImageCore", "LinearAlgebra", "OffsetArrays", "PrecompileTools", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "TiledIteration"] +git-tree-sha1 = "c371a39622dc3b941ffd7c00e6b519d63b3f3f06" uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5" -version = "0.7.5" +version = "0.7.7" [[deps.Inflate]] git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428" @@ -797,9 +785,9 @@ version = "0.7.4" [[deps.InverseFunctions]] deps = ["Test"] -git-tree-sha1 = "6667aadd1cdee2c6cd068128b3d226ebc4fb0c67" +git-tree-sha1 = "eabe3125edba5c9c10b60a160b1779a000dc8b29" uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.9" +version = "0.1.11" [[deps.IrrationalConstants]] git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" @@ -813,9 +801,9 @@ version = "1.0.0" [[deps.JLD2]] deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "42c17b18ced77ff0be65957a591d34f4ed57c631" +git-tree-sha1 = "aa6ffef1fd85657f4999030c52eaeec22a279738" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.31" +version = "0.4.33" [[deps.JLFzf]] deps = ["Pipe", "REPL", "Random", "fzf_jll"] @@ -854,16 +842,16 @@ uuid = "ef3ab10e-7fda-4108-b977-705223b18434" version = "0.4.0" [[deps.KernelAbstractions]] -deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "b48617c5d764908b5fac493cd907cf33cc11eec1" +deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "4c5875e4c228247e1c2b087669846941fb6e0118" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.6" +version = "0.9.8" [[deps.Krylov]] deps = ["LinearAlgebra", "Printf", "SparseArrays"] -git-tree-sha1 = "0356a64062656b0cbb43c504ad5de338251f4bda" +git-tree-sha1 = "6dc4ad9cd74ad4ca0a8e219e945dbd22039f2125" uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" -version = "0.9.1" +version = "0.9.2" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -885,15 +873,15 @@ version = "3.0.0+1" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "7d5788011dd273788146d40eb5b1fbdc199d0296" +git-tree-sha1 = "8695a49bfe05a2dc0feeefd06b4ca6361a018729" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.0.1" +version = "6.1.0" [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "1222116d7313cdefecf3d45a2bc1a89c4e7c9217" +git-tree-sha1 = "c35203c1e1002747da220ffc3c0762ce7754b08c" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.22+0" +version = "0.0.23+0" [[deps.LLVMOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1004,10 +992,10 @@ uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" version = "2.35.0+0" [[deps.Libtiff_jll]] -deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"] -git-tree-sha1 = "3eb79b0ca5764d4799c06699573fd8f533259713" +deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "2da088d113af58221c52828a80378e16be7d037a" uuid = "89763e89-9b03-5906-acba-b20f662cd828" -version = "4.4.0+0" +version = "4.5.1+1" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1033,9 +1021,9 @@ version = "2.5.2" [[deps.LinearSolve]] deps = ["ArrayInterface", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "LinearAlgebra", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "SuiteSparse", "UnPack"] -git-tree-sha1 = "c6a6f78167d7b7c19dfb7148161d7f1962a0b361" +git-tree-sha1 = "1b55771f2c211583ad52af5a5ca6475be374c961" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "2.2.1" +version = "2.4.1" [[deps.LogExpFunctions]] deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] @@ -1054,9 +1042,9 @@ version = "1.0.0" [[deps.LoopVectorization]] deps = ["ArrayInterface", "ArrayInterfaceCore", "CPUSummary", "ChainRulesCore", "CloseOpenIntervals", "DocStringExtensions", "ForwardDiff", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "SpecialFunctions", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "e4eed22d70ac91d7a4bf9e0f6902383061d17105" +git-tree-sha1 = "c88a4afe1703d731b1c4fdf4e3c7e77e3b176ea2" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.162" +version = "0.12.165" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] @@ -1078,9 +1066,9 @@ version = "4.1.2+0" [[deps.MPIPreferences]] deps = ["Libdl", "Preferences"] -git-tree-sha1 = "d86a788b336e8ae96429c0c42740ccd60ac0dfcc" +git-tree-sha1 = "781916a2ebf2841467cda03b6f1af43e23839d85" uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" -version = "0.1.8" +version = "0.1.9" [[deps.MPItrampoline_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] @@ -1155,10 +1143,10 @@ uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" version = "0.2.4" [[deps.NCDatasets]] -deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "NetCDF_jll", "NetworkOptions", "Printf"] -git-tree-sha1 = "4263c4220f22e20729329838bf7e94a49d1ac32f" +deps = ["CFTime", "DataStructures", "Dates", "NetCDF_jll", "Printf"] +git-tree-sha1 = "17e39eb5bbe564f48bdbefbd103bd3f49fcfcb9b" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -version = "0.12.17" +version = "0.11.9" [[deps.NLSolversBase]] deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] @@ -1191,10 +1179,10 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "1.0.2" [[deps.NetCDF_jll]] -deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg", "XML2_jll", "Zlib_jll"] -git-tree-sha1 = "072f8371f74c3b9e1b26679de7fbf059d45ea221" +deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "XML2_jll", "Zlib_jll"] +git-tree-sha1 = "7f5a03e6712f5447c9c344430b8d1927a4777483" uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3" -version = "400.902.5+1" +version = "400.902.206+0" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" @@ -1208,9 +1196,9 @@ version = "1.8.0" [[deps.OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" +git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.12.9" +version = "1.12.10" [[deps.Ogg_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1265,15 +1253,15 @@ uuid = "91d4177d-7536-5919-b921-800302f37372" version = "1.3.2+0" [[deps.OrderedCollections]] -git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.0" +version = "1.6.2" [[deps.OrdinaryDiffEq]] deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "IfElse", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "NonlinearSolve", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLNLSolve", "SciMLOperators", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "SparseDiffTools", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "2aa75defe3eb7fcb0d914c0f7df907dbb8d63d3d" +git-tree-sha1 = "4f1ab68f236fa846d0c30718c6b29c1665b019c0" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -version = "6.53.2" +version = "6.53.3" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] @@ -1286,6 +1274,12 @@ git-tree-sha1 = "67eae2738d63117a196f497d7db789821bce61d1" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" version = "0.11.17" +[[deps.PackageExtensionCompat]] +deps = ["Requires", "TOML"] +git-tree-sha1 = "32f3d52212a8d1c5d589a58851b1f04c97339110" +uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930" +version = "1.0.0" + [[deps.PaddedViews]] deps = ["OffsetArrays"] git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f" @@ -1340,15 +1334,15 @@ version = "1.3.5" [[deps.Plots]] deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Preferences", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] -git-tree-sha1 = "75ca67b2c6512ad2d0c767a7cfc55e75075f8bbc" +git-tree-sha1 = "9f8675a55b37a70aa23177ec110f6e3f4dd68466" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.38.16" +version = "1.38.17" [[deps.Polyester]] deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"] -git-tree-sha1 = "0fe4e7c4d8ff4c70bfa507f0dd96fa161b115777" +git-tree-sha1 = "3d811babe092a6e7b130beee84998fe7663348b6" uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" -version = "0.7.3" +version = "0.7.5" [[deps.PolyesterWeave]] deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] @@ -1375,10 +1369,10 @@ uuid = "21216c6a-2e73-6563-6e65-726566657250" version = "1.4.0" [[deps.PrettyTables]] -deps = ["Crayons", "Formatting", "LaTeXStrings", "Markdown", "Reexport", "StringManipulation", "Tables"] -git-tree-sha1 = "213579618ec1f42dea7dd637a42785a608b1ea9c" +deps = ["Crayons", "LaTeXStrings", "Markdown", "Printf", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "ee094908d720185ddbdc58dbe0c1cbe35453ec7a" uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" -version = "2.2.4" +version = "2.2.7" [[deps.Printf]] deps = ["Unicode"] @@ -1396,11 +1390,11 @@ git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539" uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c" version = "0.1.4" -[[deps.Qt5Base_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] -git-tree-sha1 = "0c03844e2231e12fda4d0086fd7cbe4098ee8dc5" -uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1" -version = "5.15.3+2" +[[deps.Qt6Base_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] +git-tree-sha1 = "364898e8f13f7eaaceec55fd3d08680498c0aa6e" +uuid = "c0090381-4147-56d7-9ebc-da0b1113ec56" +version = "6.4.2+3" [[deps.QuadGK]] deps = ["DataStructures", "LinearAlgebra"] @@ -1413,10 +1407,10 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[deps.RRTMGP]] -deps = ["Adapt", "CUDA", "DocStringExtensions", "GaussQuadrature", "Random", "StaticArrays"] -git-tree-sha1 = "13e066c65ac3b0b4f7913856677fd4fcaeb24ccb" +deps = ["Adapt", "CUDA", "ClimaComms", "DocStringExtensions", "GaussQuadrature", "Random", "StaticArrays"] +git-tree-sha1 = "87ddffcc9520b7d721a03e2755412d8d6b674d5f" uuid = "a01a1ee8-cea4-48fc-987c-fc7878d79da1" -version = "0.8.3" +version = "0.9.1" [[deps.Random]] deps = ["SHA", "Serialization"] @@ -1454,9 +1448,9 @@ version = "0.6.12" [[deps.RecursiveArrayTools]] deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "02ef02926f30d53b94be443bfaea010c47f6b556" +git-tree-sha1 = "7ed35fb5f831aaf09c2d7c8736d44667a1afdcb0" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "2.38.5" +version = "2.38.7" [[deps.RecursiveFactorization]] deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "SnoopPrecompile", "StrideArraysCore", "TriangularSolve"] @@ -1522,9 +1516,9 @@ version = "0.6.39" [[deps.SciMLBase]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] -git-tree-sha1 = "a22a12db91f6a921e28a7ae39a9546eed93fd92e" +git-tree-sha1 = "04370090604cd399db5bebddb636d80ab9d338e9" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "1.93.0" +version = "1.94.0" [[deps.SciMLNLSolve]] deps = ["DiffEqBase", "LineSearches", "NLsolve", "Reexport", "SciMLBase"] @@ -1534,9 +1528,9 @@ version = "0.1.8" [[deps.SciMLOperators]] deps = ["ArrayInterface", "DocStringExtensions", "Lazy", "LinearAlgebra", "Setfield", "SparseArrays", "StaticArraysCore", "Tricks"] -git-tree-sha1 = "668ed6b2265d1adb36b7239501e9369025bd6c78" +git-tree-sha1 = "745755a5b932c9a664d7e9e4beb60c692b211d4b" uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -version = "0.3.2" +version = "0.3.5" [[deps.Scratch]] deps = ["Dates"] @@ -1569,10 +1563,10 @@ uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" version = "1.1.0" [[deps.SimpleNonlinearSolve]] -deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "PrecompileTools", "Reexport", "Requires", "SciMLBase", "StaticArraysCore"] -git-tree-sha1 = "56aa73a93cdca493af5155a0338a864ed314222b" +deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "PackageExtensionCompat", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"] +git-tree-sha1 = "91fcc402c4ab978ad5759489db9a9c5a71732f2d" uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7" -version = "0.1.16" +version = "0.1.18" [[deps.SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] @@ -1630,9 +1624,9 @@ version = "0.1.1" [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "dbde6766fc677423598138a5951269432b0fcc90" +git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.7" +version = "0.8.8" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"] @@ -1647,9 +1641,9 @@ uuid = "90137ffa-7385-5640-81b9-e52037218182" version = "1.5.26" [[deps.StaticArraysCore]] -git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.0" +version = "1.4.2" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -1675,15 +1669,15 @@ version = "1.3.0" [[deps.StrideArraysCore]] deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"] -git-tree-sha1 = "26316a9a43d32ecd0c5b342cb36dc86305c73e53" +git-tree-sha1 = "f02eb61eb5c97b48c153861c72fbbfdddc607e06" uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" -version = "0.4.16" +version = "0.4.17" [[deps.StringEncodings]] deps = ["Libiconv_jll"] -git-tree-sha1 = "33c0da881af3248dafefb939a21694b97cfece76" +git-tree-sha1 = "b765e46ba27ecf6b44faf70df40c57aa3a547dcb" uuid = "69024149-9ee7-55f6-a4c4-859efe599b68" -version = "0.3.6" +version = "0.3.7" [[deps.StringManipulation]] git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123" @@ -1809,9 +1803,9 @@ version = "0.1.7" [[deps.TruncatedStacktraces]] deps = ["InteractiveUtils", "MacroTools", "Preferences"] -git-tree-sha1 = "7bc1632a4eafbe9bd94cf1a784a9a4eb5e040a91" +git-tree-sha1 = "ea3e54c2bdde39062abf5a9758a23735558705e1" uuid = "781d530d-4396-4725-bb49-402e4bee1e77" -version = "1.3.0" +version = "1.4.0" [[deps.URIs]] git-tree-sha1 = "074f993b0ca030848b897beff716d93aca60f06a" @@ -1837,10 +1831,10 @@ uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" version = "0.4.1" [[deps.Unitful]] -deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"] -git-tree-sha1 = "ba4aa36b2d5c98d6ed1f149da916b3ba46527b2b" +deps = ["ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "Random"] +git-tree-sha1 = "c4d2a349259c8eba66a00a540d550f122a3ab228" uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.14.0" +version = "1.15.0" [[deps.UnitfulLatexify]] deps = ["LaTeXStrings", "Latexify", "Unitful"] @@ -1906,6 +1900,12 @@ git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" uuid = "aed1982a-8fda-507f-9586-7b0439959a61" version = "1.1.34+0" +[[deps.XZ_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "2222b751598bd9f4885c9ce9cd23e83404baa8ce" +uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800" +version = "5.4.3+1" + [[deps.Xorg_libX11_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495" diff --git a/perf/Manifest.toml b/perf/Manifest.toml index 0501befd7..606389fb2 100644 --- a/perf/Manifest.toml +++ b/perf/Manifest.toml @@ -239,7 +239,7 @@ uuid = "d934ef94-cdd4-4710-83d6-720549644b70" version = "0.3.5" [[deps.ClimaCoupler]] -deps = ["ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] +deps = ["CLIMAParameters", "ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "StaticArrays", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"] path = ".." uuid = "4ade58fe-a8da-486c-bd89-46df092ec0c7" version = "0.1.0" diff --git a/src/FieldExchanger.jl b/src/FieldExchanger.jl index d7c152856..1f8a858ac 100644 --- a/src/FieldExchanger.jl +++ b/src/FieldExchanger.jl @@ -28,18 +28,21 @@ have to be defined for the amtospheric component model type. function import_atmos_fields!(csf, model_sims, boundary_space, turbulent_fluxes) (; atmos_sim) = model_sims - # from atmos - if turbulent_fluxes == FluxCalculator.CombinedAtmosGrid() + # turbulent fluxes + if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes Regridder.dummmy_remap!(csf.F_turb_energy, Interfacer.get_field(atmos_sim, Val(:turbulent_energy_flux))) Regridder.dummmy_remap!(csf.F_turb_moisture, Interfacer.get_field(atmos_sim, Val(:turbulent_moisture_flux))) - - # surface density is needed for q_sat and requires atmos and sfc states, so it is calculated and saved in the coupler - Regridder.dummmy_remap!(csf.ρ_sfc, FluxCalculator.calculate_surface_air_density(atmos_sim, csf.T_S)) end + # surface density - needed for q_sat and requires atmos and sfc states, so it is calculated and saved in the coupler + Regridder.dummmy_remap!(csf.ρ_sfc, FluxCalculator.calculate_surface_air_density(atmos_sim, csf.T_S)) + + # radiative fluxes Regridder.dummmy_remap!(csf.F_radiative, Interfacer.get_field(atmos_sim, Val(:radiative_energy_flux))) + + # precipitation Regridder.dummmy_remap!(csf.P_liq, Interfacer.get_field(atmos_sim, Val(:liquid_precipitation))) - Regridder.dummmy_remap!(csf.P_snow, Interfacer.get_field(atmos_sim, Val(:snow_precipitation))) # generalize + Regridder.dummmy_remap!(csf.P_snow, Interfacer.get_field(atmos_sim, Val(:snow_precipitation))) end """ @@ -67,7 +70,7 @@ function import_combined_surface_fields!(csf, model_sims, boundary_space, turbul Regridder.combine_surfaces!(combined_field, model_sims, Val(:albedo)) Regridder.dummmy_remap!(csf.albedo, combined_field) - if turbulent_fluxes !== FluxCalculator.PartitionedComponentModelGrid() + if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes Regridder.combine_surfaces!(combined_field, model_sims, Val(:roughness_momentum)) Regridder.dummmy_remap!(csf.z0m_S, combined_field) @@ -97,12 +100,15 @@ function update_sim!(atmos_sim::Interfacer.AtmosModelSimulation, csf, turbulent_ Interfacer.update_field!(atmos_sim, Val(:albedo), csf.albedo) Interfacer.update_field!(atmos_sim, Val(:surface_temperature), csf.T_S) - if turbulent_fluxes == FluxCalculator.CombinedAtmosGrid() + if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes Interfacer.update_field!(atmos_sim, Val(:roughness_momentum), csf.z0m_S) Interfacer.update_field!(atmos_sim, Val(:roughness_buoyancy), csf.z0b_S) Interfacer.update_field!(atmos_sim, Val(:beta), csf.beta) # not in this version of atmos end + if turbulent_fluxes isa FluxCalculator.PartitionedStateFluxes + Interfacer.update_field!(atmos_sim, Val(:turbulent_fluxes), csf) + end end """ @@ -114,29 +120,37 @@ Updates the surface component model cache with the current coupler fields of F_t - `sim`: [Interfacer.SurfaceModelSimulation] containing a surface model simulation object. - `csf`: [NamedTuple] containing coupler fields. """ -function update_sim!(sim::Interfacer.SurfaceModelSimulation, csf, area_fraction = nothing) +function update_sim!(sim::Interfacer.SurfaceModelSimulation, csf, turbulent_fluxes, area_fraction = nothing) FT = eltype(area_fraction) - Interfacer.update_field!(sim, Val(:air_density), csf.ρ_sfc) - Interfacer.update_field!( - sim, - Val(:turbulent_energy_flux), - Regridder.binary_mask.(area_fraction, threshold = eps(FT)) .* csf.F_turb_energy, - ) + # atmospheric surface density + Interfacer.update_field!(sim, Val(:air_density), csf.ρ_sfc) - Interfacer.update_field!( - sim, - Val(:turbulent_moisture_flux), - Regridder.binary_mask.(area_fraction, threshold = eps(FT)) .* csf.F_turb_moisture, - ) + # turbulent fluxes + # when PartitionedStateFluxes, turbulent fluxes are updated during the flux calculation + if turbulent_fluxes isa FluxCalculator.CombinedStateFluxes + + Interfacer.update_field!( + sim, + Val(:turbulent_energy_flux), + Regridder.binary_mask.(area_fraction, threshold = eps(FT)) .* csf.F_turb_energy, + ) + Interfacer.update_field!( + sim, + Val(:turbulent_moisture_flux), + Regridder.binary_mask.(area_fraction, threshold = eps(FT)) .* csf.F_turb_moisture, + ) + end + # radiative fluxes Interfacer.update_field!( sim, Val(:radiative_energy_flux), Regridder.binary_mask.(area_fraction, threshold = eps(FT)) .* csf.F_radiative, ) + # precipitation Interfacer.update_field!(sim, Val(:liquid_precipitation), .-csf.P_liq) Interfacer.update_field!(sim, Val(:snow_precipitation), .-csf.P_snow) end @@ -163,7 +177,7 @@ Iterates `update_sim!` over all component model simulations saved in `cs.model_s function update_model_sims!(model_sims, csf, turbulent_fluxes) for sim in model_sims if sim isa Interfacer.SurfaceModelSimulation - update_sim!(sim, csf, Interfacer.get_field(sim, Val(:area_fraction))) + update_sim!(sim, csf, turbulent_fluxes, Interfacer.get_field(sim, Val(:area_fraction))) else update_sim!(sim, csf, turbulent_fluxes) end diff --git a/src/FluxCalculator.jl b/src/FluxCalculator.jl index 5a9230ae7..88118899e 100644 --- a/src/FluxCalculator.jl +++ b/src/FluxCalculator.jl @@ -6,15 +6,24 @@ or to call flux calculating functions from the component models. """ module FluxCalculator +import SurfaceFluxes as SF +import Thermodynamics as TD +using StaticArrays -using ClimaCoupler.Interfacer -using ClimaCore.Fields -export PartitionedComponentModelGrid, - CombinedAtmosGrid, - compute_combined_turbulent_fluxes!, +using ClimaCoupler: Interfacer, Regridder +using ClimaCore: Fields, Spaces +export PartitionedStateFluxes, + CombinedStateFluxes, + combined_turbulent_fluxes!, TurbulentFluxPartition, - compute_atmos_turbulent_fluxes!, - calculate_surface_air_density + atmos_turbulent_fluxes!, + calculate_surface_air_density, + extrapolate_ρ_to_sfc, + MoninObukhovScheme, + BulkScheme, + partitioned_turbulent_fluxes!, + get_surface_params, + update_turbulent_fluxes_point! """ TurbulentFluxPartition @@ -24,23 +33,23 @@ Abstract type for flags that denote where and how to calculate tubulent fluxes. abstract type TurbulentFluxPartition end """ - PartitionedComponentModelGrid <: TurbulentFluxPartition + PartitionedStateFluxes <: TurbulentFluxPartition A flag indicating that the turbulent fluxes should be partitioned and calculated over each surface model and then combined. This is calculated on the coupler grid. """ -struct PartitionedComponentModelGrid <: TurbulentFluxPartition end +struct PartitionedStateFluxes <: TurbulentFluxPartition end """ - CombinedAtmosGrid <: TurbulentFluxPartition + CombinedStateFluxes <: TurbulentFluxPartition A flag indicating that the turbulent fluxes (e.g. sensible and latent heat fluxes, drag and moisture fluxes) are to be calculated on the Atmos grid, and saved in Atmos cache. """ -struct CombinedAtmosGrid <: TurbulentFluxPartition end +struct CombinedStateFluxes <: TurbulentFluxPartition end """ - compute_combined_turbulent_fluxes!(model_sims, csf, turbulent_fluxes::TurbulentFluxPartition) + combined_turbulent_fluxes!(model_sims, csf, turbulent_fluxes::TurbulentFluxPartition) Calls the method(s) which calculate turbulent surface fluxes from combined surface states in coupler fields, `csf`. @@ -50,16 +59,16 @@ Calls the method(s) which calculate turbulent surface fluxes from combined surfa - `turbulent_fluxes`: [TurbulentFluxPartition] denotes a flag for turbulent flux calculation. """ -function compute_combined_turbulent_fluxes!(model_sims, csf, turbulent_fluxes::TurbulentFluxPartition) - if turbulent_fluxes == CombinedAtmosGrid() - compute_atmos_turbulent_fluxes!(model_sims.atmos_sim, csf) +function combined_turbulent_fluxes!(model_sims, csf, turbulent_fluxes::TurbulentFluxPartition) + if turbulent_fluxes isa CombinedStateFluxes + atmos_turbulent_fluxes!(model_sims.atmos_sim, csf) else nothing # TODO: may want to add CombinedCouplerGrid end end """ - compute_atmos_turbulent_fluxes!(sim::Interfacer.ComponentModelSimulation, csf) + atmos_turbulent_fluxes!(sim::Interfacer.ComponentModelSimulation, csf) A function to calculate turbulent surface fluxes using the combined surface states. It is required that a method is defined for the given `sim` and that the fluxes are @@ -72,13 +81,13 @@ saved in that sim's cache. `csf` refers to the coupler fields. # Example: ``` -function compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) +function atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf) atmos_sim.cache.flux .= atmos_sim.c .* (csf.T_S .- atmos_sim.temperature) end ``` """ -compute_atmos_turbulent_fluxes!(sim::Interfacer.ComponentModelSimulation, _) = +atmos_turbulent_fluxes!(sim::Interfacer.ComponentModelSimulation, _) = error("calling flux calculation in " * Interfacer.name(sim) * " but no method defined") @@ -91,4 +100,273 @@ function calculate_surface_air_density(atmos_sim::Interfacer.AtmosModelSimulatio error("this function is required to be dispatched on" * Interfacer.name(atmos_sim) * ", but no method defined") end +""" + partitioned_turbulent_fluxes!(model_sims::NamedTuple, fields::NamedTuple, boundary_space::Spaces.AbstractSpace, surface_scheme, thermo_params::TD.Parameters.ThermodynamicsParameters) + +The current setup calculates the aerodynamic fluxes in the coupler (assuming no regridding is needed) +using adapter function `get_surface_fluxes_point!`, which calls `SurfaceFluxes.jl`. The coupler saves +the area-weighted sums of the fluxes. + +Args: +- `model_sims`: [NamedTuple] containing `ComponentModelSimulation`s. +- `fields`: [NamedTuple] containing coupler fields. +- `boundary_space`: [Spaces.AbstractSpace] the space of the coupler surface. +- `surface_scheme`: [AbstractSurfaceFluxScheme] the surface flux scheme. +- `thermo_params`: [TD.Parameters.ThermodynamicsParameters] the thermodynamic parameters. + +TODO: +- generalize interface for regridding and take land state out of atmos's integrator.p +- add flux accumulation +- add flux bounds + +(NB: Radiation surface fluxes are calculated by the atmosphere.) + +""" +function partitioned_turbulent_fluxes!( + model_sims::NamedTuple, + fields::NamedTuple, + boundary_space::Spaces.AbstractSpace, + surface_scheme, + thermo_params::TD.Parameters.ThermodynamicsParameters, +) + + atmos_sim = model_sims.atmos_sim + csf = fields + FT = eltype(csf[1]) + + # reset coupler fields + csf.F_turb_ρτxz .*= FT(0) + csf.F_turb_ρτyz .*= FT(0) + csf.F_turb_energy .*= FT(0) + csf.F_turb_moisture .*= FT(0) + + # iterate over all columns (when regridding, this will need to change) + Fields.bycolumn(boundary_space) do colidx + # atmos state of center level 1 + z_int = Interfacer.get_field(atmos_sim, Val(:height_int), colidx) + uₕ_int = Interfacer.get_field(atmos_sim, Val(:uv_int), colidx) + thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int), colidx) + z_sfc = Interfacer.get_field(atmos_sim, Val(:height_sfc), colidx) + + for sim in model_sims + # iterate over all surface models with non-zero area fractions + if sim isa Interfacer.SurfaceModelSimulation + # get area fraction (min = 0, max = 1) + area_fraction = Interfacer.get_field(sim, Val(:area_fraction), colidx) + # get area mask [0, 1], where area_mask = 1 if area_fraction > 0 + area_mask = Regridder.binary_mask.(area_fraction, threshold = eps()) + + if !iszero(parent(area_mask)) + + thermo_state_sfc = surface_thermo_state(sim, thermo_params, thermo_state_int, colidx) + + # set inputs based on whether the surface_scheme is `MoninObukhovScheme` or `BulkScheme` + inputs = surface_inputs( + surface_scheme, + thermo_state_sfc, + thermo_state_int, + uₕ_int, + z_int, + z_sfc, + get_scheme_properties(surface_scheme, sim, colidx), + ) + + # update fluxes in the coupler + surface_params = get_surface_params(atmos_sim) + F_turb_ρτxz, F_turb_ρτyz, F_shf, F_lhf, F_turb_moisture = + get_surface_fluxes_point!(inputs, surface_params) + + fields = (; + F_turb_ρτxz = F_turb_ρτxz, + F_turb_ρτyz = F_turb_ρτyz, + F_turb_energy = F_shf .+ F_lhf, + F_turb_moisture = F_turb_moisture, + ) + + # update the fluxes of this surface model + update_turbulent_fluxes_point!(sim, fields, colidx) + + # add the flux contributing from this surface + @. csf.F_turb_ρτxz[colidx] += F_turb_ρτxz * area_fraction * area_mask + @. csf.F_turb_ρτyz[colidx] += F_turb_ρτyz * area_fraction * area_mask + @. csf.F_turb_energy[colidx] += (F_shf .+ F_lhf) * area_fraction * area_mask + @. csf.F_turb_moisture[colidx] += F_turb_moisture * area_fraction * area_mask + + end + end + end + + end + + # TODO: add allowable bounds here, check explicitly that all fluxes are equal + +end + +abstract type AbstractSurfaceFluxScheme end +struct BulkScheme <: AbstractSurfaceFluxScheme end +struct MoninObukhovScheme <: AbstractSurfaceFluxScheme end + +""" + get_scheme_properties(scheme::AbstractSurfaceFluxScheme, sim::Interfacer.SurfaceModelSimulation, colidx::Fields.ColumnIndex) + +Returns the scheme-specific properties for the surface model simulation `sim`. +""" +function get_scheme_properties(::BulkScheme, sim::Interfacer.SurfaceModelSimulation, colidx::Fields.ColumnIndex) + Ch = Interfacer.get_field(sim, Val(:heat_transfer_coefficient), colidx) + Cd = Interfacer.get_field(sim, Val(:beta), colidx) + beta = Interfacer.get_field(sim, Val(:drag_coefficient), colidx) + FT = eltype(Ch) + return (; z0b = FT(0), z0m = FT(0), Ch = Ch, Cd = Cd, beta = beta, gustiness = FT(1)) +end +function get_scheme_properties(::MoninObukhovScheme, sim::Interfacer.SurfaceModelSimulation, colidx::Fields.ColumnIndex) + z0m = Interfacer.get_field(sim, Val(:roughness_momentum), colidx) + z0b = Interfacer.get_field(sim, Val(:roughness_buoyancy), colidx) + beta = Interfacer.get_field(sim, Val(:beta), colidx) + FT = eltype(z0m) + return (; z0b = z0b, z0m = z0m, Ch = FT(0), Cd = FT(0), beta = beta, gustiness = FT(1)) +end + +""" + surface_inputs(scheme::AbstractSurfaceFluxScheme, thermo_state_sfc, thermo_state_int, uₕ_int, z_int, z_sfc, z0b, z0m, Ch, Cd, beta, gustiness) + +Returns the inputs for the surface model simulation `sim`. +""" +function surface_inputs(::BulkScheme, thermo_state_sfc, thermo_state_int, uₕ_int, z_int, z_sfc, surface_properties) + FT = Spaces.undertype(axes(z_sfc)) + + (; z0b, z0m, Ch, Cd, beta, gustiness) = surface_properties + # wrap state values + return @. SF.Coefficients( + SF.InteriorValues(z_int, uₕ_int, thermo_state_int), # state_in + SF.SurfaceValues( # state_sfc + z_sfc, + StaticArrays.SVector(FT(0), FT(0)), + thermo_state_sfc, + ), + Cd, # Cd + Ch, # Ch + z0m, # roughness_momentum + z0b, # roughness_buoyancy + gustiness, # gustiness + beta, # beta + ) +end +function surface_inputs( + ::MoninObukhovScheme, + thermo_state_sfc, + thermo_state_int, + uₕ_int, + z_int, + z_sfc, + surface_properties, +) + FT = Spaces.undertype(axes(z_sfc)) + + (; z0b, z0m, Ch, Cd, beta, gustiness) = surface_properties + + # wrap state values + return @. SF.ValuesOnly( + SF.InteriorValues(z_int, uₕ_int, thermo_state_int), # state_in + SF.SurfaceValues( # state_sfc + z_sfc, + StaticArrays.SVector(FT(0), FT(0)), + thermo_state_sfc, + ), + z0m, # z0m + z0b, # z0b + gustiness, # gustiness + beta, # beta + ) + +end + +""" + surface_thermo_state(sim::Interfacer.SurfaceModelSimulation, thermo_params::TD.Parameters.ThermodynamicsParameters, thermo_state_int, colidx::Fields.ColumnIndex) + +Returns the surface parameters for the surface model simulation `sim`. The default is assuming saturated surfaces, unless an extension is defined for the given `SurfaceModelSimulation`. +""" +function surface_thermo_state( + sim::Interfacer.SurfaceModelSimulation, + thermo_params::TD.Parameters.ThermodynamicsParameters, + thermo_state_int, + colidx::Fields.ColumnIndex, +) + @warn("Simulation " * Interfacer.name(sim) * " uses the default thermo (saturated) surface state", maxlog = 10) + T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx) # + ρ_sfc = extrapolate_ρ_to_sfc.(thermo_params, thermo_state_int, T_sfc) # ideally the # calculate elsewhere, here just getter... + q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) # default: saturated liquid surface + @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc) +end + +# TODO: (an equivalent of) this function also lives in Atmos and Land - should move to general utilities +""" + extrapolate_ρ_to_sfc(thermo_params, ts_int, T_sfc) + +Uses the ideal gas law and hydrostatic balance to extrapolate for surface density. +""" +function extrapolate_ρ_to_sfc(thermo_params, ts_in, T_sfc) + T_int = TD.air_temperature(thermo_params, ts_in) + Rm_int = TD.gas_constant_air(thermo_params, ts_in) + ρ_air = TD.air_density(thermo_params, ts_in) + ρ_air * (T_sfc / T_int)^(TD.cv_m(thermo_params, ts_in) / Rm_int) +end + +""" + get_surface_fluxes_point!(inputs, surface_params::SF.Parameters.SurfaceFluxesParameters) + +Uses SurfaceFluxes.jl to calculate turbulent surface fluxes. It should be atmos model agnostic, and columnwise. +""" +function get_surface_fluxes_point!(inputs, surface_params::SF.Parameters.SurfaceFluxesParameters) + + # calculate all fluxes (saturated surface conditions) + outputs = @. SF.surface_conditions(surface_params, inputs) + + # drag + F_turb_ρτxz = outputs.ρτxz + F_turb_ρτyz = outputs.ρτyz + + # energy fluxes + F_shf = outputs.shf + F_lhf = outputs.lhf + + # moisture + F_turb_moisture = @. SF.evaporation(surface_params, inputs, outputs.Ch) + + return F_turb_ρτxz, F_turb_ρτyz, F_shf, F_lhf, F_turb_moisture +end + +""" + get_surface_params(atmos_sim::Interfacer.AtmosModelSimulation) + +Returns the surface parameters of type `SF.Parameters.SurfaceFluxesParameters`. + +TODO: in the future this may not need to depend on the atmos sim, but +here retaining the dependency until we know how EDMF boundary conditions will +be handled (for consistency of parameters). +""" +function get_surface_params(atmos_sim::Interfacer.AtmosModelSimulation) + return error( + "get_surface_params is required to be dispatched on" * Interfacer.name(atmos_sim) * ", but no method defined", + ) +end + +""" + update_turbulent_fluxes_point!(sim::Interfacer.SurfaceModelSimulation, fields::NamedTuple, colidx::Fields.ColumnIndex) + +Updates the fluxes in the surface model simulation `sim` with the fluxes in `fields`. +""" +function update_turbulent_fluxes_point!( + sim::Interfacer.SurfaceModelSimulation, + fields::NamedTuple, + colidx::Fields.ColumnIndex, +) + return error( + "update_turbulent_fluxes_point! is required to be dispatched on" * + Interfacer.name(sim) * + ", but no method defined", + ) +end + +update_turbulent_fluxes_point!(sim::Interfacer.SurfaceStub, fields::NamedTuple, colidx::Fields.ColumnIndex) = nothing + end # module diff --git a/src/Interfacer.jl b/src/Interfacer.jl index 6a68123e7..c1fa427bf 100644 --- a/src/Interfacer.jl +++ b/src/Interfacer.jl @@ -6,7 +6,6 @@ This modules contains abstract types, interface templates and model stubs for co module Interfacer import Thermodynamics as TD - using ClimaCore: Fields export ComponentModelSimulation, AtmosModelSimulation, diff --git a/test/Project.toml b/test/Project.toml index 6d88bc711..44c63656d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -14,7 +14,9 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" @@ -27,5 +29,7 @@ IntervalSets = "0.7" MPI = "0.20" NCDatasets = "0.11 - 0.12.16" Pkg = "1" +StaticArrays = "1.6" +SurfaceFluxes = "0.6.3" Unitful = "1" julia = "1.8.5" diff --git a/test/field_exchanger_tests.jl b/test/field_exchanger_tests.jl index 14cf2f469..0fff9a07c 100644 --- a/test/field_exchanger_tests.jl +++ b/test/field_exchanger_tests.jl @@ -3,7 +3,7 @@ using ClimaCore: Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaCoupler: Utilities, Regridder, TestHelper import ClimaCoupler.Interfacer: update_field!, AtmosModelSimulation, SurfaceModelSimulation, SurfaceStub, get_field, update_field! -import ClimaCoupler.FluxCalculator: CombinedAtmosGrid, PartitionedComponentModelGrid, calculate_surface_air_density +import ClimaCoupler.FluxCalculator: CombinedStateFluxes, PartitionedStateFluxes, calculate_surface_air_density import ClimaCoupler.FieldExchanger: import_atmos_fields!, import_combined_surface_fields!, @@ -46,7 +46,7 @@ end model_sims = (; atmos_sim = DummySimulation(atmos_fields)) - flux_types = (CombinedAtmosGrid(), PartitionedComponentModelGrid()) + flux_types = (CombinedStateFluxes(), PartitionedStateFluxes()) results = [FT(1), FT(0)] for (i, t) in enumerate(flux_types) coupler_fields = NamedTuple{coupler_names}(ntuple(i -> Fields.zeros(boundary_space), length(coupler_names))) @@ -93,9 +93,9 @@ get_field(sim::Union{TestSurfaceSimulation2, TestSurfaceSimulation2}, ::Val{:sur sims = (; a = TestSurfaceSimulation1(ones(boundary_space)), b = TestSurfaceSimulation2(ones(boundary_space))) - # test the coupler update under CombinedAtmosGrid (update all) and PartitionedComponentModelGrid (update all except + # test the coupler update under CombinedStateFluxes (update all) and PartitionedStateFluxes (update all except # surface info needed for the calculation of turbulent fluxes) - flux_types = (CombinedAtmosGrid(), PartitionedComponentModelGrid()) + flux_types = (CombinedStateFluxes(), PartitionedStateFluxes()) results = [FT(0.5), FT(0)] # 0.5 due to the area fraction weighting for (i, t) in enumerate(flux_types) coupler_fields = NamedTuple{coupler_names}(ntuple(i -> Fields.zeros(boundary_space), length(coupler_names))) @@ -159,8 +159,8 @@ end stub_sim = SurfaceStub((; area_fraction = Fields.ones(boundary_space), ρ_sfc = Fields.ones(boundary_space))), ) - # test the sim update under CombinedAtmosGrid (update all) and PartitionedComponentModelGrid (update all except turbulent fluxes) - flux_types = (CombinedAtmosGrid(), PartitionedComponentModelGrid()) + # test the sim update under CombinedStateFluxes (update all) and PartitionedStateFluxes (update all except turbulent fluxes) + flux_types = (CombinedStateFluxes(), PartitionedStateFluxes()) results = [FT(0), FT(1)] for (i, t) in enumerate(flux_types) model_sims.atmos_sim.cache.roughness_momentum .= FT(0) @@ -168,7 +168,7 @@ end # test atmos @test parent(model_sims.atmos_sim.cache.albedo)[1] == results[2] - if t == CombinedAtmosGrid() + if t isa CombinedStateFluxes @test parent(model_sims.atmos_sim.cache.roughness_momentum)[1] == results[2] else @test parent(model_sims.atmos_sim.cache.roughness_momentum)[1] == results[1] diff --git a/test/flux_calculator_tests.jl b/test/flux_calculator_tests.jl index 460f3bfbc..2f93f9f9b 100644 --- a/test/flux_calculator_tests.jl +++ b/test/flux_calculator_tests.jl @@ -3,44 +3,58 @@ using ClimaCore: Meshes, Domains, Topologies, Spaces, Fields, InputOutput using ClimaCoupler: Utilities, Regridder, TestHelper using Test import ClimaCoupler.FluxCalculator: - compute_atmos_turbulent_fluxes!, - compute_combined_turbulent_fluxes!, - CombinedAtmosGrid, - PartitionedComponentModelGrid, - calculate_surface_air_density -import ClimaCoupler.Interfacer: AtmosModelSimulation + atmos_turbulent_fluxes!, + combined_turbulent_fluxes!, + CombinedStateFluxes, + PartitionedStateFluxes, + calculate_surface_air_density, + MoninObukhovScheme, + BulkScheme, + partitioned_turbulent_fluxes!, + get_surface_params, + update_turbulent_fluxes_point!, + surface_thermo_state +import ClimaCoupler: Interfacer FT = Float64 -# test for a simple generic atmos model -struct DummySimulation{S, C} <: AtmosModelSimulation +import CLIMAParameters as CP +import Thermodynamics as TD +import Thermodynamics.Parameters as TP +import SurfaceFluxes as SF +import SurfaceFluxes.UniversalFunctions as UF + +using StaticArrays + +# # test for a simple generic atmos model +struct DummySimulation{S, C} <: Interfacer.AtmosModelSimulation state::S cache::C end -function compute_atmos_turbulent_fluxes!(sim::DummySimulation, csf) +function atmos_turbulent_fluxes!(sim::DummySimulation, csf) sim.cache.flux .= (csf.T_sfc .- sim.state.T) .* sim.cache.κ ./ sim.cache.dz # Eq. 1 end -struct DummySimulation2{C} <: AtmosModelSimulation +struct DummySimulation2{C} <: Interfacer.AtmosModelSimulation cache::C end -@testset "compute_combined_turbulent_fluxes!" begin +@testset "combined_turbulent_fluxes!" begin boundary_space = TestHelper.create_space(FT) coupler_fields = (; T_sfc = 310 .* ones(boundary_space)) sim = DummySimulation((; T = 300 .* ones(boundary_space)), (; κ = FT(0.01), dz = FT(1), flux = zeros(boundary_space))) model_sims = (; atmos_sim = sim) - flux_types = (CombinedAtmosGrid(), PartitionedComponentModelGrid()) - # the result of Eq 1 above, given these states, is 0.1 W/m2, but under PartitionedComponentModelGrid() turbulent fluxes are + flux_types = (CombinedStateFluxes(), PartitionedStateFluxes()) + # the result of Eq 1 above, given these states, is 0.1 W/m2, but under PartitionedStateFluxes() turbulent fluxes are # not calculated using this method (using combined surface properties), so the fluxes remain 0. results = [FT(0.1), FT(0.0)] for (i, t) in enumerate(flux_types) sim.cache.flux .= FT(0) - compute_combined_turbulent_fluxes!(model_sims, coupler_fields, t) + combined_turbulent_fluxes!(model_sims, coupler_fields, t) @test parent(sim.cache.flux)[1] == results[i] end sim2 = DummySimulation2((; cache = (; flux = zeros(boundary_space)))) model_sims = (; atmos_sim = sim2) - @test_throws ErrorException compute_atmos_turbulent_fluxes!(sim2, coupler_fields) + @test_throws ErrorException atmos_turbulent_fluxes!(sim2, coupler_fields) end @@ -50,3 +64,258 @@ end sim2 = DummySimulation2((; cache = (; flux = zeros(boundary_space)))) @test_throws ErrorException calculate_surface_air_density(sim2, coupler_fields.T_sfc) end + +# Test for the PartitionedStateFluxes() case +# parameter set generated CLIMAParamerers and helper functions from SurfaceFluxes.jl +function create_uf_parameters(toml_dict, ::UF.BusingerType) + FT = CP.float_type(toml_dict) + aliases = ["Pr_0_Businger", "a_m_Businger", "a_h_Businger", "ζ_a_Businger", "γ_Businger"] + + pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions") + pairs = (; pairs...) # convert to NamedTuple + + pairs = (; + Pr_0 = pairs.Pr_0_Businger, + a_m = pairs.a_m_Businger, + a_h = pairs.a_h_Businger, + ζ_a = pairs.ζ_a_Businger, + γ = pairs.γ_Businger, + ) + return UF.BusingerParams{FT}(; pairs...) +end + +function create_parameters(toml_dict, ufpt) + FT = CP.float_type(toml_dict) + + ufp = create_uf_parameters(toml_dict, ufpt) + AUFP = typeof(ufp) + + aliases = string.(fieldnames(TD.Parameters.ThermodynamicsParameters)) + pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics") + thermo_params = TD.Parameters.ThermodynamicsParameters{FT}(; pairs...) + TP = typeof(thermo_params) + + aliases = ["von_karman_const"] + pairs = CP.get_parameter_values!(toml_dict, aliases, "SurfaceFluxesParameters") + return SF.Parameters.SurfaceFluxesParameters{FT, AUFP, TP}(; pairs..., ufp, thermo_params) +end + +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +const sf_params = create_parameters(toml_dict, UF.BusingerType()) + +# atmos sim object and extensions +struct TestAtmos{P, Y, D, I} <: Interfacer.AtmosModelSimulation + params::P + Y_init::Y + domain::D + integrator::I +end +get_surface_params(::TestAtmos) = sf_params + +function get_thermo_params(sim::TestAtmos) + FT = sim.params.FT + aliases = string.(fieldnames(TP.ThermodynamicsParameters)) + toml_dict = CP.create_toml_dict(FT; dict_type = "alias") + pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics") + TP.ThermodynamicsParameters{FT}(; pairs...) +end + +Interfacer.get_field(sim::TestAtmos, ::Val{:height_int}) = sim.integrator.p.z +Interfacer.get_field(sim::TestAtmos, ::Val{:height_sfc}) = sim.integrator.p.z_sfc +Interfacer.get_field(sim::TestAtmos, ::Val{:uv_int}) = @. StaticArrays.SVector(sim.integrator.p.u, sim.integrator.p.v) +Interfacer.get_field(sim::TestAtmos, ::Val{:thermo_state_int}) = + TD.PhaseEquil_ρTq.(get_thermo_params(sim), sim.integrator.ρ, sim.integrator.T, sim.integrator.q) +Interfacer.get_field(sim::TestAtmos, ::Val{:air_density}) = sim.integrator.ρ +Interfacer.get_field(sim::TestAtmos, ::Val{:air_temperature}) = sim.integrator.T + +function update_sim!(sim::TestAtmos, fields, _) + (; F_turb_ρτxz, F_turb_energy, F_turb_moisture) = fields + ρ_int = sim.integrator.ρ + @. sim.integrator.p.energy_bc = -(F_turb_energy) + @. sim.integrator.p.ρq_tot_bc = -F_turb_moisture + @. sim.integrator.p.uₕ_bc = -(F_turb_ρτxz / ρ_int) # x-compoennt only for this test +end + +# ocean sim object and extensions +struct TestOcean{M, Y, D, I} <: Interfacer.SurfaceModelSimulation + model::M + Y_init::Y + domain::D + integrator::I +end +Interfacer.get_field(sim::TestOcean, ::Val{:surface_temperature}) = sim.integrator.T +Interfacer.get_field(sim::TestOcean, ::Val{:air_humidity}) = sim.integrator.p.q +Interfacer.get_field(sim::TestOcean, ::Val{:roughness_momentum}) = sim.integrator.p.z0m +Interfacer.get_field(sim::TestOcean, ::Val{:roughness_buoyancy}) = sim.integrator.p.z0b +Interfacer.get_field(sim::TestOcean, ::Val{:beta}) = sim.integrator.p.beta +Interfacer.get_field(sim::TestOcean, ::Val{:area_fraction}) = sim.integrator.p.area_fraction +Interfacer.get_field(sim::TestOcean, ::Val{:heat_transfer_coefficient}) = sim.integrator.p.Ch +Interfacer.get_field(sim::TestOcean, ::Val{:drag_coefficient}) = sim.integrator.p.Cd +Interfacer.get_field(sim::TestOcean, ::Val{:albedo}) = sim.integrator.p.α + +function surface_thermo_state( + sim::TestOcean, + thermo_params::TD.Parameters.ThermodynamicsParameters, + thermo_state_int, + colidx::Fields.ColumnIndex, +) + T_sfc = Interfacer.get_field(sim, Val(:surface_temperature), colidx) + ρ_sfc = thermo_state_int.ρ # arbitrary + q_sfc = Interfacer.get_field(sim, Val(:air_humidity), colidx) # read from cache + @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc) +end + +function update_turbulent_fluxes_point!(sim::TestOcean, fields::NamedTuple, colidx::Fields.ColumnIndex) + (; F_turb_energy) = fields + @. sim.integrator.p.F_aero[colidx] = F_turb_energy +end + +@testset "calculate correct fluxes: dry" begin + surface_scheme_list = (MoninObukhovScheme(), BulkScheme()) + for scheme in surface_scheme_list + boundary_space = TestHelper.create_space(FT) + + params = (; surface_scheme = scheme, FT = FT) + + # atmos + p = (; + energy_bc = zeros(boundary_space), + ρq_tot_bc = zeros(boundary_space), + uₕ_bc = ones(boundary_space), + z = ones(boundary_space), + z_sfc = zeros(boundary_space), + u = ones(boundary_space), + v = ones(boundary_space), + ) + Y_init = (; ρ = ones(boundary_space) .* 1.2, T = ones(boundary_space) .* 310, q = zeros(boundary_space)) + integrator = (; Y_init..., p = p) + atmos_sim = TestAtmos(params, Y_init, nothing, integrator) + + # ocean + p = (; + F_aero = zeros(boundary_space), + z0m = FT(0.01), + z0b = FT(0.01), + beta = ones(boundary_space), + α = ones(boundary_space) .* 0.5, + q = zeros(boundary_space), + Cd = FT(0.01), + Ch = FT(0.01), + area_fraction = ones(boundary_space) .* 0.5, + ) + Y_init = (; T = ones(boundary_space) .* 300.0) + integrator = (; Y_init..., p = p) + ocean_sim = TestOcean(nothing, Y_init, nothing, integrator) + + # ocean + ocean_sim2 = TestOcean(nothing, Y_init, nothing, integrator) + + model_sims = (; atmos_sim, ocean_sim, ocean_sim2) + + coupler_cache_names = ( + :T_S, + :albedo, + :F_R_sfc, + :F_R_toa, + :P_liq, + :P_snow, + :P_net, + :F_turb_energy, + :F_turb_ρτxz, + :F_turb_ρτyz, + :F_turb_moisture, + ) + fields = NamedTuple{coupler_cache_names}(ntuple(i -> Fields.zeros(boundary_space), length(coupler_cache_names))) + + # calculate turbulent fluxes + thermo_params = get_thermo_params(atmos_sim) + partitioned_turbulent_fluxes!(model_sims, fields, boundary_space, scheme, thermo_params) + + # calculating the fluxes twice ensures that no accumulation occurred (i.e. fluxes are reset to zero each time) + # TODO: this will need to be extended once flux accumulation is re-enabled + partitioned_turbulent_fluxes!(model_sims, fields, boundary_space, scheme, thermo_params) + + windspeed = @. hypot(atmos_sim.integrator.p.u, atmos_sim.integrator.p.v) + + thermo_params = get_thermo_params(atmos_sim) + thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int)) + + surface_thermo_states = similar(thermo_state_int) + Fields.bycolumn(boundary_space) do colidx + surface_thermo_states[colidx] .= + surface_thermo_state(ocean_sim, thermo_params, thermo_state_int[colidx], colidx) + end + + # analytical solution is possible for the BulkScheme() case + if scheme isa BulkScheme + ρ_sfc = Interfacer.get_field(atmos_sim, Val(:air_density)) + cpm = TD.cv_m.(thermo_params, thermo_state_int) .+ TD.gas_constant_air.(thermo_params, thermo_state_int) # cp = R + cv + gz = + ( + Interfacer.get_field(atmos_sim, Val(:height_int)) .- + Interfacer.get_field(atmos_sim, Val(:height_sfc)) + ) .* 9.81 + shf_analytical = @. (cpm * (ocean_sim.integrator.T - atmos_sim.integrator.T) - gz) * + ocean_sim.integrator.p.Ch * + ρ_sfc * + windspeed #-ρ_sfc * Ch * windspeed(sc) * (cp_m * ΔT + ΔΦ) + + colidx = Fields.ColumnIndex{2}((1, 1), 73) # arbitrary index + # check the coupler field update + @test isapprox(parent(shf_analytical[colidx]), parent(fields.F_turb_energy[colidx]), rtol = 1e-6) + + # test the surface field update + @test parent(fields.F_turb_energy[colidx]) == parent(ocean_sim.integrator.p.F_aero[colidx]) + + # test the atmos field update + update_sim!(atmos_sim, fields, nothing) + @test parent(fields.F_turb_energy[colidx]) == -parent(atmos_sim.integrator.p.energy_bc[colidx]) + + end + @test parent(fields.F_turb_moisture)[1] ≈ FT(0) + + end +end +@testset "get_surface_params" begin + @test get_surface_params(TestAtmos([], [], [], [])) == sf_params + sim = DummySimulation([], []) + @test_throws ErrorException( + "get_surface_params is required to be dispatched on" * Interfacer.name(sim) * ", but no method defined", + ) get_surface_params(DummySimulation([], [])) +end + +struct DummySurfaceSimulation3{M, Y, D, I} <: Interfacer.SurfaceModelSimulation + model::M + Y_init::Y + domain::D + integrator::I +end +Interfacer.get_field(sim::DummySurfaceSimulation3, ::Val{:surface_temperature}) = sim.integrator.T + +@testset "update_turbulent_fluxes_point!" begin + sim = Interfacer.SurfaceStub([]) + sim2 = DummySurfaceSimulation3([], [], [], []) + colidx = Fields.ColumnIndex{2}((1, 1), 73) # arbitrary index + @test update_turbulent_fluxes_point!(sim, (;), colidx) == nothing + @test_throws ErrorException( + "update_turbulent_fluxes_point! is required to be dispatched on" * + Interfacer.name(sim2) * + ", but no method defined", + ) update_turbulent_fluxes_point!(sim2, (;), colidx) == ErrorException +end + +@testset "surface_thermo_state" begin + boundary_space = TestHelper.create_space(FT) + _ones = Fields.ones(boundary_space) + sim = DummySurfaceSimulation3( + [], + [], + [], + (; T = _ones .* FT(300), ρ = _ones .* FT(1.2), p = (; q = _ones .* FT(0.01))), + ) + atmos_sim = TestAtmos((; FT = FT), [], [], (; T = _ones .* FT(300), ρ = _ones .* FT(1.2), q = _ones .* FT(0.01))) + thermo_params = get_thermo_params(atmos_sim) + thermo_state_int = Interfacer.get_field(atmos_sim, Val(:thermo_state_int)) + colidx = Fields.ColumnIndex{2}((1, 1), 73) # arbitrary index + @test surface_thermo_state(sim, thermo_params, thermo_state_int[colidx], colidx).ρ == thermo_state_int[colidx].ρ +end