From 2fff9761b76216063e5f4fa5345cd492ce45fd16 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 12 Feb 2025 12:42:45 -0800 Subject: [PATCH 1/3] add ClimaLandSimulation struct --- NEWS.md | 6 + .../components/land/climaland_bucket.jl | 74 +-- .../components/land/climaland_helpers.jl | 74 +++ .../components/land/climaland_integrated.jl | 531 ++++++++++++++++++ .../component_model_tests/climaland_tests.jl | 48 ++ experiments/ClimaEarth/test/runtests.jl | 3 + src/Interfacer.jl | 5 + 7 files changed, 668 insertions(+), 73 deletions(-) create mode 100644 experiments/ClimaEarth/components/land/climaland_helpers.jl create mode 100644 experiments/ClimaEarth/components/land/climaland_integrated.jl create mode 100644 experiments/ClimaEarth/test/component_model_tests/climaland_tests.jl diff --git a/NEWS.md b/NEWS.md index ea5a664e44..b9d8e32b00 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,6 +6,12 @@ ClimaCoupler.jl Release Notes ### ClimaCoupler features +#### Add `ClimaLandSimulation` object PR[#1199](https://github.com/CliMA/ClimaCoupler.jl/pull/1199) +Add methods to support running `ClimaLand.LandModel` in a coupled simulation. +Also add tests to verify the constructor setup and taking a step. +This type is not yet tested within a coupled simulation, but much +of the necessary software infrastructure is added in this PR. + #### Add default `get_field` methods for surface models PR[#1210](https://github.com/CliMA/ClimaCoupler.jl/pull/1210) Add default methods for `get_field` methods that are commonly not extended for surface models. These return reasonable default diff --git a/experiments/ClimaEarth/components/land/climaland_bucket.jl b/experiments/ClimaEarth/components/land/climaland_bucket.jl index c34bf7266e..76150f0cad 100644 --- a/experiments/ClimaEarth/components/land/climaland_bucket.jl +++ b/experiments/ClimaEarth/components/land/climaland_bucket.jl @@ -10,6 +10,7 @@ import ClimaLand.Parameters as LP import ClimaDiagnostics as CD import ClimaCoupler: Checkpointer, FluxCalculator, Interfacer using NCDatasets +include("climaland_helpers.jl") include("../shared/restore.jl") @@ -29,21 +30,6 @@ struct BucketSimulation{M, D, I, A} <: Interfacer.LandModelSimulation end Interfacer.name(::BucketSimulation) = "BucketSimulation" -""" - get_new_cache(p, Y, energy_check) - -Returns a new `p` with an updated field to store e_per_area if energy conservation - checks are turned on. -""" -function get_new_cache(p, Y, energy_check) - if energy_check - e_per_area_field = CC.Fields.zeros(axes(Y.bucket.W)) - return merge(p, (; e_per_area = e_per_area_field)) - else - return p - end -end - """ bucket_init @@ -352,64 +338,6 @@ function Checkpointer.get_model_prog_state(sim::BucketSimulation) return sim.integrator.u.bucket end -""" - temp_anomaly_aquaplanet(coord) - -Introduce a temperature IC anomaly for the aquaplanet case. -The values for this case follow the moist Held-Suarez setup of Thatcher & -Jablonowski (2016, eq. 6), consistent with ClimaAtmos aquaplanet. -""" -temp_anomaly_aquaplanet(coord) = 29 * exp(-coord.lat^2 / (2 * 26^2)) - -""" - temp_anomaly_amip(coord) - -Introduce a temperature IC anomaly for the AMIP case. -The values used in this case have been tuned to align with observed temperature -and result in stable simulations. -""" -temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4 - -""" - make_land_domain( - atmos_boundary_space::CC.Spaces.SpectralElementSpace2D, - zlim::Tuple{FT, FT}, - nelements_vert::Int,) where {FT} - -Creates the land model domain from the horizontal space of the atmosphere, and information -about the number of elements and extent of the vertical domain. -""" -function make_land_domain( - atmos_boundary_space::CC.Spaces.SpectralElementSpace2D, - zlim::Tuple{FT, FT}, - nelements_vert::Int, -) where {FT} - @assert zlim[1] < zlim[2] - depth = zlim[2] - zlim[1] - - mesh = CC.Spaces.topology(atmos_boundary_space).mesh - - radius = mesh.domain.radius - nelements_horz = mesh.ne - npolynomial = CC.Spaces.Quadratures.polynomial_degree(CC.Spaces.quadrature_style(atmos_boundary_space)) - nelements = (nelements_horz, nelements_vert) - vertdomain = CC.Domains.IntervalDomain( - CC.Geometry.ZPoint(FT(zlim[1])), - CC.Geometry.ZPoint(FT(zlim[2])); - boundary_names = (:bottom, :top), - ) - - vertmesh = CC.Meshes.IntervalMesh(vertdomain, CC.Meshes.Uniform(), nelems = nelements[2]) - verttopology = CC.Topologies.IntervalTopology(vertmesh) - vert_center_space = CC.Spaces.CenterFiniteDifferenceSpace(verttopology) - subsurface_space = CC.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space) - space = (; surface = atmos_boundary_space, subsurface = subsurface_space) - - fields = CL.Domains.get_additional_coordinate_field_data(subsurface_space) - - return CL.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space, fields) -end - function Checkpointer.get_model_cache(sim::BucketSimulation) return sim.integrator.p end diff --git a/experiments/ClimaEarth/components/land/climaland_helpers.jl b/experiments/ClimaEarth/components/land/climaland_helpers.jl new file mode 100644 index 0000000000..12e255cf4b --- /dev/null +++ b/experiments/ClimaEarth/components/land/climaland_helpers.jl @@ -0,0 +1,74 @@ +import ClimaCore as CC +import ClimaLand +""" + temp_anomaly_aquaplanet(coord) + +Introduce a temperature IC anomaly for the aquaplanet case. +The values for this case follow the moist Held-Suarez setup of Thatcher & +Jablonowski (2016, eq. 6), consistent with ClimaAtmos aquaplanet. +""" +temp_anomaly_aquaplanet(coord) = 29 * exp(-coord.lat^2 / (2 * 26^2)) + +""" + temp_anomaly_amip(coord) + +Introduce a temperature IC anomaly for the AMIP case. +The values used in this case have been tuned to align with observed temperature +and result in stable simulations. +""" +temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4 + +""" + get_new_cache(p, Y, energy_check) + +Returns a new `p` with an updated field to store e_per_area if energy conservation + checks are turned on. +""" +function get_new_cache(p, Y, energy_check) + if energy_check + e_per_area_field = CC.Fields.zeros(axes(Y.bucket.W)) + return merge(p, (; e_per_area = e_per_area_field)) + else + return p + end +end + +""" + make_land_domain( + atmos_boundary_space::CC.Spaces.SpectralElementSpace2D, + zlim::Tuple{FT, FT}, + nelements_vert::Int,) where {FT} + +Creates the land model domain from the horizontal space of the atmosphere, and information +about the number of elements and extent of the vertical domain. +""" +function make_land_domain( + atmos_boundary_space::CC.Spaces.SpectralElementSpace2D, + zlim::Tuple{FT, FT}, + nelements_vert::Int, +) where {FT} + @assert zlim[1] < zlim[2] + depth = zlim[2] - zlim[1] + + mesh = CC.Spaces.topology(atmos_boundary_space).mesh + + radius = mesh.domain.radius + nelements_horz = mesh.ne + npolynomial = CC.Spaces.Quadratures.polynomial_degree(CC.Spaces.quadrature_style(atmos_boundary_space)) + nelements = (nelements_horz, nelements_vert) + vertdomain = CC.Domains.IntervalDomain( + CC.Geometry.ZPoint(FT(zlim[1])), + CC.Geometry.ZPoint(FT(zlim[2])); + boundary_names = (:bottom, :top), + ) + + vertmesh = CC.Meshes.IntervalMesh(vertdomain, CC.Meshes.Uniform(), nelems = nelements[2]) + verttopology = CC.Topologies.IntervalTopology(vertmesh) + vert_center_space = CC.Spaces.CenterFiniteDifferenceSpace(verttopology) + subsurface_space = CC.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space) + space = (; surface = atmos_boundary_space, subsurface = subsurface_space) + + fields = ClimaLand.Domains.get_additional_coordinate_field_data(subsurface_space) + + return ClimaLand.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space, fields) +end diff --git a/experiments/ClimaEarth/components/land/climaland_integrated.jl b/experiments/ClimaEarth/components/land/climaland_integrated.jl new file mode 100644 index 0000000000..5809258a17 --- /dev/null +++ b/experiments/ClimaEarth/components/land/climaland_integrated.jl @@ -0,0 +1,531 @@ +import ClimaParams +import ClimaLand as CL +import Dates +import ClimaUtilities.TimeVaryingInputs: LinearInterpolation, PeriodicCalendar, TimeVaryingInput +import ClimaCoupler: Checkpointer, FieldExchanger, FluxCalculator, Interfacer, Utilities +import ClimaCore as CC +import SciMLBase +import ClimaTimeSteppers as CTS +import ClimaDiagnostics as CD +import ClimaUtilities.TimeManager: ITime + +include("climaland_helpers.jl") + +""" + ClimaLandSimulation{M, P, D, I, A} + +The integrated ClimaLand model simulation object. +""" +struct ClimaLandSimulation{M, D, I, A} <: Interfacer.LandModelSimulation + model::M + domain::D + integrator::I + area_fraction::A +end + +""" + ClimaLandSimulation( + ::Type{FT}, + dt::TT, + tspan::Tuple{TT, TT}, + start_date::Dates.DateTime, + output_dir::String, + boundary_space, + area_fraction; + saveat::Vector{TT} = [tspan[1], tspan[2]], + domain_type::String = "sphere", + land_temperature_anomaly::String = "amip", + use_land_diagnostics::Bool = true, + ) + +Creates a ClimaLandSimulation object containing a land domain, +a ClimaLand.LandModel, and an integrator. + +This type of model contains a canopy model, soil model, snow model, and +soil CO2 model. Specific details about the complexity of the model +can be found in the ClimaLand.jl documentation. +""" +function ClimaLandSimulation( + ::Type{FT}, + dt::TT, + tspan::Tuple{TT, TT}, + start_date::Dates.DateTime, + output_dir::String, + boundary_space, + area_fraction; + saveat::Vector{TT} = [tspan[1], tspan[2]], + domain_type::String = "sphere", + land_temperature_anomaly::String = "amip", + use_land_diagnostics::Bool = true, +) where {FT, TT <: Union{Float64, ITime}} + @assert domain_type == "sphere" "Currently only spherical shell domains are supported; single column may be supported in the future." + + # Set up domain + depth = FT(10) # soil depth + n_vertical_elements = 10 + domain = make_land_domain(boundary_space, (-depth, FT(0.0)), n_vertical_elements) + surface_space = domain.space.surface + subsurface_space = domain.space.subsurface + + # Set up spatially-varying parameters + earth_param_set = CL.Parameters.LandParameters(FT) + spatially_varying_soil_params = CL.default_spatially_varying_soil_parameters(subsurface_space, surface_space, FT) + # Set up the soil model args + soil_model_type = CL.Soil.EnergyHydrology{FT} + soil_args = create_soil_args(FT, domain, spatially_varying_soil_params) + # Set up the soil microbes model args + soilco2_type = CL.Soil.Biogeochemistry.SoilCO2Model{FT} + soilco2_args = create_soilco2_args(FT, domain) + # Set up the canopy model args + canopy_model_args, canopy_component_args = create_canopy_args(FT, domain, earth_param_set, start_date) + canopy_component_types = (; + autotrophic_respiration = CL.Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = CL.Canopy.TwoStreamModel{FT}, + photosynthesis = CL.Canopy.FarquharModel{FT}, + conductance = CL.Canopy.MedlynConductanceModel{FT}, + hydraulics = CL.Canopy.PlantHydraulicsModel{FT}, + energy = CL.Canopy.BigLeafEnergyModel{FT}, + ) + # setup snow model args + snow_args = create_snow_args(FT, domain, earth_param_set, dt) + snow_model_type = CL.Snow.SnowModel + + # Setup overall land model + f_over = FT(3.28) # 1/m + R_sb = FT(1.484e-4 / 1000) # m/s + runoff_model = + CL.Soil.Runoff.TOPMODELRunoff{FT}(; f_over = f_over, f_max = spatially_varying_soil_params.f_max, R_sb = R_sb) + + Csom = CL.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5)) + + land_input = ( + atmos = CL.CoupledAtmosphere{FT}(), + radiation = CL.CoupledRadiativeFluxes{FT}(), + runoff = runoff_model, + soil_organic_carbon = Csom, + ) + + model = CL.LandModel{FT}(; + soilco2_type = soilco2_type, + soilco2_args = soilco2_args, + land_args = land_input, + soil_model_type = soil_model_type, + soil_args = soil_args, + canopy_component_types = canopy_component_types, + canopy_component_args = canopy_component_args, + canopy_model_args = canopy_model_args, + snow_args = snow_args, + snow_model_type = snow_model_type, + ) + + Y, p, coords = CL.initialize(model) + + # Set initial conditions + (; θ_r, ν) = spatially_varying_soil_params + T_sfc0 = FT(276.85) + @. Y.soil.ϑ_l = θ_r + (ν - θ_r) / 2 + Y.soil.θ_i .= FT(0.0) + ρc_s = CL.Soil.volumetric_heat_capacity.(Y.soil.ϑ_l, Y.soil.θ_i, soil_args.parameters.ρc_ds, earth_param_set) + Y.soil.ρe_int .= CL.Soil.volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_sfc0, earth_param_set) + Y.soilco2.C .= FT(0.000412) + + plant_ν = FT(1.44e-4) + Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν + + # Get temperature anomaly function + T_functions = Dict("aquaplanet" => temp_anomaly_aquaplanet, "amip" => temp_anomaly_amip) + haskey(T_functions, land_temperature_anomaly) || + error("land temp anomaly function $land_temperature_anomaly not supported") + temp_anomaly = T_functions[land_temperature_anomaly] + @. Y.canopy.energy.T = T_sfc0 + temp_anomaly(coords.surface) + Y.snow.S .= 0.0 + Y.snow.S_l .= 0.0 + Y.snow.U .= 0.0 + + set_initial_cache! = CL.make_set_initial_cache(model) + exp_tendency! = CL.make_exp_tendency(model) + imp_tendency! = CL.make_imp_tendency(model) + jacobian! = CL.make_jacobian(model) + set_initial_cache!(p, Y, tspan[1]) + + # set up jacobian info + jac_kwargs = (; jac_prototype = CL.FieldMatrixWithSolver(Y), Wfact = jacobian!) + + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction( + T_exp! = exp_tendency!, + T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), + dss! = CL.dss!, + ), + Y, + tspan, + p, + ) + + # Set up diagnostics + if use_land_diagnostics + netcdf_writer = CD.Writers.NetCDFWriter(subsurface_space, output_dir; start_date) + scheduled_diagnostics = CL.default_diagnostics( + model, + start_date, + output_writer = netcdf_writer, + output_vars = :short, + average_period = :monthly, + ) + diagnostic_handler = CD.DiagnosticsHandler(scheduled_diagnostics, Y, p, tspan[1]; dt = dt) + diag_cb = CD.DiagnosticsCallback(diagnostic_handler) + else + diag_cb = nothing + end + + # Set up time stepper and integrator + stepper = CTS.ARS111() + ode_algo = + CTS.IMEXAlgorithm(stepper, CTS.NewtonsMethod(max_iters = 3, update_j = CTS.UpdateEvery(CTS.NewNewtonIteration))) + integrator = SciMLBase.init( + prob, + ode_algo; + dt = dt, + saveat = saveat, + adaptive = false, + callback = SciMLBase.CallbackSet(diag_cb), + ) + + sim = ClimaLandSimulation(model, (; domain = domain, soil_depth = depth), integrator, area_fraction) + + # DSS state to ensure we have continuous fields + dss_state!(sim) + return sim +end + +############################################################################### +# Helper functions for constructing the land model +############################################################################### +""" + create_canopy_args( + ::Type{FT}, + domain, + earth_param_set, + start_date, + ) + +Creates the arguments for the canopy model. +""" +function create_canopy_args(::Type{FT}, domain, earth_param_set, start_date) where {FT} + surface_space = domain.space.surface + (; Ω, rooting_depth, is_c3, Vcmax25, g1, G_Function, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf) = + CL.clm_canopy_parameters(surface_space) + + # Energy Balance model + ac_canopy = FT(2.5e3) + # Plant Hydraulics and general plant parameters + SAI = FT(0.0) # m2/m2 + f_root_to_shoot = FT(3.5) + RAI = FT(1.0) + K_sat_plant = FT(5e-9) # m/s # seems much too small? + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa + Weibull_param = FT(4) # unitless, Holtzman's original c param value + a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity + conductivity_model = CL.Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = CL.Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a) + plant_ν = FT(1.44e-4) + plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m + n_stem = 0 + n_leaf = 1 + h_stem = FT(0.0) + h_leaf = FT(1.0) + zmax = FT(0.0) + h_canopy = h_stem + h_leaf + compartment_midpoints = n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2] + compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] + + z0_m = FT(0.13) * h_canopy + z0_b = FT(0.1) * z0_m + + # Individual Component arguments + # Set up autotrophic respiration + autotrophic_respiration_args = (; parameters = CL.Canopy.AutotrophicRespirationParameters(FT)) + # Set up radiative transfer + radiative_transfer_args = (; + parameters = CL.Canopy.TwoStreamParameters(FT; Ω, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf, G_Function) + ) + # Set up conductance + conductance_args = (; parameters = CL.Canopy.MedlynConductanceParameters(FT; g1)) + # Set up photosynthesis + photosynthesis_args = (; parameters = CL.Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25)) + # Set up plant hydraulics + modis_lai_artifact_path = CL.Artifacts.modis_lai_forcing_data2008_path() + modis_lai_ncdata_path = joinpath(modis_lai_artifact_path, "Yuan_et_al_2008_1x1.nc") + LAIfunction = CL.prescribed_lai_modis( + modis_lai_ncdata_path, + surface_space, + start_date; + time_interpolation_method = LinearInterpolation(PeriodicCalendar()), + ) + ai_parameterization = CL.Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) + + plant_hydraulics_ps = CL.Canopy.PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + rooting_depth = rooting_depth, + conductivity_model = conductivity_model, + retention_model = retention_model, + ) + plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_midpoints = compartment_midpoints, + compartment_surfaces = compartment_surfaces, + ) + + energy_args = (parameters = CL.Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) + + canopy_component_args = (; + autotrophic_respiration = autotrophic_respiration_args, + radiative_transfer = radiative_transfer_args, + photosynthesis = photosynthesis_args, + conductance = conductance_args, + hydraulics = plant_hydraulics_args, + energy = energy_args, + ) + + shared_params = CL.Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}(z0_m, z0_b, earth_param_set) + + canopy_model_args = (; parameters = shared_params, domain = CL.obtain_surface_domain(domain)) + return canopy_model_args, canopy_component_args +end + +""" + create_soilco2_args(::Type{FT}, domain) + +Creates the arguments for the soil CO2 model. +""" +function create_soilco2_args(::Type{FT}, domain) where {FT} + soilco2_ps = CL.Soil.Biogeochemistry.SoilCO2ModelParameters(FT) + soilco2_top_bc = CL.Soil.Biogeochemistry.AtmosCO2StateBC() + soilco2_bot_bc = CL.Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux + soilco2_sources = (CL.Soil.Biogeochemistry.MicrobeProduction{FT}(),) + soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc) + return (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = domain, + parameters = soilco2_ps, + ) +end + +""" + create_soil_args(::Type{FT}, domain, spatially_varying_soil_params) + +Creates the arguments for the soil model. +""" +function create_soil_args(::Type{FT}, domain, spatially_varying_soil_params) where {FT} + (; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo_dry, + NIR_albedo_dry, + PAR_albedo_wet, + NIR_albedo_wet, + f_max, + ) = spatially_varying_soil_params + + soil_params = CL.Soil.EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo_dry = PAR_albedo_dry, + NIR_albedo_dry = NIR_albedo_dry, + PAR_albedo_wet = PAR_albedo_wet, + NIR_albedo_wet = NIR_albedo_wet, + ) + return (domain = domain, parameters = soil_params) +end + +""" + create_snow_args(::Type{FT}, domain, earth_param_set, dt) + +Creates the arguments for the snow model. +""" +function create_snow_args(::Type{FT}, domain, earth_param_set, dt) where {FT} + snow_parameters = CL.Snow.SnowParameters{FT}(dt; earth_param_set = earth_param_set) + return (; parameters = snow_parameters, domain = CL.obtain_surface_domain(domain)) +end + +############################################################################### +### Functions required by ClimaCoupler.jl for a SurfaceModelSimulation +############################################################################### + +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:area_fraction}) = sim.area_fraction +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:beta}) = + CL.surface_evaporative_scaling(sim.model, sim.integrator.u, sim.integrator.p) +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:emissivity}) = sim.integrator.p.ϵ_sfc +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:energy}) = CL.total_energy(sim.integrator.u, sim.integrator.p) +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:surface_direct_albedo}) = + CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p) +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:surface_diffuse_albedo}) = + CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p) +Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:water}) = CL.total_water(sim.integrator.u, sim.integrator.p) +function Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:height_disp}) + d_canopy = CL.displacement_height(sim.model.canopy, sim.integrator.u, sim.integrator.p) + d_snow = CL.displacement_height(sim.model.snow, sim.integrator.u, sim.integrator.p) + d_soil = CL.displacement_height(sim.model.soil, sim.integrator.u, sim.integrator.p) + return NamedTuple{(:canopy, :snow, :soil)}((d_canopy, d_snow, d_soil)) +end +function Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:roughness_buoyancy}) + z_0b_canopy = sim.model.canopy.parameters.z_0b + z_0b_snow = sim.model.snow.parameters.z_0b + z_0b_soil = sim.model.soil.parameters.z_0b + return NamedTuple{(:canopy, :snow, :soil)}((z_0b_canopy, z_0b_snow, z_0b_soil)) +end +function Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:roughness_momentum}) + z_0m_canopy = sim.model.canopy.parameters.z_0m + z_0m_snow = sim.model.snow.parameters.z_0m + z_0m_soil = sim.model.soil.parameters.z_0m + return NamedTuple{(:canopy, :snow, :soil)}((z_0m_canopy, z_0m_snow, z_0m_soil)) +end +function Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:surface_temperature}) + T_eff = sim.integrator.p.T_sfc + T_canopy = CL.surface_temperature(sim.model.canopy, sim.integrator.u, sim.integrator.p, sim.integrator.t) + T_snow = CL.surface_temperature(sim.model.snow, sim.integrator.u, sim.integrator.p, sim.integrator.t) + T_soil = CL.surface_temperature(sim.model.soil, sim.integrator.u, sim.integrator.p, sim.integrator.t) + return NamedTuple{(:effective, :canopy, :snow, :soil)}((T_eff, T_canopy, T_snow, T_soil)) +end + +# Update fields stored in land model sub-components +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:area_fraction}, field) + parent(sim.area_fraction) .= parent(field) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:diffuse_fraction}, field) + p = sim.integrator.p + parent(p.canopy.radiative_transfer.diffuse_fraction) .= parent(field) +end + +""" + update_field!(sim::ClimaLandSimulation, ::Val{:turbulent_energy_flux}, fields::NamedTuple) + +Takes in a NamedTuple having fields `canopy`, `snow`, and `soil`, which in turn +each contain fields `lhf` and `shf`, and updates each component of the land +simulation with the corresponding LHF and SHF fields. +""" +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:turbulent_energy_flux}, fields::NamedTuple) + turb_energy_flux_canopy = fields.canopy + turb_energy_flux_snow = fields.snow + turb_energy_flux_soil = fields.soil + + p = sim.integrator.p + parent(p.canopy.energy.turbulent_fluxes.lhf) .= parent(turb_energy_flux_canopy.lhf) + parent(p.canopy.energy.turbulent_fluxes.shf) .= parent(turb_energy_flux_canopy.shf) + parent(p.snow.turbulent_fluxes.lhf) .= parent(turb_energy_flux_snow.lhf) + parent(p.snow.turbulent_fluxes.shf) .= parent(turb_energy_flux_snow.shf) + parent(p.soil.turbulent_fluxes.lhf) .= parent(turb_energy_flux_soil.lhf) + parent(p.soil.turbulent_fluxes.shf) .= parent(turb_energy_flux_soil.shf) +end + +""" + update_field!(sim::ClimaLandSimulation, ::Val{:turbulent_moisture_flux}, fields::NamedTuple) + +Takes in a NamedTuple having fields `canopy`, `snow`, and `soil`, and updates +each component of the land simulation with the corresponding field. + +The soil moisture flux is split into liquid and ice components, `Ẽ_l` and `Ẽ_i`. +""" +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:turbulent_moisture_flux}, fields::NamedTuple) + turb_moisture_flux_canopy = fields.canopy + turb_moisture_flux_snow = fields.snow + turb_moisture_flux_soil_liq = fields.soil.Ẽ_l + turb_moisture_flux_soil_ice = fields.soil.Ẽ_i + + p = sim.integrator.p + parent(p.canopy.energy.turbulent_fluxes.transpiration) .= parent(turb_moisture_flux_canopy) + parent(p.snow.vapor_flux) .= parent(turb_moisture_flux_snow) + parent(p.soil.vapor_flux_liq) .= parent(turb_moisture_flux_soil_liq) + parent(p.soil.vapor_flux_ice) .= parent(turb_moisture_flux_soil_ice) +end + +# Update fields stored in land drivers +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:air_temperature}, field) + parent(sim.integrator.p.drivers.T) .= parent(field) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:air_pressure}, field) + parent(sim.integrator.p.drivers.P) .= parent(field) +end +function Interfacer.update_field!(::ClimaLandSimulation, ::Val{:air_humidity}, field) + parent(sim.integrator.p.drivers.q) .= parent(field) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:c_co2}, field) + sim.integrator.p.drivers.c_co2 .= field +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:liquid_precipitation}, field) + ρ_liq = (LP.ρ_cloud_liq(sim.model.soil.parameters.earth_param_set)) + parent(sim.integrator.p.drivers.P_liq) .= parent(field ./ ρ_liq) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:snow_precipitation}, field) + ρ_liq = (LP.ρ_cloud_liq(sim.model.parameters.earth_param_set)) + parent(sim.integrator.p.drivers.P_snow) .= parent(field ./ ρ_liq) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:lw_d}, field) + parent(sim.integrator.p.drivers.LW_d) .= parent(field) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:sw_d}, field) + parent(sim.integrator.p.drivers.SW_d) .= parent(field) +end +function Interfacer.update_field!(sim::ClimaLandSimulation, ::Val{:zenith_angle}, field) + parent(sim.integrator.p.drivers.θs) .= parent(field) +end + +Interfacer.step!(sim::ClimaLandSimulation, t) = Interfacer.step!(sim.integrator, t - sim.integrator.t, true) +Interfacer.reinit!(sim::ClimaLandSimulation, t) = Interfacer.reinit!(sim.integrator, t) + +""" +Extend Interfacer.add_coupler_fields! to add the fields required for ClimaLandSimulation. + +The fields added are: +- `:SW_d` (for radiative transfer) +- `:LW_d` (for radiative transfer) +- `:cos_zenith_angle` (for radiative transfer) +- `:diffuse_fraction` (for radiative transfer) +- `:c_co2` (for photosynthesis, biogeochemistry) +- `:d_sfc` (for turbulent fluxes) +- `:P_air` (for canopy conductance) +- `:T_air` (for canopy conductance) +- `:q_air` (for canopy conductance) +""" +function Interfacer.add_coupler_fields!(coupler_field_names, ::ClimaLandSimulation) + land_coupler_fields = [:SW_d, :LW_d, :cos_zenith_angle, :diffuse_fraction, :c_co2, :d_sfc, :P_air, :T_air, :q_air] + push!(coupler_field_names, land_coupler_fields...) +end + +function Checkpointer.get_model_prog_state(sim::ClimaLandSimulation) + error("get_model_prog_state not implemented") +end + +Interfacer.name(::ClimaLandSimulation) = "ClimaLandSimulation" + + +""" + dss_state!(sim::ClimaLandSimulation)) + +Perform DSS on the state of a component simulation, intended to be used +before the initial step of a run. This method acts on ClimaLand simulations. +The `dss!` function of ClimaLand must be called because it uses either the 2D +or 3D dss buffer stored in the cache depending on space of each variable in +`sim.integrator.u`. +""" +function dss_state!(sim::ClimaLandSimulation) + CL.dss!(sim.integrator.u, sim.integrator.p, sim.integrator.t) +end diff --git a/experiments/ClimaEarth/test/component_model_tests/climaland_tests.jl b/experiments/ClimaEarth/test/component_model_tests/climaland_tests.jl new file mode 100644 index 0000000000..3ff0d599e1 --- /dev/null +++ b/experiments/ClimaEarth/test/component_model_tests/climaland_tests.jl @@ -0,0 +1,48 @@ +using Test +import ClimaCore as CC +import ClimaCoupler +import Dates + +include(joinpath("..", "TestHelper.jl")) +import .TestHelper +include(joinpath("..", "..", "components", "land", "climaland_integrated.jl")) + +FT = Float32 + +@testset "ClimaLandSimulation constructor" begin + dt = Float64(450) + tspan = (Float64(0), 3.0dt) + start_date = Dates.DateTime(2008) + output_dir = pwd() + boundary_space = TestHelper.create_space(FT) + area_fraction = CC.Fields.ones(boundary_space) + + # Construct simulation object + land_sim = ClimaLandSimulation(FT, dt, tspan, start_date, output_dir, boundary_space, area_fraction) + + # Try taking a timestep + Interfacer.step!(land_sim, dt) + + # Check that the simulation object is correctly initialized + @test Interfacer.name(land_sim) == "ClimaLandSimulation" + @test land_sim.area_fraction == area_fraction + + # Check that the state is correctly initialized + state_names = propertynames(land_sim.integrator.u) + @test :canopy in state_names + @test :soil in state_names + @test :snow in state_names + @test :soilco2 in state_names + + # Check that the cache is correctly initialized + cache_names = propertynames(land_sim.integrator.p) + @test :canopy in cache_names + @test :soil in cache_names + @test :snow in cache_names + @test :soilco2 in cache_names + @test :drivers in cache_names + + # Check that the drivers are correctly initialized + driver_names = propertynames(land_sim.integrator.p.drivers) + @test driver_names == (:P_liq, :P_snow, :c_co2, :T, :P, :q, :SW_d, :LW_d, :cosθs, :frac_diff, :soc) +end diff --git a/experiments/ClimaEarth/test/runtests.jl b/experiments/ClimaEarth/test/runtests.jl index ccd36cb1bc..be15787df9 100644 --- a/experiments/ClimaEarth/test/runtests.jl +++ b/experiments/ClimaEarth/test/runtests.jl @@ -7,6 +7,9 @@ gpu_broken = ClimaComms.device() isa ClimaComms.CUDADevice @safetestset "component model test: ClimaAtmos" begin include("component_model_tests/climaatmos_tests.jl") end +@safetestset "component model test: ClimaLand integrated model" begin + include("component_model_tests/climaland_tests.jl") +end @safetestset "component model test: prescr. sea ice" begin include("component_model_tests/prescr_seaice_tests.jl") end diff --git a/src/Interfacer.jl b/src/Interfacer.jl index a5d6c5bd47..ce265ae871 100644 --- a/src/Interfacer.jl +++ b/src/Interfacer.jl @@ -106,6 +106,7 @@ default_coupler_fields() = [ :z0m_sfc, :z0b_sfc, :beta, + :emissivity, :F_turb_energy, :F_turb_moisture, :F_turb_ρτxz, @@ -247,13 +248,17 @@ update_field!( val::Union{ Val{:air_density}, Val{:area_fraction}, + Val{:air_temperature}, Val{:liquid_precipitation}, Val{:radiative_energy_flux_sfc}, Val{:snow_precipitation}, + Val{:thermo_state_int}, Val{:turbulent_energy_flux}, Val{:turbulent_moisture_flux}, Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}, + Val{:surface_pressure}, + Val{:zenith_angle}, }, _, ) = update_field_warning(sim, val) From cd5fd51501d950648023f7712fa4520bb9a21425 Mon Sep 17 00:00:00 2001 From: Julia Sloan Date: Thu, 27 Feb 2025 14:45:08 -0800 Subject: [PATCH 2/3] delete bucket get_new_cache --- .../components/land/climaland_bucket.jl | 10 ++-------- .../components/land/climaland_helpers.jl | 15 --------------- 2 files changed, 2 insertions(+), 23 deletions(-) diff --git a/experiments/ClimaEarth/components/land/climaland_bucket.jl b/experiments/ClimaEarth/components/land/climaland_bucket.jl index 76150f0cad..cc9712ba1c 100644 --- a/experiments/ClimaEarth/components/land/climaland_bucket.jl +++ b/experiments/ClimaEarth/components/land/climaland_bucket.jl @@ -38,7 +38,7 @@ Initializes the bucket model variables. function BucketSimulation( ::Type{FT}, tspan::Tuple{TT, TT}, - config::String, + domain_type::String, albedo_type::String, land_initial_condition::String, land_temperature_anomaly::String, @@ -54,12 +54,7 @@ function BucketSimulation( surface_elevation, use_land_diagnostics::Bool, ) where {FT, TT <: Union{Float64, ITime}} - if config != "sphere" - println( - "Currently only spherical shell domains are supported; single column set-up will be addressed in future PR.", - ) - @assert config == "sphere" - end + @assert domain_type == "sphere" "Currently only spherical shell domains are supported; single column may be supported in the future." α_snow = FT(0.8) # snow albedo if albedo_type == "map_static" # Read in albedo from static data file (default type) @@ -109,7 +104,6 @@ function BucketSimulation( # Initial conditions with no moisture Y, p, coords = CL.initialize(model) - p = get_new_cache(p, Y, energy_check) # Get temperature anomaly function T_functions = Dict("aquaplanet" => temp_anomaly_aquaplanet, "amip" => temp_anomaly_amip) diff --git a/experiments/ClimaEarth/components/land/climaland_helpers.jl b/experiments/ClimaEarth/components/land/climaland_helpers.jl index 12e255cf4b..852c08dabf 100644 --- a/experiments/ClimaEarth/components/land/climaland_helpers.jl +++ b/experiments/ClimaEarth/components/land/climaland_helpers.jl @@ -18,21 +18,6 @@ and result in stable simulations. """ temp_anomaly_amip(coord) = 40 * cosd(coord.lat)^4 -""" - get_new_cache(p, Y, energy_check) - -Returns a new `p` with an updated field to store e_per_area if energy conservation - checks are turned on. -""" -function get_new_cache(p, Y, energy_check) - if energy_check - e_per_area_field = CC.Fields.zeros(axes(Y.bucket.W)) - return merge(p, (; e_per_area = e_per_area_field)) - else - return p - end -end - """ make_land_domain( atmos_boundary_space::CC.Spaces.SpectralElementSpace2D, From 780d984c71607b0b22b827de2068a8e89b6df77d Mon Sep 17 00:00:00 2001 From: Julia Sloan Date: Thu, 27 Feb 2025 16:07:08 -0800 Subject: [PATCH 3/3] use CL branch tr/support-coupling-LandModel --- experiments/ClimaEarth/Manifest-v1.11.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/experiments/ClimaEarth/Manifest-v1.11.toml b/experiments/ClimaEarth/Manifest-v1.11.toml index 407fb2fecd..a6bd517233 100644 --- a/experiments/ClimaEarth/Manifest-v1.11.toml +++ b/experiments/ClimaEarth/Manifest-v1.11.toml @@ -349,7 +349,9 @@ version = "0.2.12" [[deps.ClimaLand]] deps = ["ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"] -git-tree-sha1 = "84d93ae74deda75c29a69f0660556e82cb5e45de" +git-tree-sha1 = "9657340ac78e5eee92f85f5515c3721658b45e55" +repo-rev = "tr/support-coupling-LandModel" +repo-url = "https://github.com/CliMA/ClimaLand.jl.git" uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532" version = "0.15.10"