Spatial parameter gap filling; soil longer run; bug with texture parameters (#1165)

Author: kmdeck

.buildkite/longruns_gpu/pipeline.yml

@@ -49,7 +49,7 @@ steps:
           - "snowy_land_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 15:00:00
+          slurm_time: 3:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
@@ -69,7 +69,7 @@ steps:
         artifact_paths: "soil_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 15:00:00
+          slurm_time: 3:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
@@ -83,18 +83,29 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
-  - group: "Very long run of snowy land"
+  - group: "Longer runs of Global Land Models"
     if: build.env("LONGER_RUN") != null
     steps:
-      - label: ":snow_capped_mountain: Snowy Land"
+      - label: "Snowy Land, 10 years"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/snowy_land.jl
         artifact_paths:
           - "snowy_land_longrun_gpu/*png"
           - "snowy_land_longrun_gpu/*pdf"
         agents:
           slurm_gpus: 1
-          slurm_time: 24:00:00
+          slurm_time: 15:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
           LONGER_RUN: ""
+
+      - label: "Soil, 10 years"
+        command:
+          - julia --color=yes --project=.buildkite experiments/long_runs/soil.jl
+        artifact_paths: "soil_longrun_gpu/*pdf"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 15:00:00
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+          LONGER_RUN: ""

experiments/long_runs/soil.jl

@@ -43,7 +43,12 @@ import NCDatasets
 using Poppler_jll: pdfunite
 
 const FT = Float64;
-time_interpolation_method = LinearInterpolation(PeriodicCalendar())
+# If you want to do a very long run locally, you can enter `export
+# LONGER_RUN=""` in the terminal and run this script. If you want to do a very
+# long run on Buildkite manually, then make a new build and pass `LONGER_RUN=""`
+# as an environment variable. In both cases, the value of `LONGER_RUN` does not
+# matter.
+const LONGER_RUN = haskey(ENV, "LONGER_RUN") ? true : false
 context = ClimaComms.context()
 ClimaComms.init(context)
 device = ClimaComms.device()
@@ -53,12 +58,23 @@ diagnostics_outdir = joinpath(root_path, "global_diagnostics")
 outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
 
-function setup_model(FT, start_date, domain, earth_param_set)
+function setup_model(FT, start_date, stop_date, domain, earth_param_set)
+    time_interpolation_method =
+        LONGER_RUN ? LinearInterpolation() :
+        LinearInterpolation(PeriodicCalendar())
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
     # Forcing data
-    era5_ncdata_path =
-        ClimaLand.Artifacts.era5_land_forcing_data2008_path(; context)
+    if LONGER_RUN
+        era5_ncdata_path = ClimaLand.Artifacts.find_era5_year_paths(
+            start_date,
+            stop_date;
+            context,
+        )
+    else
+        era5_ncdata_path =
+            ClimaLand.Artifacts.era5_land_forcing_data2008_path(; context)
+    end
     atmos, radiation = ClimaLand.prescribed_forcing_era5(
         era5_ncdata_path,
         surface_space,
@@ -133,8 +149,10 @@ function setup_model(FT, start_date, domain, earth_param_set)
 end
 
 function setup_simulation(; greet = false)
-    start_date = DateTime(2008)
-    stop_date = DateTime(2010)
+    # If not LONGER_RUN, run for 2 years; note that the forcing from 2008 is repeated.
+    # If LONGER run, run for 10 years, with the correct forcing each year.
+    start_date = LONGER_RUN ? DateTime(2004) : DateTime(2008)
+    stop_date = LONGER_RUN ? DateTime(2014) : DateTime(2010)
     Δt = 450.0
     nelements = (101, 15)
     if greet
@@ -147,7 +165,7 @@ function setup_simulation(; greet = false)
     domain =
         ClimaLand.ModelSetup.global_domain(FT; comms_ctx = context, nelements)
     params = LP.LandParameters(FT)
-    model = setup_model(FT, start_date, domain, params)
+    model = setup_model(FT, start_date, stop_date, domain, params)
     diagnostics = ClimaLand.default_diagnostics(
         model,
         start_date;

src/simulations/spatial_parameters.jl

@@ -280,14 +280,17 @@ function default_spatially_varying_soil_parameters(
         file_reader_kwargs = (; preprocess_func = (data) -> data > 0,),
     )
     # If the parameter mask =  0, set to physical value
+    # (equal to the mean where we have data; the mean of α is in log space)
     # This function in applied in **simulation mask** aware
     # manner.
+    # That is, we replace values in the simulation, but without data, with the mean
+    # over the data.
     masked_to_value(field, mask, value) =
         mask == 1.0 ? field : eltype(field)(value)
 
-    μ = FT(0.27)
+    μ = FT(0.33)
     vg_α .= masked_to_value.(vg_α, soil_params_mask, 10.0^μ)
-    μ = FT(1.65)
+    μ = FT(1.74)
     vg_n .= masked_to_value.(vg_n, soil_params_mask, μ)
 
     vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
@@ -326,20 +329,15 @@ function default_spatially_varying_soil_parameters(
         regridder_kwargs = (; extrapolation_bc,),
     )
 
-    μ = FT(-6)
+    # Set missing values to the mean. For Ksat, we use the mean in log space.
+    μ = FT(-5.08)
     K_sat .= masked_to_value.(K_sat, soil_params_mask, 10.0^μ)
 
-    ν .= masked_to_value.(ν, soil_params_mask, 1)
+    ν .= masked_to_value.(ν, soil_params_mask, 0.47)
 
-    θ_r .= masked_to_value.(θ_r, soil_params_mask, 0)
+    θ_r .= masked_to_value.(θ_r, soil_params_mask, 0.09)
 
-
-    S_s =
-        masked_to_value.(
-            ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3),
-            soil_params_mask,
-            1,
-        )
+    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
 
     soilgrids_artifact_path =
         ClimaLand.Artifacts.soil_grids_params_artifact_path(;
@@ -372,10 +370,11 @@ function default_spatially_varying_soil_parameters(
     # we require that the sum of these be less than 1 and equal to or bigger than zero.
     # The input should satisfy this almost exactly, but the regridded values may not.
     # Values of zero are OK here, so we dont need to apply any masking
-    texture_norm = @. min(ν_ss_gravel + ν_ss_quartz + ν_ss_om, 1)
-    @. ν_ss_gravel = ν_ss_gravel / max(texture_norm, eps(FT))
-    @. ν_ss_om = ν_ss_om / max(texture_norm, eps(FT))
-    @. ν_ss_quartz = ν_ss_quartz / max(texture_norm, eps(FT))
+    # if sum > 1, normalize by sum, else "normalize" by 1 (i.e., do not normalize)
+    texture_norm = @. max(ν_ss_gravel + ν_ss_quartz + ν_ss_om, 1)
+    @. ν_ss_gravel = ν_ss_gravel / texture_norm
+    @. ν_ss_om = ν_ss_om / texture_norm
+    @. ν_ss_quartz = ν_ss_quartz / texture_norm
 
     # Read in f_max data and topmodel params land sea mask
     infile_path = ClimaLand.Artifacts.topmodel_data_path()