force bucket and richards model with ERA5

experiments/Artifacts.toml

@@ -0,0 +1,15 @@
+[lwp_pallardy_etal2018]
+git-tree-sha1 = "29ed37c288601ab353f4f5dc13cf652ec4aec305"
+
+["processed_topographic_index 2.5 deg"]
+git-tree-sha1 = "31d4110d008fcfb1ecfd77eb9675c983b650eedf"
+
+[era5_land_forcing_data2021]
+git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a"
+
+[soil_params_Gupta2020_2022]
+git-tree-sha1 = "ffc4c3fa76ea632ba46ae5d58671d4161124af29"
+
+    [[soil_params_Gupta2020_2022.download]]
+    sha256 = "636db48ecb387a9aa59afbff8af86f0ccc092690c688f0cbb52041c133361e22"
+    url = "https://caltech.box.com/shared/static/f2y23qx0lggjskftzgh7ht7fsbh36gmm.gz"

experiments/Manifest.toml

@@ -336,10 +336,10 @@ uuid = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 version = "0.7.19"
 
 [[deps.ClimaUtilities]]
-deps = ["Artifacts", "CFTime", "Dates"]
-git-tree-sha1 = "d1b30a07829248b325ffffbb739d0f6fe79dfe7c"
+deps = ["Artifacts", "Dates"]
+git-tree-sha1 = "617ef474d78e2fd713626e22db8e9fb4881e90ee"
 uuid = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
-version = "0.1.3"
+version = "0.1.5"
 weakdeps = ["Adapt", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "Interpolations", "NCDatasets"]
 
     [deps.ClimaUtilities.extensions]

experiments/standalone/Bucket/global_bucket_staticmap.jl

@@ -20,9 +20,10 @@ using CairoMakie
 using Dates
 using DelimitedFiles
 using Statistics
-
+using ClimaUtilities.ClimaArtifacts
+import Interpolations
 import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
-
+import ClimaUtilities: @clima_artifact
 import ClimaTimeSteppers as CTS
 import NCDatasets
 using ClimaCore
@@ -55,7 +56,14 @@ function compute_extrema(v)
     return (minimum(mins), maximum(maxes))
 end
 
-anim_plots = false
+function compute_clims(v)
+    means = [mean(u) for u in v]
+    sigmas = [std(u) for u in v]
+    return (minimum(means) - maximum(sigmas), maximum(means) + maximum(sigmas))
+end
+
+anim_plots = true
+regridder_type = :InterpolationsRegridder
 FT = Float64;
 context = ClimaComms.context()
 earth_param_set = LP.LandParameters(FT);
@@ -70,11 +78,11 @@ soil_depth = FT(3.5);
 bucket_domain = ClimaLand.Domains.SphericalShell(;
     radius = FT(6.3781e6),
     depth = soil_depth,
-    nelements = (50, 10),
+    nelements = (30, 10),
     npolynomial = 1,
     dz_tuple = FT.((1.0, 0.05)),
 );
-ref_time = DateTime(2005);
+ref_time = DateTime(2021);
 
 # Set up parameters
 σS_c = FT(0.2);
@@ -85,43 +93,90 @@ z_0b = FT(1e-3);
 ρc_soil = FT(2e6);
 τc = FT(3600);
 t0 = 0.0;
-tf = 7 * 86400;
-Δt = 3600.0;
+tf = 14 * 86400;
+Δt = 3600.0 / 3;
 
 # Construct albedo parameter object using static map
 # Use separate regridding directory for CPU and GPU runs to avoid race condition
 device_suffix =
     typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ?
     "cpu" : "gpu"
-regrid_dir = joinpath(
-    pkgdir(ClimaLand),
-    "experiments/standalone/Bucket/$device_suffix/regrid-static/",
-)
-!ispath(regrid_dir) && mkpath(regrid_dir)
+t_start = t0
 surface_space = bucket_domain.space.surface
 α_snow = FT(0.8)
 albedo = PrescribedBaregroundAlbedo{FT}(α_snow, surface_space);
 
 bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc);
 
+# Forcing data
+era5_artifact_path = @clima_artifact("era5_land_forcing_data2021", context)
+# Below, the preprocess_func argument is used to
+# 1. Convert precipitation to be negative (as it is downwards)
+# 2. Convert accumulations over an hour to a rate per second
 # Precipitation:
-precip = (t) -> 0;
-snow_precip = (t) -> -5e-7 * (t < 1 * 86400);
-# Diurnal temperature variations:
-T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2);
-# Constant otherwise:
-u_atmos = (t) -> 3.0;
-q_atmos = (t) -> 0.001;
-h_atmos = FT(2);
-P_atmos = (t) -> 101325;
+precip = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
+    "rf",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+    file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+)
+
+snow_precip = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
+    "sf",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+    file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+)
+
+u_atmos = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
+    "ws",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+)
+q_atmos = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
+    "q",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+)
+P_atmos = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
+    "sp",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+)
+
+T_atmos = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
+    "t2m",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+)
+h_atmos = FT(10);
+
 
 bucket_atmos = PrescribedAtmosphere(
-    TimeVaryingInput(precip),
-    TimeVaryingInput(snow_precip),
-    TimeVaryingInput(T_atmos),
-    TimeVaryingInput(u_atmos),
-    TimeVaryingInput(q_atmos),
-    TimeVaryingInput(P_atmos),
+    precip,
+    snow_precip,
+    T_atmos,
+    u_atmos,
+    q_atmos,
+    P_atmos,
     ref_time,
     h_atmos,
 );
@@ -130,14 +185,26 @@ bucket_atmos = PrescribedAtmosphere(
 # peak at local noon, and a prescribed downwelling LW radiative
 # flux, assuming the air temperature is on average 275 degrees
 # K with a diurnal amplitude of 5 degrees K:
-SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0);
-LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4;
-bucket_rad = PrescribedRadiativeFluxes(
-    FT,
-    TimeVaryingInput(SW_d),
-    TimeVaryingInput(LW_d),
-    ref_time,
-);
+SW_d = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
+    "ssrd",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+    file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+)
+LW_d = TimeVaryingInput(
+    joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
+    "strd",
+    surface_space;
+    reference_date = ref_time,
+    t_start,
+    regridder_type,
+    file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+)
+
+bucket_rad = PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time);
 
 model = BucketModel(
     parameters = bucket_parameters,
@@ -214,18 +281,22 @@ F_sfc = [
     ) for k in 1:length(sol.t)
 ];
 
+sw_forcing = [
+    Remapping.interpolate(remapper, saved_values.saveval[k].drivers.SW_d)
+    for k in 1:length(sol.t)
+];
+
 # save prognostic state to CSV - for comparison between GPU and CPU output
 open(joinpath(outdir, "tf_state_$device_suffix.txt"), "w") do io
     writedlm(io, hcat(T_sfc[end][:], W[end][:], Ws[end][:], σS[end][:]), ',')
 end;
 # animation settings
-nframes = length(W)
-framerate = 2
+nframes = length(T_sfc) # hourly data
 fig_ts = Figure(size = (1000, 1000))
 for (i, (field_ts, field_name)) in enumerate(
     zip(
-        [W, σS, T_sfc, evaporation, F_sfc],
-        ["W", "σS", "T_sfc", "evaporation", "F_sfc"],
+        [W, σS, T_sfc, evaporation, F_sfc, sw_forcing],
+        ["W", "σS", "T_sfc", "evaporation", "F_sfc", "SW forcing"],
     ),
 )
     if anim_plots
@@ -236,13 +307,18 @@ for (i, (field_ts, field_name)) in enumerate(
             ylabel = "Latitude",
             title = field_name,
         )
-        clims = compute_extrema(field_ts)
+        clims = compute_clims(field_ts)
         CairoMakie.Colorbar(fig[:, end + 1], colorrange = clims)
         outfile = joinpath(
             outdir,
             string("anim_$(device_suffix)_", field_name, ".mp4"),
         )
-        record(fig, outfile, 1:nframes; framerate = framerate) do frame
+        record(
+            fig,
+            outfile,
+            (nframes - 7 * 24):2:nframes;
+            framerate = 3,
+        ) do frame
             CairoMakie.heatmap!(
                 longpts,
                 latpts,
@@ -251,18 +327,7 @@ for (i, (field_ts, field_name)) in enumerate(
             )
         end
     end
-    # Plot the timeseries of the mean value as well.
-    xlabel = i == 5 ? "Time (days)" : ""
-    ax2 = Axis(
-        fig_ts[i, 1],
-        xlabel = xlabel,
-        ylabel = field_name,
-        title = "Global bucket with static map albedo",
-    )
-    CairoMakie.lines!(ax2, sol.t ./ 3600 ./ 24, [mean(x) for x in field_ts])
+
 end
 outfile = joinpath(outdir, string("ts_$device_suffix.png"))
 CairoMakie.save(outfile, fig_ts)
-
-# delete regrid directory
-rm(regrid_dir, recursive = true)