wip [skip ci]

experiments/integrated/global/global_parameters.jl

@@ -0,0 +1,201 @@
+# Read in f_max data and land sea mask
+topmodel_dataset = ArtifactWrapper(
+    @__DIR__,
+    "processed_topographic_index 2.5 deg",
+    ArtifactFile[ArtifactFile(
+        url = "https://caltech.box.com/shared/static/dwa7g0uzhxd50a2z3thbx3a12n0r887s.nc",
+        filename = "means_2.5_new.nc",
+    ),],
+)
+infile_path = joinpath(get_data_folder(topmodel_dataset), "means_2.5_new.nc")
+outfile_root =
+    joinpath(pkgdir(ClimaLand), "experiments/standalone/Soil/static_data_cgll")
+
+f_max = SpaceVaryingInput(
+    infile_path,
+    "fmax",
+    surface_space;
+    regridder_type = :InterpolationsRegridder,
+)
+mask = SpaceVaryingInput(
+    infile_path,
+    "landsea_mask",
+    surface_space;
+    regridder_type = :InterpolationsRegridder,
+)
+
+oceans_to_zero(field, mask) = mask > 0.5 ? field : eltype(field)(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,
+)
+soil_params_artifact_path = "/groups/esm/ClimaArtifacts/artifacts/soil_params_Gupta2020_2022"# @clima_artifact("soil_params_Gupta2020_2022")
+function mask_vg(var, value)
+    if var < 1e-8
+        return value
+    else
+        return var
+    end
+end
+vg_α = SpaceVaryingInput(
+    joinpath(
+        soil_params_artifact_path,
+        "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc",
+    ),
+    "α",
+    subsurface_space;
+    regridder_type = :InterpolationsRegridder,
+    regridder_kwargs = (;
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    ),
+)
+vg_α .= mask_vg.(vg_α, 1e-3)
+function mask_vg_n(var, value)
+    if var < 1
+        return value
+    else
+        return var
+    end
+end
+vg_n = SpaceVaryingInput(
+    joinpath(soil_params_artifact_path, "vGn_map_gupta_etal2020_2.5x2.5x4.nc"),
+    "n",
+    subsurface_space;
+    regridder_type = :InterpolationsRegridder,
+    regridder_kwargs = (;
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    ),
+)
+vg_n .= mask_vg_n.(vg_n, 1.001)
+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_2.5x2.5x4.nc",
+    ),
+    "θ_r",
+    subsurface_space;
+    regridder_type = :InterpolationsRegridder,
+    regridder_kwargs = (;
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    ),
+)
+
+ν = SpaceVaryingInput(
+    joinpath(
+        soil_params_artifact_path,
+        "porosity_map_gupta_etal2020_2.5x2.5x4.nc",
+    ),
+    "ν",
+    subsurface_space;
+    regridder_type = :InterpolationsRegridder,
+    regridder_kwargs = (;
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    ),
+)
+ν .= mask_vg.(ν, 1.0)
+K_sat = SpaceVaryingInput(
+    joinpath(soil_params_artifact_path, "ksat_map_gupta_etal2020_2.5x2.5x4.nc"),
+    "Ksat",
+    subsurface_space;
+    regridder_type = :InterpolationsRegridder,
+    regridder_kwargs = (;
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+    ),
+)
+S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
+ν_ss_om = ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1)
+ν_ss_quartz = ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1)
+ν_ss_gravel = ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1)
+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,
+    NIR_albedo,
+);
+
+
+# TwoStreamModel parameters
+Ω = FT(0.69)
+ld = FT(0.5)
+α_PAR_leaf = FT(0.1)
+τ_PAR_leaf = FT(0.05)
+α_NIR_leaf = FT(0.45)
+τ_NIR_leaf = FT(0.25)
+
+# Energy Balance model
+ac_canopy = FT(2.5e3)
+
+# Conductance Model
+g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300.
+
+#Photosynthesis model
+Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5
+
+# Plant Hydraulics and general plant parameters
+SAI = FT(1.0) # m2/m2 or: estimated from Wang et al, FT(0.00242) ?
+f_root_to_shoot = FT(3.5)
+RAI = f_root_to_shoot * SAI
+K_sat_plant = 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 = FT(0.5) # from Natan
+n_stem = 1
+n_leaf = 1
+h_stem = 1.0
+h_leaf = 1.0
+zmax = 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;
+    ν = 1.0,# INCORRECT!
+)