diff --git a/experiments/long_runs/land_leaderboard.jl b/experiments/long_runs/land_leaderboard.jl
new file mode 100644
index 0000000000..4c99d109af
--- /dev/null
+++ b/experiments/long_runs/land_leaderboard.jl
@@ -0,0 +1,596 @@
+# # Global run of land model
+
+# The code sets up and runs the soil/canopy model for 6 hours on a spherical domain,
+# using ERA5 data. In this simulation, we have
+# turned lateral flow off because horizontal boundary conditions and the
+# land/sea mask are not yet supported by ClimaCore.
+
+# Simulation Setup
+# Number of spatial elements: 101 in horizontal, 15 in vertical
+# Soil depth: 50 m
+# Simulation duration: 2.1 * 365 d
+# Timestep: 450 s
+# Timestepper: ARS111
+# Fixed number of iterations: 3
+# Jacobian update: every new Newton iteration
+# Atmos forcing update: every 3 hours
+import SciMLBase
+import ClimaComms
+ClimaComms.@import_required_backends
+import ClimaTimeSteppers as CTS
+using ClimaCore
+using ClimaUtilities.ClimaArtifacts
+import Interpolations
+
+import ClimaDiagnostics
+import ClimaUtilities
+
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
+import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
+import ClimaUtilities.Regridders: InterpolationsRegridder
+import ClimaUtilities.ClimaArtifacts: @clima_artifact
+import ClimaParams as CP
+
+using ClimaLand
+using ClimaLand.Soil
+using ClimaLand.Canopy
+import ClimaLand
+import ClimaLand.Parameters as LP
+
+using Statistics
+using Dates
+import NCDatasets
+
+const FT = Float64;
+time_interpolation_method = LinearInterpolation(PeriodicCalendar())
+regridder_type = :InterpolationsRegridder
+context = ClimaComms.context()
+device = ClimaComms.device()
+device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu"
+root_path = "land_longrun_$(device_suffix)"
+diagnostics_outdir = joinpath(root_path, "global_diagnostics")
+outdir =
+    ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
+
+function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
+
+    earth_param_set = LP.LandParameters(FT)
+    radius = FT(6378.1e3)
+    depth = FT(50)
+    domain = ClimaLand.Domains.SphericalShell(;
+        radius = radius,
+        depth = depth,
+        nelements = nelements,
+        npolynomial = 1,
+        dz_tuple = FT.((10.0, 0.05)),
+    )
+    surface_space = domain.space.surface
+    subsurface_space = domain.space.subsurface
+
+    start_date = DateTime(2012)
+    # Forcing data
+    era5_artifact_path =
+        ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context)
+    era5_ncdata_path = joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc")
+    atmos, radiation = ClimaLand.prescribed_forcing_era5(
+        era5_ncdata_path,
+        surface_space,
+        start_date,
+        earth_param_set,
+        FT;
+        time_interpolation_method = time_interpolation_method,
+        regridder_type = regridder_type,
+    )
+
+    soil_params_artifact_path =
+        ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
+    extrapolation_bc = (
+        Interpolations.Periodic(),
+        Interpolations.Flat(),
+        Interpolations.Flat(),
+    )
+    soil_params_mask = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "α",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
+    )
+    oceans_to_value(field, mask, value) =
+        mask == 1.0 ? field : eltype(field)(value)
+
+    vg_α = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "α",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    vg_n = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "vGn_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "n",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    x = parent(vg_α)
+    μ = mean(log10.(x[x .> 0]))
+    vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ)
+
+    x = parent(vg_n)
+    μ = mean(x[x .> 0])
+    vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ)
+
+    vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
+        ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...)
+    hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n)
+
+    θ_r = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "residual_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "θ_r",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+
+    ν = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "porosity_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "ν",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    K_sat = SpaceVaryingInput(
+        joinpath(
+            soil_params_artifact_path,
+            "ksat_map_gupta_etal2020_1.0x1.0x4.nc",
+        ),
+        "Ksat",
+        subsurface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+
+    x = parent(K_sat)
+    μ = mean(log10.(x[x .> 0]))
+    K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ)
+
+    ν .= oceans_to_value.(ν, soil_params_mask, 1)
+
+    θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0)
+
+
+    S_s =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
+            soil_params_mask,
+            1,
+        )
+    ν_ss_om =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    ν_ss_quartz =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    ν_ss_gravel =
+        oceans_to_value.(
+            ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1),
+            soil_params_mask,
+            0,
+        )
+    PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2)
+    soil_params = Soil.EnergyHydrologyParameters(
+        FT;
+        ν,
+        ν_ss_om,
+        ν_ss_quartz,
+        ν_ss_gravel,
+        hydrology_cm,
+        K_sat,
+        S_s,
+        θ_r,
+        PAR_albedo = PAR_albedo,
+        NIR_albedo = NIR_albedo,
+    )
+
+    soil_params_mask_sfc =
+        ClimaLand.Domains.top_center_to_surface(soil_params_mask)
+
+    # Read in f_max data and land sea mask
+    infile_path = ClimaLand.Artifacts.topmodel_data_path()
+    f_max =
+        SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type)
+    mask = SpaceVaryingInput(
+        infile_path,
+        "landsea_mask",
+        surface_space;
+        regridder_type,
+    )
+    # Unsure how to handle two masks
+    f_max = oceans_to_value.(f_max, mask, FT(0.0))
+    f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
+    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 = f_max,
+        R_sb = R_sb,
+    )
+
+    # Energy Balance model
+    ac_canopy = FT(2.5e3)
+
+    #clm_data is used for g1, vcmax, rooting, and two_stream param maps
+    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
+    # Foliage clumping index data derived from MODIS
+    modis_ci_artifact_path =
+        ClimaLand.Artifacts.modis_ci_data_folder_path(; context)
+
+    # TwoStreamModel parameters
+    nans_to_one(x) = isnan(x) ? eltype(x)(1) : x
+    Ω = SpaceVaryingInput(
+        joinpath(modis_ci_artifact_path, "He_et_al_2012_1x1.nc"),
+        "ci",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = nans_to_one,),
+    )
+    χl = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "xl",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    G_Function = CLMGFunction(χl)
+    α_PAR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rholvis",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    τ_PAR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "taulvis",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    α_NIR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rholnir",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    τ_NIR_leaf = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "taulnir",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    # Conductance Model
+    # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa)
+    g1 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "medlynslope",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),),
+    )
+
+    #Photosynthesis model
+    # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s
+    Vcmax25 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "vcmx25",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
+    )
+    # photosynthesis mechanism is read as a float, where 1.0 indicates c3 and 0.0 c4
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    # 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 =
+        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
+    rooting_depth = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rooting_depth",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    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
+
+
+    soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
+
+    soil_args = (domain = domain, parameters = soil_params)
+    soil_model_type = Soil.EnergyHydrology{FT}
+
+    # Soil microbes model
+    soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
+
+    # soil microbes args
+    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
+
+    # Set the soil CO2 BC to being atmospheric CO2
+    soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
+    soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux
+    soilco2_sources = (Soil.Biogeochemistry.MicrobeProduction{FT}(),)
+
+    soilco2_boundary_conditions =
+        (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
+
+    soilco2_args = (;
+        boundary_conditions = soilco2_boundary_conditions,
+        sources = soilco2_sources,
+        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
+
+    LAIfunction = ClimaLand.prescribed_lai_era5(
+        era5_ncdata_path,
+        surface_space,
+        start_date;
+        time_interpolation_method = time_interpolation_method,
+        regridder_type = regridder_type,
+    )
+    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),
+    )
+
+    # Integrated plant hydraulics and soil model
+    land_input = (
+        atmos = atmos,
+        radiation = radiation,
+        runoff = runoff_model,
+        soil_organic_carbon = Csom,
+    )
+    land = SoilCanopyModel{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,
+    )
+
+    Y, p, cds = initialize(land)
+
+    init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2
+    Y.soil.ϑ_l .= init_soil.(ν, θ_r)
+    Y.soil.θ_i .= FT(0.0)
+    T = FT(276.85)
+    ρc_s =
+        Soil.volumetric_heat_capacity.(
+            Y.soil.ϑ_l,
+            Y.soil.θ_i,
+            soil_params.ρc_ds,
+            soil_params.earth_param_set,
+        )
+    Y.soil.ρe_int .=
+        Soil.volumetric_internal_energy.(
+            Y.soil.θ_i,
+            ρc_s,
+            T,
+            soil_params.earth_param_set,
+        )
+    Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
+    Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν
+    evaluate!(Y.canopy.energy.T, atmos.T, t0)
+
+    set_initial_cache! = make_set_initial_cache(land)
+    exp_tendency! = make_exp_tendency(land)
+    imp_tendency! = ClimaLand.make_imp_tendency(land)
+    jacobian! = ClimaLand.make_jacobian(land)
+    set_initial_cache!(p, Y, t0)
+
+    # set up jacobian info
+    jac_kwargs =
+        (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!)
+
+    prob = SciMLBase.ODEProblem(
+        CTS.ClimaODEFunction(
+            T_exp! = exp_tendency!,
+            T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
+            dss! = ClimaLand.dss!,
+        ),
+        Y,
+        (t0, tf),
+        p,
+    )
+
+    updateat = Array(t0:(3600 * 3):tf)
+    drivers = ClimaLand.get_drivers(land)
+    updatefunc = ClimaLand.make_update_drivers(drivers)
+
+    # ClimaDiagnostics
+
+    nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, outdir)
+
+    diags = ClimaLand.default_diagnostics(
+        land,
+        start_date;
+        output_writer = nc_writer,
+        output_vars = :short,
+        average_period = :monthly,
+    )
+
+    diagnostic_handler =
+        ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt)
+
+    diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler)
+
+    driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+    return prob, SciMLBase.CallbackSet(driver_cb, diag_cb)
+end
+
+function setup_and_solve_problem(; greet = false)
+
+    t0 = 0.0
+    tf = 2.1 * 60 * 60.0 * 24 * 365
+    Δt = 450.0
+    nelements = (101, 15)
+    if greet
+        @info "Run: Global Soil-Canopy Model"
+        @info "Resolution: $nelements"
+        @info "Timestep: $Δt s"
+        @info "Duration: $(tf - t0) s"
+    end
+
+    prob, cb = setup_prob(t0, tf, Δt; nelements)
+
+    # Define timestepper and ODE algorithm
+    stepper = CTS.ARS111()
+    ode_algo = CTS.IMEXAlgorithm(
+        stepper,
+        CTS.NewtonsMethod(
+            max_iters = 3,
+            update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+        ),
+    )
+    SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb, adaptive = false)
+    return nothing
+end
+
+setup_and_solve_problem(; greet = true);
+
+# read in diagnostics and make some plots!
+#### ClimaAnalysis ####
+include("leaderboard/leaderboard.jl")
+diagnostics_folder_path = outdir
+leaderboard_base_path = root_path
+compute_leaderboard(leaderboard_base_path, diagnostics_folder_path)
diff --git a/experiments/long_runs/leaderboard/data_sources.jl b/experiments/long_runs/leaderboard/data_sources.jl
new file mode 100644
index 0000000000..55575975b4
--- /dev/null
+++ b/experiments/long_runs/leaderboard/data_sources.jl
@@ -0,0 +1,130 @@
+import ClimaAnalysis
+
+# To add a variable for the leaderboard, add a key-value pair to the dictionary
+# `sim_var_dict` whose key is the short name of the variable and the value is a function
+# that returns a `OutputVar` Any preprocessing is done in the function which includes unit
+# conversion and shifting the dates.
+#
+# Then, add a key-value pair to the dictionary `obs_var_dict` whose key is the same short
+# name as before and the value is a function that takes in a start date and returns a
+# `OutputVar`. Any preprocessing is done in the function.
+#
+# Next, add a key-value pair to the dictionary `mask_dict` whose key is the same short name
+# as before and the value is a function that takes in a OutputVar representing simulation
+# data and a OutputVar representing observational data and returns a masking function or
+# `nothing` if no masking function is needed. The masking function is used to correctly
+# normalize the global bias and global RMSE.
+#
+# Finally, add a key-value pair to the dictionary `compare_vars_biases_plot_extrema` whose
+# key is the same short name as before and the value is a tuple of floats which determine
+# the range of the bias plots.
+
+# For plotting bias plots
+compare_vars_biases_plot_extrema =
+    Dict("et" => (-0.00001, 0.00001), "gpp" => (-8.0, 8.0), "lwu" => (-40.0, 40.0))
+
+# Dict for loading in simulation data
+sim_var_dict = Dict{String,Any}()
+
+for short_name in ["et", "lwu"]
+    sim_var_dict[short_name] =
+        () -> begin
+            sim_var = get(
+                ClimaAnalysis.SimDir(diagnostics_folder_path),
+                short_name = short_name,
+            )
+            # Remove the line below later when start_date is added to the diagnostics (see issue #879 in ClimaLand.jl)
+            haskey(sim_var.attributes, "start_date") || (sim_var.attributes["start_date"] = "2012-01-01T00:00:00")
+            sim_var = ClimaAnalysis.shift_to_start_of_previous_month(sim_var)
+            return sim_var
+        end
+end
+
+for short_name in ["gpp"]
+    sim_var_dict[short_name] =
+        () -> begin
+            sim_var = get(
+                ClimaAnalysis.SimDir(diagnostics_folder_path),
+                short_name = short_name,
+            )
+            # Remove the line below later when start_date is added to the diagnostics (see issue #879 in ClimaLand.jl)
+            haskey(sim_var.attributes, "start_date") || (sim_var.attributes["start_date"] = "2012-01-01T00:00:00")
+            sim_var = ClimaAnalysis.shift_to_start_of_previous_month(sim_var)
+            sim_var = ClimaAnalysis.convert_units(
+                sim_var,
+                "g m-2 day-1",
+                conversion_function = units -> units * 86400.0 * 12.011,
+            )
+            return sim_var
+        end
+end
+
+# Dict for loading in observational data
+obs_var_dict = Dict{String,Any}()
+obs_var_dict["et"] =
+    (start_date) -> begin
+        obs_var = ClimaAnalysis.OutputVar(
+            ClimaLand.Artifacts.ilamb_dataset_path(; context = "evspsbl_MODIS_et_0.5x0.5.nc"),
+            "et",
+            new_start_date = start_date,
+            shift_by = Dates.firstdayofmonth,
+        )
+        ClimaAnalysis.units(obs_var) == "kg/m2/s" &&
+            (obs_var = ClimaAnalysis.set_units(obs_var, "kg m^-2 s^-1"))
+        obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
+        return obs_var
+    end
+
+obs_var_dict["gpp"] =
+    (start_date) -> begin
+        obs_var = ClimaAnalysis.OutputVar(
+            ClimaLand.Artifacts.ilamb_dataset_path(; context = "gpp_FLUXCOM_gpp.nc"),
+            "gpp",
+            new_start_date = start_date,
+            shift_by = Dates.firstdayofmonth,
+        )
+        ClimaAnalysis.dim_units(obs_var, "lon") == "degree" &&
+            (obs_var.dim_attributes["lon"]["units"] = "degrees_east")
+        ClimaAnalysis.dim_units(obs_var, "lat") == "degree" &&
+            (obs_var.dim_attributes["lat"]["units"] = "degrees_north")
+        obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
+        return obs_var
+    end
+
+obs_var_dict["lwu"] =
+    (start_date) -> begin
+        obs_var = ClimaAnalysis.OutputVar(
+            ClimaLand.Artifacts.ilamb_dataset_path(; context = "rlus_CERESed4.2_rlus.nc"),
+            "rlus",
+            new_start_date = start_date,
+            shift_by = Dates.firstdayofmonth,
+        )
+        ClimaAnalysis.units(obs_var) == "W m-2" &&
+            (obs_var = ClimaAnalysis.set_units(obs_var, "W m^-2"))
+        return obs_var
+    end
+
+# Dict for loading in masks
+mask_dict = Dict{String,Any}()
+
+mask_dict["et"] =
+    (sim_var, obs_var) -> begin
+        return ClimaAnalysis.make_lonlat_mask(
+            ClimaAnalysis.slice(obs_var, time = ClimaAnalysis.times(obs_var) |> first);
+            set_to_val = isnan,
+            true_val = 0.0
+        )
+    end
+
+mask_dict["gpp"] =
+    (sim_var, obs_var) -> begin
+        return ClimaAnalysis.make_lonlat_mask(
+            ClimaAnalysis.slice(obs_var, time = ClimaAnalysis.times(obs_var) |> first);
+            set_to_val = isnan,
+            true_val = 0.0
+        )
+    end
+
+mask_dict["lwu"] = (sim_var, obs_var) -> begin
+    return nothing
+end
diff --git a/experiments/long_runs/leaderboard/leaderboard.jl b/experiments/long_runs/leaderboard/leaderboard.jl
new file mode 100644
index 0000000000..8ae72e16c1
--- /dev/null
+++ b/experiments/long_runs/leaderboard/leaderboard.jl
@@ -0,0 +1,137 @@
+import ClimaAnalysis
+import GeoMakie
+import CairoMakie
+import Dates
+
+include("data_sources.jl")
+
+"""
+    compute_leaderboard(leaderboard_base_path, diagnostics_folder_path)
+
+Plot the biases against observations and simulation data.
+
+The paramter `leaderboard_base_path` is the path to save the leaderboards and bias plots and
+`diagnostics_folder_path` is the path to the simulation data.
+"""
+function compute_leaderboard(leaderboard_base_path, diagnostics_folder_path)
+    @info "Error against observations"
+
+    # Set up dict for storing simulation and observational data after processing
+    sim_obs_comparsion_dict = Dict()
+
+    # Print dates for debugging
+    _, var_func = first(sim_var_dict)
+    var = var_func()
+    output_dates =
+        Dates.DateTime(var.attributes["start_date"]) .+
+        Dates.Second.(ClimaAnalysis.times(var))
+    @info "Working with dates:"
+    @info output_dates
+
+    short_names = keys(sim_var_dict)
+    for short_name in short_names
+        @info short_name
+        # Simulation data
+        sim_var = sim_var_dict[short_name]()
+
+        # Observational data
+        obs_var = obs_var_dict[short_name](sim_var.attributes["start_date"])
+
+        # Remove first spin_up_months months from simulation
+        spin_up_months = 12
+        spinup_cutoff = spin_up_months * 30 * 86400.0
+        ClimaAnalysis.times(sim_var)[end] >= spinup_cutoff &&
+            (sim_var = ClimaAnalysis.window(sim_var, "time", left = spinup_cutoff))
+
+        # Get 12 or less months of data
+        num_times = length(ClimaAnalysis.times(sim_var))
+        num_times > 12 && (sim_var = ClimaAnalysis.window(sim_var, "time", right = ClimaAnalysis.times(sim_var)[12]))
+
+        obs_var = ClimaAnalysis.resampled_as(obs_var, sim_var)
+
+        sim_obs_comparsion_dict[short_name] = (sim_var, obs_var)
+    end
+
+    # Get the right number of months for plotting
+    _, sim_obs_tuple = first(sim_obs_comparsion_dict)
+    num_times = length(ClimaAnalysis.times(sim_obs_tuple[1]))
+    months = Dates.monthname.(1:num_times)
+
+    # Plot monthly comparsions
+    for short_name in short_names
+        fig = CairoMakie.Figure(size = (650 * ceil(num_times / 2), 450 * 2))
+        sim_var, obs_var = sim_obs_comparsion_dict[short_name]
+        mask = mask_dict[short_name](sim_var, obs_var)
+        times = ClimaAnalysis.times(sim_var) |> copy
+        times = vcat(times, Array{Union{Missing, eltype(times)}}(missing, 12 - num_times))
+        times = reshape(times, (2, 6))
+        for ((indices, t), month) in zip(pairs(times), months)
+            layout = fig[Tuple(indices)...] = CairoMakie.GridLayout()
+            ClimaAnalysis.Visualize.plot_bias_on_globe!(
+                layout,
+                ClimaAnalysis.slice(sim_var, time = t),
+                ClimaAnalysis.slice(obs_var, time = t),
+                cmap_extrema = compare_vars_biases_plot_extrema[short_name],
+                mask = mask,
+            )
+            CairoMakie.Label(layout[0, 1], month, tellwidth = false, fontsize = 30)
+        end
+        CairoMakie.save(joinpath(leaderboard_base_path, "$(short_name)_bias_plot.png"), fig)
+    end
+
+    # Plot month (x-axis) and global bias and global RMSE (y-axis)
+    fig = CairoMakie.Figure(size = (450 * length(keys(sim_var_dict)), 600))
+    for (col, short_name) in enumerate(keys(sim_var_dict))
+        sim_var, obs_var = sim_obs_comparsion_dict[short_name]
+        mask = mask_dict[short_name](sim_var, obs_var)
+        # Compute globas bias and global rmse
+        rmse_vec = []
+        bias_vec = []
+        times = ClimaAnalysis.times(sim_var) |> copy
+        for t in times
+            g_rmse = ClimaAnalysis.global_rmse(
+                ClimaAnalysis.slice(sim_var, time = t),
+                ClimaAnalysis.slice(obs_var, time = t),
+                mask = mask,
+            )
+            g_bias = ClimaAnalysis.global_bias(
+                ClimaAnalysis.slice(sim_var, time = t),
+                ClimaAnalysis.slice(obs_var, time = t),
+                mask = mask,
+            )
+            push!(rmse_vec, g_rmse)
+            push!(bias_vec, g_bias)
+        end
+        ax_rmse = CairoMakie.Axis(
+            fig[1, col],
+            title = "Global RMSE for $short_name",
+            xlabel = "Month",
+            ylabel = "Global RMSE ($(ClimaAnalysis.units(sim_var)))",
+            xticks = (1:num_times, Dates.monthname.(1:num_times) .|> x -> x[1:3]),
+        )
+        CairoMakie.lines!(ax_rmse, 1:num_times |> collect, rmse_vec)
+        CairoMakie.scatter!(ax_rmse, 1:num_times |> collect, rmse_vec)
+
+        ax_bias = CairoMakie.Axis(
+            fig[2, col],
+            title = "Global Bias for $short_name",
+            xlabel = "Month",
+            ylabel = "Global Bias ($(ClimaAnalysis.units(sim_var)))",
+            xticks = (1:num_times, Dates.monthname.(1:num_times) .|> x -> x[1:3]),
+        )
+        CairoMakie.lines!(ax_bias, 1:num_times |> collect, bias_vec)
+        CairoMakie.scatter!(ax_bias, 1:num_times |> collect, bias_vec)
+    end
+    CairoMakie.save(joinpath(leaderboard_base_path, "global_rmse_and_bias_graphs.png"), fig)
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    if length(ARGS) < 2
+        error(
+            "Usage: julia leaderboard.jl <Filepath to save leaderboard and bias plots> <Filepath to simulation data>",
+        )
+    end
+    leaderboard_base_path = ARGS[begin]
+    diagnostics_folder_path = ARGS[begin+1]
+    compute_leaderboard(leaderboard_base_path, diagnostics_folder_path)
+end