Calibrate global bucket LHF and SHF

.buildkite/Project.toml

@@ -6,6 +6,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
+ClimaCalibrate = "4347a170-ebd6-470c-89d3-5c705c0cacc2"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
 ClimaCoreMakie = "908f55d8-4145-4867-9c14-5dad1a479e4d"
@@ -17,6 +18,7 @@ ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"

Artifacts.toml

@@ -141,3 +141,6 @@ git-tree-sha1 = "4cc329093b0dbc3443db36e7c01133267ac94d37"
     [[bedrock_depth_60arcseconds.download]]
     sha256 = "6eb359e900141de837c41c0082257f7cf082e06e1371d45b07b1078f07e96c53"
     url = "https://caltech.box.com/shared/static/k1zq0w1psqu8dk0s9jotrfuu6hm8kxfp.gz"
+
+[era5_surface_fluxes]
+git-tree-sha1 = "f2cb69e41cba84991054162e794292f98bf0ae34"

derecho_job.pbs

@@ -0,0 +1,18 @@
+#!/bin/bash
+#PBS -N derecho_calibration
+#PBS -o output.txt
+#PBS -e error.txt
+#PBS -l walltime=10:00:00
+#PBS -l select=1:ncpus=64:ngpus=4
+
+## Account number for CliMA
+#PBS -A UCIT0011
+#PBS -q main
+
+module load climacommon
+# Run your command
+export CLIMACOMMS_DEVICE="CUDA"
+export CLIMACOMMS_CONTEXT="SINGLETON"
+julia --project=.buildkite -e 'using Pkg; Pkg.instantiate(;verbose=true)'
+
+julia --project=.buildkite/ experiments/calibration/global_land/calibrate_land.jl

experiments/calibration/global_bucket/calibrate_bucket.jl

@@ -0,0 +1,71 @@
+# Use the .buildkite environment
+# bucket_turbulent_fluxes(params) returns lhf and shf as a function of 6 parameters
+include("calibrate_bucket_function.jl")
+
+using Random
+
+#= Target (perfect model)
+params = (;
+    κ_soil = FT(1.5),
+    ρc_soil = FT(2e6),
+    f_bucket = FT(0.75),
+    W_f = FT(0.2),
+    p = FT(1),
+    z_0m = FT(1e-2),
+)
+target = bucket_turbulent_fluxes(params)[1]
+=#
+
+target = full_obs_era5
+
+# Parameters prior
+prior_κ_soil = EKP.constrained_gaussian("κ_soil", 2, 1, 0, Inf);
+prior_ρc_soil = EKP.constrained_gaussian("ρc_soil", 4e6, 2e6, 0, Inf);
+prior_f_bucket = EKP.constrained_gaussian("f_bucket", 0.5, 0.3, 0, 1);
+prior_W_f = EKP.constrained_gaussian("W_f", 0.4, 0.4, 0, Inf);
+prior_p = EKP.constrained_gaussian("p", 2, 1, 1, Inf);
+prior_z_0m = EKP.constrained_gaussian("z_0m", 0.01, 0.1, 0, Inf);
+prior = EKP.combine_distributions([
+    prior_κ_soil,
+    prior_ρc_soil,
+    prior_f_bucket,
+    prior_W_f,
+    prior_p,
+    prior_z_0m,
+]);
+
+rng_seed = 2
+rng = Random.MersenneTwister(rng_seed)
+
+N_ensemble = 10
+N_iterations = 5
+
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble);
+ensemble_kalman_process = EKP.EnsembleKalmanProcess(
+    initial_ensemble,
+    target,
+    EKP.Inversion();
+    rng = rng,
+);
+
+for i in 1:N_iterations
+    params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process)
+
+    G_ens = hcat(
+        [
+            bucket_turbulent_fluxes((;
+                κ_soil = params_i[1, j],
+                ρc_soil = params_i[2, j],
+                f_bucket = params_i[3, j],
+                W_f = params_i[4, j],
+                p = params_i[5, j],
+                z_0m = params_i[6, j],
+            ))[2] for j in 1:N_ensemble
+        ]...,
+    )
+
+    EKP.update_ensemble!(ensemble_kalman_process, G_ens)
+end
+
+final_ensemble = EKP.get_ϕ_final(prior, ensemble_kalman_process)
+println(final_ensemble)