prep to calibrate model with new IC [skip ci]

Author: kmdeck

Artifacts.toml

@@ -1,3 +1,10 @@
+[saturated_land_ic]
+git-tree-sha1 = "cc2888b590f5158113b54c75eee3e2f797be8956"
+
+    [[saturated_land_ic.download]]
+    sha256 = "d73ec9bf8a20f1c63a0133d483a08e801db9125fda43ab51d080bc0192343f84"
+    url = "https://caltech.box.com/shared/static/99j3juuj3vnnzy16ye31dsfcx25s71z6.gz"
+
 [landsea_mask_30arcseconds]
 git-tree-sha1 = "517c925535981f9d17ac9e7a65e2c35c3ce9597d"
 

experiments/calibration/models/snowy_land.jl

@@ -63,8 +63,24 @@ function setup_model(
         scf,
     )
 
+    ground = ClimaLand.PrognosticGroundConditions{FT}()
+    canopy_forcing = (; atmos, radiation, ground)
+    photosynthesis = PModel{FT}(domain, toml_dict)
+    conductance = PModelConductance{FT}(toml_dict)
+    soil_moisture_stress =
+        ClimaLand.Canopy.PiecewiseMoistureStressModel{FT}(domain, toml_dict)
+    canopy = ClimaLand.Canopy.CanopyModel{FT}(
+        surface_domain,
+        canopy_forcing,
+        LAI,
+        toml_dict;
+        prognostic_land_components = (:canopy, :snow, :soil, :soilco2),
+        photosynthesis,
+        conductance,
+        soil_moisture_stress,
+    )
     # Construct the land model with all default components except for snow
-    land = LandModel{FT}(forcing, LAI, toml_dict, domain, Δt; snow)
+    land = LandModel{FT}(forcing, LAI, toml_dict, domain, Δt; snow, canopy)
     return land
 end
 

experiments/calibration/run_calibration.jl

@@ -11,24 +11,30 @@ import JLD2
 include(joinpath(pkgdir(ClimaLand), "experiments/calibration/api.jl"))
 
 const CALIBRATE_CONFIG = CalibrateConfig(;
-    short_names = ["lwu"],
+    short_names = ["lwu", "shf", "lhf"],
     minibatch_size = 1,
-    n_iterations = 1,
-    sample_date_ranges = [("2007-12-1", "2007-12-1")],
+    n_iterations = 10,
+    sample_date_ranges = [
+        ("$(2000 + 2*i)-12-1", "$(2002 + 2*i)-9-1") for i in 0:9
+    ], # 2000 to 2020
     extend = Dates.Month(3),
-    spinup = Dates.Month(0),
+    spinup = Dates.Month(3),
     nelements = (180, 360, 15),
     output_dir = "experiments/calibration/land_model",
     rng_seed = 42,
-    obs_vec_filepath = "experiments/calibration/land_observation_vector.jld2",
+    obs_vec_filepath = "/glade/derecho/scratch/kdeck/recalibrate_saturated",
     model_type = ClimaLand.LandModel,
 )
 
 
 if abspath(PROGRAM_FILE) == @__FILE__
-    # true solution is at 0.96
-    priors =
-        [EKP.constrained_gaussian("emissivity_bare_soil", 0.82, 0.12, 0.0, 2.0)]
+    priors = [
+        EKP.constrained_gaussian("moisture_stress_c", 1.0, 0.5, 0, 2),
+        EKP.constrained_gaussian("pmodel_cstar", 0.41, 0.11, 0, Inf),
+        EKP.constrained_gaussian("pmodel_β", 146, 10, 0, Inf),
+        EKP.constrained_gaussian("canopy_u0", 100, 90, 0, Inf),
+        EKP.constrained_gaussian("ac_canopy", 3.5e4, 2e4, 2.5e3, 2.5e5),
+    ]
     prior = EKP.combine_distributions(priors)
 
     observation_vector = JLD2.load_object(CALIBRATE_CONFIG.obs_vec_filepath)

experiments/long_runs/snowy_land_pmodel.jl

@@ -69,6 +69,7 @@ function setup_model(
     surface_domain = ClimaLand.Domains.obtain_surface_domain(domain)
     surface_space = domain.space.surface
     # Forcing data - high resolution
+    use_lowres_forcing = true
     atmos, radiation = ClimaLand.prescribed_forcing_era5(
         start_date,
         stop_date,
@@ -77,6 +78,7 @@ function setup_model(
         FT;
         max_wind_speed = 25.0,
         context,
+        use_lowres_forcing,
     )
     forcing = (; atmos, radiation)
 

src/Artifacts.jl

@@ -6,6 +6,11 @@ import ClimaUtilities.ClimaArtifacts: @clima_artifact
 
 import LazyArtifacts
 
+function saturated_land_ic_path(; context = nothing)
+    dir = @clima_artifact("saturated_land_ic", context)
+    return joinpath(dir, "saturated_land_ic.nc")
+end
+
 """
     soil_ic_2008_50m_path(; context)
 

src/shared_utilities/Domains.jl

@@ -1247,8 +1247,8 @@ end
         apply_mask = true,
         mask_threshold = 0.5,
         nelements = (101, 15),
-        dz_tuple = (10.0, 0.05),
-        depth = 50.0,
+        dz_tuple = (3.0, 0.05),
+        depth = 15.0,
         npolynomial = 0,
         context = ClimaComms.context(),
         filepath = landseamask_file_path(;context),

src/simulations/Simulations.jl

@@ -70,7 +70,7 @@ end
         model;
         outdir = ".",
         set_ic! = make_set_initial_state_from_file(
-            ClimaLand.Artifacts.soil_ic_2008_50m_path(;
+            ClimaLand.Artifacts.saturated_land_ic_path(;
                 context = ClimaComms.context(model),
             ),
             model,
@@ -116,7 +116,7 @@ function LandSimulation(
     model;
     outdir = ".",
     set_ic! = make_set_initial_state_from_file(
-        ClimaLand.Artifacts.soil_ic_2008_50m_path(;
+        ClimaLand.Artifacts.saturated_land_ic_path(;
             context = ClimaComms.context(model),
         ),
         model,

src/simulations/initial_conditions.jl

@@ -49,11 +49,11 @@ function set_soil_initial_conditions!(
         regridder_kwargs = (; extrapolation_bc, interpolation_method),
     )
 
-    Y.soil.ϑ_l .= enforce_residual_constraint.(Y.soil.ϑ_l, θ_r)
-    Y.soil.ϑ_l .= enforce_porosity_constraint.(Y.soil.ϑ_l, ν)
-    Y.soil.θ_i .=
-        enforce_residual_constraint.(Y.soil.θ_i, eltype(Y.soil.θ_i)(0))
-    Y.soil.θ_i .= enforce_porosity_constraint.(Y.soil.ϑ_l, Y.soil.θ_i, ν)
+#    Y.soil.ϑ_l .= enforce_residual_constraint.(Y.soil.ϑ_l, θ_r)
+#    Y.soil.ϑ_l .= enforce_porosity_constraint.(Y.soil.ϑ_l, ν)
+#    Y.soil.θ_i .=
+#        enforce_residual_constraint.(Y.soil.θ_i, eltype(Y.soil.θ_i)(0))
+#    Y.soil.θ_i .= enforce_porosity_constraint.(Y.soil.ϑ_l, Y.soil.θ_i, ν)
     ρc_s =
         ClimaLand.Soil.volumetric_heat_capacity.(
             Y.soil.ϑ_l,
@@ -68,21 +68,21 @@ function set_soil_initial_conditions!(
         regridder_type,
         regridder_kwargs = (; extrapolation_bc, interpolation_method),
     )
-    T =
-        ClimaLand.Soil.temperature_from_ρe_int.(
-            Y.soil.ρe_int,
-            Y.soil.θ_i,
-            ρc_s,
-            soil.parameters.earth_param_set,
-        )
-    T .= clip_to_bounds.(T, T_bounds[1], T_bounds[2])
-    Y.soil.ρe_int .=
-        ClimaLand.Soil.volumetric_internal_energy.(
-            Y.soil.θ_i,
-            ρc_s,
-            T,
-            soil.parameters.earth_param_set,
-        )
+#    T =
+#        ClimaLand.Soil.temperature_from_ρe_int.(
+#            Y.soil.ρe_int,
+#            Y.soil.θ_i,
+#            ρc_s,
+#            soil.parameters.earth_param_set,
+#        )
+#    T .= clip_to_bounds.(T, T_bounds[1], T_bounds[2])
+#    Y.soil.ρe_int .=
+#        ClimaLand.Soil.volumetric_internal_energy.(
+#            Y.soil.θ_i,
+#            ρc_s,
+#            T,
+#            soil.parameters.earth_param_set,
+#        )
     return nothing
 end
 
@@ -247,20 +247,12 @@ function make_set_initial_state_from_file(
         # to soil potential (soil moisture), averaged over the soil layers,
         # which would correspond to approximate steady state
         if land.canopy.hydraulics isa ClimaLand.Canopy.PlantHydraulicsModel
-            @. p.soil.ψ = ClimaLand.Soil.pressure_head(
-                land.soil.parameters.hydrology_cm,
-                land.soil.parameters.θ_r,
-                Y.soil.ϑ_l,
-                land.soil.parameters.ν - Y.soil.θ_i,
-                land.soil.parameters.S_s,
-            )
-            ψ_roots = ClimaCore.Fields.zeros(axes(Y.canopy.hydraulics.ϑ_l.:1))
-            z = land.soil.domain.fields.z
-            tmp = @. ClimaLand.Canopy.root_distribution(
-                z,
-                land.canopy.biomass.rooting_depth,
-            ) * p.soil.ψ / land.soil.domain.fields.depth
-            ClimaCore.Operators.column_integral_definite!(ψ_roots, tmp)
+            ψ_roots = SpaceVaryingInput(
+                ic_path,
+                "lwp",
+                land.snow.domain.space.surface;
+                regridder_type,
+                regridder_kwargs = (; extrapolation_bc, interpolation_method))
             Y.canopy.hydraulics.ϑ_l.:1 .=
                 ClimaLand.Canopy.PlantHydraulics.inverse_water_retention_curve.(
                     land.canopy.hydraulics.parameters.retention_model,

src/standalone/Soil/Soil.jl

@@ -278,7 +278,10 @@ function EnergyHydrology{FT}(
         runoff;
         prognostic_land_components,
     )
-    bottom_bc = EnergyWaterFreeDrainage()
+    bottom_water_flux = WaterFluxBC((p, t) -> 0.0)
+    bottom_heat_flux = HeatFluxBC((p, t) -> 0.0)
+    bottom_bc =
+        WaterHeatBC(; water = bottom_water_flux, heat = bottom_heat_flux)
     boundary_conditions = (; top = top_bc, bottom = bottom_bc)
     # sublimation and subsurface runoff are added automatically
     sources = (additional_sources..., PhaseChange{FT}())