use ClimaParams for shared parameters in coupled sims

Author: juliasloan25

experiments/ClimaEarth/components/atmosphere/climaatmos.jl

@@ -5,6 +5,7 @@ import LinearAlgebra
 import ClimaAtmos as CA
 import ClimaAtmos: set_surface_albedo!
 import ClimaAtmos.Parameters as CAP
+import ClimaParams as CP
 import ClimaCore as CC
 import ClimaCore.Geometry: ⊗
 import SurfaceFluxes as SF
@@ -611,14 +612,17 @@ The TOML parameter file to use is chosen using the following priority:
 If a coupler TOML file is provided, it is used.
 Otherwise we use an atmos TOML file if it's provided.
 If neither file is provided, the default ClimaAtmos parameters are used.
+This is analogous to the logic in `get_coupler_config_dict` for selecting
+the coupler configuration TOML file.
 """
-function get_atmos_config_dict(coupler_config::Dict, atmos_output_dir)
+function get_atmos_config_dict(
+    coupler_config::Dict,
+    atmos_output_dir,
+    coupled_param_dict::CP.ParamDict,
+    comms_ctx,
+)
     atmos_config = deepcopy(coupler_config)
 
-    # Rename parameters for consistency with ClimaAtmos.jl
-    # Rename atmosphere config file from ClimaCoupler convention to ClimaAtmos convention
-    atmos_config["config_file"] = coupler_config["atmos_config_file"]
-
     # Ensure Atmos's own checkpoints are synced up with ClimaCoupler, so that we
     # can pick up from where we have left. NOTE: This should not be needed, but
     # there is no easy way to initialize ClimaAtmos with a different t_start
@@ -631,18 +635,25 @@ function get_atmos_config_dict(coupler_config::Dict, atmos_output_dir)
         coupler_config["dt"]
     atmos_config["dt"] = dt_atmos
 
-    # Use atmos toml if coupler toml is not defined
-    # If we can't find the file at the relative path, prepend pkgdir(ClimaAtmos)
+    # Set up the atmosphere parameter dictionary (`toml_dict`)
+    # If we can't find the atmos TOML file at the relative path, prepend pkgdir(ClimaAtmos)
     atmos_toml = map(atmos_config["toml"]) do file
         isfile(file) ? file : joinpath(pkgdir(CA), file)
     end
+    # Use atmos toml only if coupler toml is not defined
     coupler_toml = atmos_config["coupler_toml"]
     toml = isempty(coupler_toml) ? atmos_toml : coupler_toml
     if !isempty(toml)
         @info "Overwriting Atmos parameters from input TOML file(s): $toml"
         atmos_config["toml"] = toml
     end
 
+    # Override atmos parameters with coupled parameters
+    FT = atmos_config["FLOAT_TYPE"] == "Float64" ? Float64 : Float32
+    override_file = CP.merge_toml_files(atmos_config["toml"], override = true)
+    atmos_toml = CP.create_toml_dict(FT; override_file) # TODO use ClimaParams 1.0.3
+    toml_dict = CP.merge_override_default_values(coupled_param_dict, atmos_toml)
+
     # Specify atmos output directory to be inside the coupler output directory
     atmos_config["output_dir_style"] = "RemovePreexisting"
     atmos_config["output_dir"] = atmos_output_dir
@@ -659,7 +670,25 @@ function get_atmos_config_dict(coupler_config::Dict, atmos_output_dir)
         atmos_config["restart_file"] =
             replace(atmos_config["restart_file"], "active" => "0000")
     end
-    return CA.AtmosConfig(atmos_config)
+
+    # Construct the AtmosConfig object
+    toml_dict = coupled_param_dict
+    config_files = atmos_config["toml"]
+    job_id = atmos_config["job_id"]
+    config = CA.override_default_config(atmos_config)
+
+    TD = typeof(toml_dict)
+    PA = typeof(config)
+    C = typeof(comms_ctx)
+    CF = typeof(config_files)
+    return CA.AtmosConfig{FT, TD, PA, C, CF}(
+        toml_dict,
+        config,
+        comms_ctx,
+        config_files,
+        job_id,
+    )
+    # return CA.AtmosConfig(atmos_config)
 end
 
 """

experiments/ClimaEarth/components/land/climaland_bucket.jl

@@ -66,16 +66,17 @@ function BucketSimulation(
     stepper = CTS.RK4(),
     albedo_type::String = "map_static",
     bucket_initial_condition::String = "",
-    energy_check::Bool = false,
-    parameter_files = [],
+    coupled_param_dict = CP.create_toml_dict(FT),
 ) where {FT, TT <: Union{Float64, ITime}}
-    # Construct the parameter dictionary based on the float type and any parameter files
-    toml_dict = LP.create_toml_dict(FT; override_files = parameter_files)
+    # Get default land parameters from ClimaLand.LandParameters
+    land_toml_dict = LP.create_toml_dict(FT)
+    # Override land parameters with coupled parameters
+    toml_dict = CP.merge_override_default_values(coupled_param_dict, land_toml_dict)
 
     # Note that this does not take into account topography of the surface, which is OK for this land model.
     # But it must be taken into account when computing surface fluxes, for Δz.
     if isnothing(shared_surface_space)
-        domain = make_land_domain(depth; nelements, dz_tuple)
+        domain = make_land_domain(depth, toml_dict; nelements, dz_tuple)
     else
         domain = make_land_domain(
             shared_surface_space,
@@ -379,11 +380,9 @@ function FieldExchanger.update_sim!(sim::BucketSimulation, csf)
     p = sim.integrator.p
     t = sim.integrator.t
 
-    # TODO: get sigma from parameters
-    FT = eltype(Y)
-    σ = FT(5.67e-8)
     # Note: here we add negative signs to account for a difference in sign convention
     #  between ClimaLand.jl and ClimaCoupler.jl in SW_d and LW_d.
+    σ = model.parameters.earth_param_set.Stefan
     sim.integrator.p.bucket.R_n .=
         .-(1 .- CL.surface_albedo(model, Y, p)) .* sim.radiative_fluxes.SW_d .-
         Interfacer.get_field(sim, Val(:emissivity)) .*

experiments/ClimaEarth/components/land/climaland_helpers.jl

@@ -1,5 +1,7 @@
 import ClimaCore as CC
 import ClimaLand as CL
+import ClimaParams as CP
+
 """
     temp_anomaly_aquaplanet(coord)
 
@@ -19,17 +21,22 @@ and result in stable simulations.
 temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4
 
 """
-    make_land_domain(nelements::Tuple{Int, Int}, depth::FT) where {FT}
-
+    make_land_domain(
+        depth::FT,
+        toml_dict::CP.ParamDict;
+        nelements::Tuple{Int, Int} = (101, 15),
+        dz_tuple::Tuple{FT, FT} = FT.((10.0, 0.05)),
+    ) where {FT}
 Creates a land model domain with the given number of vertical elements
 and vertical depth.
 """
 function make_land_domain(
-    depth::FT;
+    depth::FT,
+    toml_dict::CP.ParamDict;
     nelements::Tuple{Int, Int} = (101, 15),
     dz_tuple::Tuple{FT, FT} = FT.((10.0, 0.05)),
 ) where {FT}
-    radius = FT(6378.1e3)
+    radius = toml_dict["planet_radius"] # in meters
     npolynomial = 0
     domain = CL.Domains.SphericalShell(; radius, depth, nelements, npolynomial, dz_tuple)
     return domain

experiments/ClimaEarth/components/land/climaland_integrated.jl

@@ -86,13 +86,19 @@ function ClimaLandSimulation(
     atmos_h,
     land_temperature_anomaly::String = "amip",
     use_land_diagnostics::Bool = true,
-    parameter_files = [],
+    coupled_param_dict = CP.create_toml_dict(FT),
     land_ic_path::Union{Nothing, String} = nothing,
 ) where {FT, TT <: Union{Float64, ITime}}
+    # Get default land parameters from ClimaLand.LandParameters
+    land_toml_dict = LP.create_toml_dict(FT)
+    # Override land parameters with coupled parameters
+    toml_dict = CP.merge_override_default_values(coupled_param_dict, land_toml_dict)
+    earth_param_set = CL.Parameters.LandParameters(toml_dict)
+
     # Note that this does not take into account topography of the surface, which is OK for this land model.
     # But it must be taken into account when computing surface fluxes, for Δz.
     if isnothing(shared_surface_space)
-        domain = make_land_domain(depth; nelements, dz_tuple)
+        domain = make_land_domain(depth, toml_dict; nelements, dz_tuple)
     else
         domain = make_land_domain(
             shared_surface_space,
@@ -114,10 +120,6 @@ function ClimaLandSimulation(
     #  since that's where we compute fluxes for this land model
     atmos_h = Interfacer.remap(atmos_h, surface_space)
 
-    # Set up spatially-varying parameters
-    toml_dict = LP.create_toml_dict(FT; override_files = parameter_files)
-    earth_param_set = CL.Parameters.LandParameters(toml_dict)
-
     # Set up atmosphere and radiation forcing
     forcing = (;
         atmos = CL.CoupledAtmosphere{FT}(surface_space, atmos_h),

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -4,6 +4,7 @@ import ClimaCoupler: Checkpointer, FieldExchanger, FluxCalculator, Interfacer, U
 import ClimaComms
 import ClimaCore as CC
 import Thermodynamics as TD
+import ClimaParams as CP
 import ClimaOcean.EN4: download_dataset
 using KernelAbstractions: @kernel, @index, @inbounds
 
@@ -46,6 +47,7 @@ function OceananigansSimulation(
     stop_date;
     output_dir,
     comms_ctx = ClimaComms.context(),
+    coupled_param_dict = CP.create_toml_dict(eltype(area_fraction)),
 )
     arch = comms_ctx.device isa ClimaComms.CUDADevice ? OC.GPU() : OC.CPU()
 
@@ -154,10 +156,13 @@ function OceananigansSimulation(
 
     remapping = (; remapper_cc, scratch_cc1, scratch_cc2, scratch_arr1, scratch_arr2)
 
+    # Get some ocean properties and parameters
     ocean_properties = (;
         ocean_reference_density = 1020,
         ocean_heat_capacity = 3991,
         ocean_fresh_water_density = 999.8,
+        σ = coupled_param_dict["stefan_boltzmann_constant"],
+        C_to_K = coupled_param_dict["temperature_water_freeze"],
     )
 
     # Before version 0.96.22, the NetCDFWriter was broken on GPU
@@ -241,7 +246,7 @@ Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_diffuse_albedo}
 
 # NOTE: This is 3D, but it will be remapped to 2D
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_temperature}) =
-    273.15 + sim.ocean.model.tracers.T
+    sim.ocean.model.tracers.T + sim.ocean_properties.C_to_K # convert from Celsius to Kelvin
 
 """
     FluxCalculator.update_turbulent_fluxes!(sim::OceananigansSimulation, fields)
@@ -365,16 +370,15 @@ function FieldExchanger.update_sim!(sim::OceananigansSimulation, csf)
     # Update only the part due to radiative fluxes. For the full update, the component due
     # to latent and sensible heat is missing and will be updated in update_turbulent_fluxes.
     oc_flux_T = surface_flux(sim.ocean.model.tracers.T)
-    # TODO: get sigma from parameters
-    σ = 5.67e-8
+    (; σ, C_to_K) = sim.ocean_properties
     α = Interfacer.get_field(sim, Val(:surface_direct_albedo)) # scalar
     ϵ = Interfacer.get_field(sim, Val(:emissivity)) # scalar
     OC.interior(oc_flux_T, :, :, 1) .=
         (
             -(1 - α) .* remapped_SW_d .-
             ϵ * (
                 remapped_LW_d .-
-                σ .* (273.15 .+ OC.interior(sim.ocean.model.tracers.T, :, :, 1)) .^ 4
+                σ .* (C_to_K .+ OC.interior(sim.ocean.model.tracers.T, :, :, 1)) .^ 4
             )
         ) ./ (ocean_reference_density * ocean_heat_capacity)
 

experiments/ClimaEarth/components/ocean/prescr_ocean.jl

@@ -66,6 +66,7 @@ function PrescribedOceanSimulation(
     start_date,
     t_start,
     area_fraction,
+    coupled_param_dict,
     thermo_params,
     comms_ctx;
     z0m = FT(5.8e-5),
@@ -92,12 +93,13 @@ function PrescribedOceanSimulation(
         end : sst_path
     @info "PrescribedOcean: using SST file" sst_data
 
+    C_to_K = coupled_param_dict["temperature_water_freeze"]
     SST_timevaryinginput = TimeVaryingInput(
         sst_data,
         "SST",
         space,
         reference_date = start_date,
-        file_reader_kwargs = (; preprocess_func = (data) -> data + FT(273.15),), ## convert to Kelvin
+        file_reader_kwargs = (; preprocess_func = (data) -> data + C_to_K,), ## convert Celsius to Kelvin
     )
 
     SST_init = zeros(space)

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -34,16 +34,68 @@ end
 
 # sea-ice parameters
 Base.@kwdef struct IceSlabParameters{FT <: AbstractFloat}
-    h::FT = 2               # ice thickness [m]
-    ρ::FT = 900             # density of sea ice [kg / m3]
-    c::FT = 2100            # specific heat of sea ice [J / kg / K]
-    T_base::FT = 271.2      # temperature of sea water at the ice base
-    z0m::FT = 1e-4          # roughness length for momentum [m]
-    z0b::FT = 1e-4          # roughness length for tracers [m]
-    T_freeze::FT = 271.2    # freezing temperature of sea water [K]
-    k_ice::FT = 2           # thermal conductivity of sea ice [W / m / K] (less in HM71)
-    α::FT = 0.65            # albedo of sea ice, roughly tuned to match observations
-    ϵ::FT = 1               # emissivity of sea ice
+    h::FT                   # ice thickness [m]
+    ρ::FT                   # density of sea ice [kg / m3]
+    c::FT                   # specific heat of sea ice [J / kg / K]
+    T_base::FT              # temperature of sea water at the ice base
+    z0m::FT                 # roughness length for momentum [m]
+    z0b::FT                 # roughness length for tracers [m]
+    T_freeze::FT            # freezing temperature of sea water [K]
+    k_ice::FT               # thermal conductivity of sea ice [W / m / K] (less in HM71)
+    α::FT                   # albedo of sea ice, roughly tuned to match observations
+    ϵ::FT                   # emissivity of sea ice
+    σ::FT                   # Stefan-Boltzmann constant [W / m2 / K4]
+end
+
+"""
+    IceSlabParameters{FT}(coupled_param_dict; h = FT(2), ρ = FT(900), c = FT(2100),
+                          T_base = FT(271.2), z0m = FT(1e-4), z0b = FT(1e-4),
+                          T_freeze = FT(271.2), k_ice = FT(2), α = FT(0.65), ϵ = FT(1))
+
+Initialize the `IceSlabParameters` object with the coupled parameters.
+
+# Arguments
+- `coupled_param_dict`: a dictionary of coupled parameters (required)
+- `h`: ice thickness [m] (default: 2)
+- `ρ`: density of sea ice [kg / m3] (default: 900)
+- `c`: specific heat of sea ice [J / kg / K] (default: 2100)
+- `T_base`: temperature of sea water at the ice base [K] (default: 271.2)
+- `z0m`: roughness length for momentum [m] (default: 1e-4)
+- `z0b`: roughness length for tracers [m] (default: 1e-4)
+- `T_freeze`: freezing temperature of sea water [K] (default: 271.2)
+- `k_ice`: thermal conductivity of sea ice [W / m / K] (default: 2)
+- `α`: albedo of sea ice (default: 0.65)
+- `ϵ`: emissivity of sea ice (default: 1)
+
+# Returns
+- `IceSlabParameters{FT}`: an `IceSlabParameters` object
+"""
+function IceSlabParameters{FT}(
+    coupled_param_dict;
+    h = FT(2),
+    ρ = FT(900),
+    c = FT(2100),
+    T_base = FT(271.2),
+    z0m = FT(1e-4),
+    z0b = FT(1e-4),
+    T_freeze = FT(271.2),
+    k_ice = FT(2),
+    α = FT(0.65),
+    ϵ = FT(1),
+) where {FT}
+    return IceSlabParameters{FT}(;
+        h,
+        ρ,
+        c,
+        T_base,
+        z0m,
+        z0b,
+        T_freeze,
+        k_ice,
+        α,
+        ϵ,
+        σ = coupled_param_dict["stefan_boltzmann_constant"],
+    )
 end
 
 # init simulation
@@ -86,6 +138,7 @@ function PrescribedIceSimulation(
     dt,
     saveat,
     space,
+    coupled_param_dict,
     thermo_params,
     comms_ctx,
     start_date,
@@ -127,7 +180,7 @@ function PrescribedIceSimulation(
     ice_fraction = @. max(min(SIC_init, FT(1) - land_fraction), FT(0))
     ice_fraction = ifelse.(ice_fraction .> FT(0.5), FT(1), FT(0))
 
-    params = IceSlabParameters{FT}()
+    params = IceSlabParameters{FT}(coupled_param_dict)
 
     Y = slab_ice_space_init(FT, space, params)
     cache = (;
@@ -266,7 +319,7 @@ for temperature (curbed at the freezing point).
 """
 function ice_rhs!(dY, Y, p, t)
     FT = eltype(Y)
-    (; k_ice, h, T_base, ρ, c, ϵ, α, T_freeze) = p.params
+    (; k_ice, h, T_base, ρ, c, ϵ, α, T_freeze, σ) = p.params
 
     # Update the cached area fraction with the current SIC
     evaluate!(p.area_fraction, p.SIC_timevaryinginput, t)
@@ -278,9 +331,6 @@ function ice_rhs!(dY, Y, p, t)
 
     # Calculate the conductive flux, and set it to zero if the area fraction is zero
     F_conductive = @. k_ice / (h) * (T_base - Y.T_bulk) # fluxes are defined to be positive when upward
-
-    # TODO: get sigma from parameters
-    σ = FT(5.67e-8)
     rhs = @. (
         -p.F_turb_energy +
         (1 - α) * p.SW_d +