Extend historical solar forcing to cover the years from 1700 to 2300 (#99)

* Extend historical solar forcing from 1700 to 2300

* Fix an import statement

* Calculate tsi_mean correctly

* Use pre-industrial TSI for years before historical period

* Use a variable for the historical year range

Author: penguian

CMIP7/esm1p6/atmosphere/cmip7_ancil_constants.py

@@ -1,4 +1,5 @@
 from datetime import datetime
 
 ANCIL_TODAY = datetime.now().strftime("%Y.%m.%d")
+REAL_MISSING_DATA_INDICATOR = -32768 * 32768.0
 UM_VERSION = "7.3"

CMIP7/esm1p6/atmosphere/solar/cmip7_HI_solar_generate.py

@@ -1,14 +1,20 @@
 from argparse import ArgumentParser
 
 import iris
+import numpy as np
 from cmip7_ancil_argparse import common_parser
+from cmip7_ancil_constants import REAL_MISSING_DATA_INDICATOR
 from cmip7_HI import (
     CMIP7_HI_BEG_YEAR,
     CMIP7_HI_END_YEAR,
     esm_hi_forcing_save_dirpath,
 )
+from cmip7_PI import CMIP7_PI_YEAR
 from solar.cmip7_solar import cmip7_solar_dirpath, load_cmip7_solar_cube
 
+SOLAR_ARRAY_BEG_YEAR = 1700
+SOLAR_ARRAY_END_YEAR = 2300
+
 
 def parse_args():
     parser = ArgumentParser(
@@ -25,19 +31,52 @@ def parse_args():
     return parser.parse_args()
 
 
+def cmip7_hi_solar_year_mean(cube):
+    """
+    Calculate mean TSI values for each year and save them into an array.
+    """
+    NBR_YEARS = SOLAR_ARRAY_END_YEAR - SOLAR_ARRAY_BEG_YEAR + 1
+    solar_array = np.zeros(NBR_YEARS)
+    # Calculate and save the mean annual TSI for each CMIP7 historical year.
+    historical_year_range = range(CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR + 1)
+    assert CMIP7_PI_YEAR in historical_year_range
+    for year in historical_year_range:
+        year_cons = iris.Constraint(time=lambda cell: cell.point.year == year)
+        year_cube = cube.extract(year_cons)
+        year_mean = year_cube.collapsed("time", iris.analysis.MEAN).data
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = year_mean
+        # Save the year mean for the pre-industrial year.
+        if year == CMIP7_PI_YEAR:
+            pi_year_mean = year_mean
+
+    # For the years from SOLAR_ARRAY_BEG_YEAR to CMIP7_HI_BEG_YEAR - 1,
+    # set the saved TSI value to the pre-industrial year mean TSI.
+    for year in range(SOLAR_ARRAY_BEG_YEAR, CMIP7_HI_BEG_YEAR):
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = pi_year_mean
+
+    # For the years from CMIP7_HI_END_YEAR + 1 to SOLAR_ARRAY_END_YEAR,
+    # set the saved TSI value to the real missing data indicator
+    for year in range(CMIP7_HI_END_YEAR + 1, SOLAR_ARRAY_END_YEAR + 1):
+        solar_array[year - SOLAR_ARRAY_BEG_YEAR] = REAL_MISSING_DATA_INDICATOR
+    return solar_array
+
+
 def cmip7_hi_solar_save(args, cube):
+    """
+    Save the TSI values for each year into a text file.
+    """
+    solar_array = cmip7_hi_solar_year_mean(cube)
     save_dirpath = esm_hi_forcing_save_dirpath(args)
     # Ensure that the save directory exists.
     save_dirpath.mkdir(mode=0o755, parents=True, exist_ok=True)
     save_filepath = save_dirpath / args.save_filename
     with open(save_filepath, "w") as save_file:
-        for year in range(CMIP7_HI_BEG_YEAR, CMIP7_HI_END_YEAR + 1):
-            year_cons = iris.Constraint(
-                time=lambda cell: cell.point.year == year
-            )
-            year_cube = cube.extract(year_cons)
-            year_mean = year_cube.collapsed("time", iris.analysis.MEAN).data
-            print(year, f"{year_mean:.3f}", file=save_file)
+        for year in range(SOLAR_ARRAY_BEG_YEAR, SOLAR_ARRAY_END_YEAR + 1):
+            year_mean = solar_array[year - SOLAR_ARRAY_BEG_YEAR]
+            if year_mean == REAL_MISSING_DATA_INDICATOR:
+                print(year, f"{year_mean:.1f}", file=save_file)
+            else:
+                print(year, f"{year_mean:.3f}", file=save_file)
 
 
 if __name__ == "__main__":