use nearest neighbor interpolation (#1290)

Author: kmdeck

.buildkite/longruns_gpu/pipeline.yml

@@ -53,6 +53,18 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
+      - label: ":snow_capped_mountain: Low-Res Snowy Land"
+        command:
+          - julia --color=yes --project=.buildkite experiments/long_runs/low_res_snowy_land.jl
+        artifact_paths:
+          - "lowres_snowy_land_longrun_gpu/*png"
+          - "lowres_snowy_land_longrun_gpu/*pdf"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 3:00:00
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+
       - label: ":sunglasses: California regional simulation"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/land_region.jl

NEWS.md

@@ -3,6 +3,7 @@ ClimaLand.jl Release Notes
 
 main
 -------
+- Use nearest neighbor spatial interpolation PR[#1290](https://github.com/CliMA/ClimaLand.jl/pull/1290)
 - Use global MODIS LAI for all simulations; remove `modis_lai_fluxnet_sites` artifact PR[#1282](https://github.com/CliMA/ClimaLand.jl/pull/1282)
 - Use ClimaUtilities v0.1.25 to read in spatial data to a point using lat/lon PR[#1279](https://github.com/CliMA/ClimaLand.jl/pull/1279)
 - Add data handling tools to FluxnetSimulationsExt and use throughout docs and experiments PR[#1238](https://github.com/CliMA/ClimaLand.jl/pull/1238)

experiments/long_runs/land_region.jl

@@ -51,8 +51,6 @@ outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
 
 function setup_model(FT, context, start_date, Δt, domain, earth_param_set)
-
-    time_interpolation_method = LinearInterpolation(PeriodicCalendar())
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
 
@@ -66,7 +64,7 @@ function setup_model(FT, context, start_date, Δt, domain, earth_param_set)
         earth_param_set,
         FT;
         max_wind_speed = 25.0,
-        time_interpolation_method,
+        time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
     )
 
     (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) =
@@ -179,14 +177,14 @@ function setup_model(FT, context, start_date, Δt, domain, earth_param_set)
     # Set up plant hydraulics
     modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(;
         context = context,
-        start_date = start_date + Second(t0),
-        end_date = start_date + Second(t0) + Second(tf),
+        start_date,
+        end_date = stop_date,
     )
     LAIfunction = ClimaLand.prescribed_lai_modis(
         modis_lai_ncdata_path,
         surface_space,
         start_date;
-        time_interpolation_method,
+        time_interpolation_method = LinearInterpolation(),
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)

experiments/long_runs/low_res_snowy_land.jl

@@ -0,0 +1,322 @@
+# # Global run of land model at low resolution
+
+# The code sets up and runs ClimaLand v1, which
+# includes soil, canopy, and snow, on a spherical domain,
+# using low resolution ERA5 data as forcing.
+
+# Simulation Setup
+# Number of spatial elements: 30 in horizontal, 15 in vertical
+# Soil depth: 50 m
+# Simulation duration: 730 d
+# Timestep: 450 s
+# Timestepper: ARS111
+# Fixed number of iterations: 3
+# Jacobian update: every new Newton iteration
+# Atmos forcing update: every 3 hours
+
+import ClimaComms
+ClimaComms.@import_required_backends
+using ClimaUtilities.ClimaArtifacts
+import ClimaUtilities.TimeManager: ITime, date
+
+import ClimaDiagnostics
+import ClimaUtilities
+
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
+import ClimaUtilities.ClimaArtifacts: @clima_artifact
+import ClimaParams as CP
+using ClimaCore
+using ClimaLand
+using ClimaLand.Snow
+using ClimaLand.Soil
+using ClimaLand.Canopy
+import ClimaLand
+import ClimaLand.Parameters as LP
+import ClimaLand.Simulations: LandSimulation, solve!
+
+using Dates
+
+using CairoMakie, GeoMakie, Poppler_jll, ClimaAnalysis
+LandSimVis =
+    Base.get_extension(
+        ClimaLand,
+        :LandSimulationVisualizationExt,
+    ).LandSimulationVisualizationExt;
+
+const FT = Float64;
+context = ClimaComms.context()
+ClimaComms.init(context)
+device = ClimaComms.device()
+device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
+root_path = "lowres_snowy_land_longrun_$(device_suffix)"
+diagnostics_outdir = joinpath(root_path, "global_diagnostics")
+outdir =
+    ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
+
+function setup_model(FT, start_date, stop_date, Δt, domain, earth_param_set)
+    surface_space = domain.space.surface
+    subsurface_space = domain.space.subsurface
+    # Forcing data
+    era5_ncdata_path = ClimaLand.Artifacts.era5_land_forcing_data2008_path(;
+        context,
+        lowres = true,
+    )
+    atmos, radiation = ClimaLand.prescribed_forcing_era5(
+        era5_ncdata_path,
+        surface_space,
+        start_date,
+        earth_param_set,
+        FT;
+        max_wind_speed = 25.0,
+        time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+    )
+
+    (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) =
+        ClimaLand.Soil.soil_composition_parameters(subsurface_space, FT)
+    (; ν, hydrology_cm, K_sat, θ_r) =
+        ClimaLand.Soil.soil_vangenuchten_parameters(subsurface_space, FT)
+    soil_albedo = Soil.CLMTwoBandSoilAlbedo{FT}(;
+        ClimaLand.Soil.clm_soil_albedo_parameters(surface_space)...,
+    )
+    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
+    soil_params = Soil.EnergyHydrologyParameters(
+        FT;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        albedo = soil_albedo,
+    )
+    f_over = FT(3.28) # 1/m
+    R_sb = FT(1.484e-4 / 1000) # m/s
+    runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
+        f_over = f_over,
+        f_max = ClimaLand.Soil.topmodel_fmax(surface_space, FT),
+        R_sb = R_sb,
+    )
+
+    # Spatially varying canopy parameters from CLM
+    g1 = ClimaLand.Canopy.clm_medlyn_g1(surface_space)
+    rooting_depth = ClimaLand.Canopy.clm_rooting_depth(surface_space)
+    (; is_c3, Vcmax25) =
+        ClimaLand.Canopy.clm_photosynthesis_parameters(surface_space)
+    (; Ω, G_Function, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf) =
+        ClimaLand.Canopy.clm_canopy_radiation_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(7e-8) # m/s
+    ψ63 = FT(-4 / 0.0098) # / MPa to m
+    Weibull_param = FT(4) # unitless
+    a = FT(0.2 * 0.0098) # 1/m
+    conductivity_model =
+        Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
+    retention_model = 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
+
+    soil_args = (domain = domain, parameters = soil_params)
+    soil_model_type = Soil.EnergyHydrology{FT}
+
+    # Soil microbes model
+    soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
+    soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
+    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
+
+    soilco2_args = (; domain = domain, parameters = soilco2_ps)
+
+    # Now we set up the canopy model, which we set up by component:
+    # Component Types
+    canopy_component_types = (;
+        autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT},
+        radiative_transfer = Canopy.TwoStreamModel{FT},
+        photosynthesis = Canopy.FarquharModel{FT},
+        conductance = Canopy.MedlynConductanceModel{FT},
+        hydraulics = Canopy.PlantHydraulicsModel{FT},
+        energy = Canopy.BigLeafEnergyModel{FT},
+    )
+    # Individual Component arguments
+    # Set up autotrophic respiration
+    autotrophic_respiration_args =
+        (; parameters = Canopy.AutotrophicRespirationParameters(FT))
+    # Set up radiative transfer
+    radiative_transfer_args = (;
+        parameters = Canopy.TwoStreamParameters(
+            FT;
+            Ω,
+            α_PAR_leaf,
+            τ_PAR_leaf,
+            α_NIR_leaf,
+            τ_NIR_leaf,
+            G_Function,
+        )
+    )
+    # Set up conductance
+    conductance_args =
+        (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
+    # Set up photosynthesis
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
+    # Set up plant hydraulics
+    modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(;
+        context = nothing,
+        start_date,
+        end_date = stop_date,
+    )
+    LAIfunction = ClimaLand.prescribed_lai_modis(
+        modis_lai_ncdata_path,
+        surface_space,
+        start_date;
+        time_interpolation_method = LinearInterpolation(),
+    )
+    ai_parameterization =
+        Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
+
+    plant_hydraulics_ps = 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 = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),)
+
+    # Canopy component args
+    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,
+    )
+
+    # Other info needed
+    shared_params = Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}(
+        z0_m,
+        z0_b,
+        earth_param_set,
+    )
+
+    canopy_model_args = (;
+        parameters = shared_params,
+        domain = ClimaLand.obtain_surface_domain(domain),
+    )
+
+    # Snow model
+    #    α_snow = Snow.ConstantAlbedoModel(FT(0.7))
+    # Set β = 0 in order to regain model without density dependence
+    α_snow = Snow.ZenithAngleAlbedoModel(
+        FT(0.64),
+        FT(0.06),
+        FT(2);
+        β = FT(0.4),
+        x0 = FT(0.2),
+    )
+    horz_degree_res =
+        sum(ClimaLand.Domains.average_horizontal_resolution_degrees(domain)) / 2 # mean of resolution in latitude and longitude, in degrees
+    scf = Snow.WuWuSnowCoverFractionModel(
+        FT(0.08),
+        FT(1.77),
+        FT(1.0),
+        horz_degree_res,
+    )
+    snow_parameters = SnowParameters{FT}(
+        Δt;
+        earth_param_set = earth_param_set,
+        α_snow = α_snow,
+        scf = scf,
+    )
+    snow_args = (;
+        parameters = snow_parameters,
+        domain = ClimaLand.obtain_surface_domain(domain),
+    )
+    snow_model_type = Snow.SnowModel
+
+    land_input = (
+        atmos = atmos,
+        radiation = radiation,
+        runoff = runoff_model,
+        soil_organic_carbon = Csom,
+    )
+    land = 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,
+    )
+    return land
+end
+# Note that since the Northern hemisphere's winter season is defined as DJF,
+# we simulate from and until the beginning of
+# March so that a full season is included in seasonal metrics.
+start_date = DateTime("2008-03-01")
+stop_date = DateTime("2010-03-02")
+Δt = 450.0
+nelements = (30, 15)
+domain = ClimaLand.Domains.global_domain(
+    FT;
+    context,
+    nelements,
+    mask_threshold = FT(0.99),
+)
+params = LP.LandParameters(FT)
+model = setup_model(FT, start_date, stop_date, Δt, domain, params)
+user_callbacks = (
+    ClimaLand.NaNCheckCallback(
+        Dates.Month(6),
+        start_date,
+        ITime(Δt, epoch = start_date),
+        mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)),
+    ),
+    ClimaLand.ReportCallback(10000),
+)
+simulation =
+    LandSimulation(FT, start_date, stop_date, Δt, model; user_callbacks, outdir)
+@info "Run: Global Soil-Canopy-Snow Model"
+@info "Resolution: $nelements"
+@info "Timestep: $Δt s"
+@info "Start Date: $start_date"
+@info "Stop Date: $stop_date"
+ClimaLand.Simulations.solve!(simulation)
+
+LandSimVis.make_annual_timeseries(simulation; savedir = root_path)
+LandSimVis.make_heatmaps(simulation; savedir = root_path, date = stop_date)
+LandSimVis.make_leaderboard_plots(simulation; savedir = root_path)