Merge pull request #123 from chrisroadmap/eclipse

Eclipse

input/fair-2.1.3/v1.4/AR6-updated_no-contrails_eclipse/calibration/01_tabulate-4xCO2.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+"""Combine Hege-Beate's data into one file."""
+
+import glob
+import os
+from pathlib import PurePath
+
+import pandas as pd
+from dotenv import load_dotenv
+from fair import __version__
+
+load_dotenv()
+
+print("Making nice 4xCO2 data...")
+
+cal_v = os.getenv("CALIBRATION_VERSION")
+fair_v = os.getenv("FAIR_VERSION")
+constraint_set = os.getenv("CONSTRAINT_SET")
+
+assert fair_v == __version__
+
+available_files = glob.glob(
+    "../../../../../data/cmip6-hbf/cmip_data/*/abrupt-4xCO2/"
+    "*_abrupt-4xCO2_*_anomalies.txt"
+)
+
+maxlen = 0
+
+models = []
+runs = []
+lines = []
+for file in available_files:
+    model = PurePath(file).parts[8]
+    run = PurePath(file).parts[10].split("_")[2]
+    models.append(model)
+    runs.append(run)
+    df = pd.read_csv(file, index_col=0)
+    vars = {}
+    for var in ["tas", "rlut", "rsut", "rsdt"]:
+        vars[var] = df[var].values.squeeze()
+        if len(vars[var]) > maxlen:
+            maxlen = len(vars[var])
+        line = [
+            "CMIP",
+            model,
+            run,
+            "CMIP6",
+            "unspecified",
+            "World",
+            "abrupt-4xCO2",
+            "W m^-2",
+            var,
+        ]
+        line.extend(vars[var][:150])
+        lines.append(line)
+    vars["rndt"] = vars["rsdt"] - vars["rsut"] - vars["rlut"]
+    line = [
+        "CMIP",
+        model,
+        run,
+        "CMIP6",
+        "unspecified",
+        "World",
+        "abrupt-4xCO2",
+        "W m^-2",
+        "rndt",
+    ]
+    line.extend(vars["rndt"][:150])
+    lines.append(line)
+
+df = pd.DataFrame(
+    lines,
+    columns=(
+        [
+            "activity_id",
+            "climate_model",
+            "member_id",
+            "mip_era",
+            "model",
+            "region",
+            "scenario",
+            "unit",
+            "variable",
+        ]
+        + ["X%d" % year for year in range(1850, 2000)]
+    ),
+)
+
+to_remove = [f"X{year}" for year in range(1850, 2000)]
+df.dropna(subset=to_remove, inplace=True)
+
+os.makedirs(
+    f"../../../../../output/fair-{fair_v}/v{cal_v}/{constraint_set}/calibrations/",
+    exist_ok=True,
+)
+
+df.to_csv(
+    f"../../../../../output/fair-{fair_v}/v{cal_v}/{constraint_set}/"
+    "calibrations/4xCO2_cmip6.csv",
+    index=False,
+)

input/fair-2.1.3/v1.4/AR6-updated_no-contrails_eclipse/calibration/02_calibrate_cummins_3layer.r

@@ -0,0 +1,148 @@
+#!/usr/bin/env Rscript
+
+# Goes through each of the models in turn and tunes the parameters of the
+# Cummins three layer model.
+#
+# It will produce an output csv table of Cummins three layer parameters, which
+# we will then put in impulse-response form for FaIR.
+#
+# References:
+# Cummins, Cummins, D. P., Stephenson, D. B., & Stott, P. A. (2020). Optimal
+# Estimation of Stochastic Energy Balance Model Parameters, Journal of Climate,
+# 33(18), 7909-7926, https://doi.org/10.1175/JCLI-D-19-0589.1
+#
+# donaldcummins. (2021). donaldcummins/EBM: Optional quadratic penalty
+# (v1.1.0). Zenodo. https://doi.org/10.5281/zenodo.5217975
+
+message("Running R script for 3 layer model calibrations...")
+
+# Use Donald Cummins' package
+library(EBM)
+
+# Get environment variable describing calibration version
+readRenviron("../../../../../.env")
+cal_v = paste("v", Sys.getenv("CALIBRATION_VERSION"), sep="")
+fair_v = paste("fair-", Sys.getenv("FAIR_VERSION"), sep="")
+constraint_set = Sys.getenv("CONSTRAINT_SET")
+
+# Get the precalculated 4xCO2 N and T data
+input_data = read.csv(file.path("..", "..", "..", "..", "..", "output", fair_v, cal_v, constraint_set, "calibrations", "4xCO2_cmip6.csv"))
+
+# Initial guess for parameter values
+inits3 <- list(
+	gamma = 2,
+	C = c(4, 15, 80),
+	kappa = c(1, 2, 1),
+	epsilon = 1.1,
+	sigma_eta = 0.5,
+	sigma_xi = 0.5,
+	F_4xCO2 = 8
+)
+
+# prepare empty output dataframe
+output <- data.frame(matrix(ncol = 15, nrow = 0))
+
+# grab models
+models = unique(input_data$climate_model)
+
+# iterate through models
+for (model in models) {
+	runs <- input_data[
+		(input_data$climate_model == model) &
+		(input_data$variable == "tas"), 3
+	]
+
+	# iterate through runs of the same model for now, though we will probably
+	# want to combine/downweight this in the final edition
+	for (run in runs) {
+		rndt <- input_data[
+			(input_data$climate_model == model) &
+			(input_data$variable == "rndt") &
+			(input_data$member_id == run), 10:159
+		]
+		rndt <- unname(rndt)
+
+		tas <- input_data[
+			(input_data$climate_model == model) &
+			(input_data$variable ==  "tas") &
+			(input_data$member_id == run), 10:159
+		]
+		tas <- unname(tas)
+
+		tas <- unlist(tas)
+		rndt <- unlist(rndt)
+
+        # Try and achieve convergence: if it fails, use a different alpha and try again.
+		# If it still fails, give up and move on.
+		attempt <- 0
+		success <- FALSE
+		while (attempt < 5) {
+			tryCatch(
+				{
+					check <- capture.output(
+						result <- FitKalman(
+							inits3,
+							T1 = tas,
+							N = rndt,
+							alpha = 1e-05 * 10^attempt,
+							maxeval=20000)$p
+						)
+					success <- TRUE
+					print(check)
+				},
+				error = function(c) {
+					message(paste(model, "did not converge or ran out of iterations"))
+					message("Here's the original error message from FitKalman:")
+					message(c)
+					if (attempt < 4) {
+						message(paste("Trying again with more liberal quadratic penalty (alpha) = ", 1e-05 * 10^attempt))
+					}
+					else {
+						message(paste("This isn't going to work, is it? I'm giving up for", model))
+					}
+				}
+			)
+			if (success) { break }
+			attempt <- attempt + 1
+		}
+		if (!success) {
+			message(paste("I am excluding" , model, " from my table of results."))
+			next
+		}
+
+		conv <- TRUE
+		nit <- strsplit(check, " ")[[1]][5]
+
+		row_out <- c(
+			model,
+			run,
+			conv,
+			nit,
+			result$gamma,
+			result$C[1],
+			result$C[2],
+			result$C[3],
+			result$kappa[1],
+			result$kappa[2],
+			result$kappa[3],
+			result$epsilon,
+			result$sigma_eta,
+			result$sigma_xi,
+			result$F_4xCO2
+		)
+
+		output <- rbind(output, row_out)
+	}
+}
+
+# rename columns away from random defaults
+names = c("model", "run", "conv", "nit", "gamma", "C1", "C2", "C3", "kappa1",
+	"kappa2", "kappa3", "epsilon", "sigma_eta", "sigma_xi", "F_4xCO2")
+colnames(output) <- names
+
+# save output
+write.csv(
+	output,
+	file.path("..", "..", "..", "..", "..", "output", fair_v, cal_v, constraint_set, "calibrations", "4xCO2_cummins_ebm3_cmip6.csv"),
+	row.names=FALSE
+)