Merge pull request #47 from GEOS-ESM/develop

GitFlow: Merge develop into main

CMakeLists.txt

@@ -10,5 +10,5 @@ set (alldirs
 esma_add_library (${this}
   SRCS GEOS_RadiationGridComp.F90
   SUBCOMPONENTS ${alldirs}
-  DEPENDENCIES MAPL GEOS_Shared esmf)
+  DEPENDENCIES MAPL GEOS_Shared ESMF::ESMF)
 

GEOS_RadiationShared/NRLSSI2/Mg_SB_from_daily_file.py

@@ -0,0 +1,157 @@
+'''
+Get Solar Mg and SB indices from NRLSSI2 daily auxilliary text file
+IMPORTANT: The filenames list below will be read in order, and any date that
+already has final data will not be overwriten by data for that date in a later
+file. For this reason, the filenames list should ONLY be APPENDED to, which
+will ensure that any older final data is not overwritten, thus ensuring
+historical reproducability.
+'''
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+from calendar import monthrange
+from datetime import datetime, timedelta
+
+# because datetime.strftime doesnt work before 1900
+def _yyyymmdd(t):
+  return('%04d%02d%02d' % (t.year, t.month, t.day))
+
+class Mg_SB_Daily:
+  ''' Daily time series for Solar Mg and SB indices from NRLSSI2 '''
+
+  def __init__(self, DATADIR,
+      filenames=[
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20170630_c20170928.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20191231_c20200225.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20201231_c20210204.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20220630_c20220803.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20230331_c20230411.txt',
+        'tsi-ssi_v02r01_model-input-time-series_s18820101_e20231231_c20240221.txt',
+        ],
+      verbose=False):
+
+    self.filenames = filenames
+    self.verbose   = verbose
+
+    if verbose:
+      print('reading Mg and SB indicies from', filenames)
+
+    # read the files
+    self._data = {}
+    self.nfinal = 0
+    first = True
+    for filename in filenames:
+      with open(os.sep.join((DATADIR,filename))) as f:
+        lines = f.readlines()
+
+      # process the data, ignoring the header (top line)
+      for line in lines[1:]:
+        datestr, SB, Mg, status = line.strip().split(',')
+        # yyyy-mm-dd to yyyymmdd
+        yyyymmdd = datestr.replace('-','')
+        # make data real
+        SB = float(SB)
+        Mg = float(Mg)
+        # validate status to rc
+        if status == 'final':
+          rc = 0
+        elif status == 'prelim':
+          rc = -1
+        else:
+          raise RuntimeError('invalid status detected: %s' % status)
+        # if a date is already finalized, ignore further records
+        #   for that date to ensure historical reproducability
+        if yyyymmdd in self._data:
+          rc_existing = self._data[yyyymmdd][0]
+          if rc_existing == 0: continue
+        # load data dictionary
+        self._data[yyyymmdd] = (rc, Mg, SB)
+        # keep track of min and max final times
+        if rc == 0:
+          self.nfinal += 1
+          d = datetime.strptime(yyyymmdd,'%Y%m%d').date()
+          if first:
+            self.date_min_final = d
+            self.date_max_final = d
+            first = False
+          else:
+            if d < self.date_min_final: self.date_min_final = d
+            if d > self.date_max_final: self.date_max_final = d
+
+  def keys(self):
+    return self._data.keys()
+
+  def getday(self, yyyymmdd):
+    '''
+    get daily Mg and SB values for daily string yyyymmdd
+    returns (rc, Mg, SB)
+    rc = 0 (final), -1 (prelim), -2 (unavailable)
+    '''
+    if yyyymmdd not in self._data:
+      return -2, np.nan, np.nan
+    else:
+      return self._data[yyyymmdd]
+
+  def gettime(self, t):
+    '''
+    get time-interpolated Mg and SB values for datetime t
+    returns (rc, Mg, SB)
+    rc = 0 (final), -1 (prelim), -2 (unavailable)
+    '''
+
+    # assume daily average valid at noon
+    tnoon = datetime(t.year,t.month,t.day,12)
+    vnoon = self.getday(_yyyymmdd(tnoon))
+    if t == tnoon: return vnoon
+
+    # other noon bracketing t
+    tother = tnoon + timedelta(days=(-1 if t < tnoon else +1))
+    vother = self.getday(_yyyymmdd(tother))
+
+    # fraction that the other daily average contributes
+    fother = abs((t-tnoon).total_seconds()) / 86400.
+
+    # only interpolate if both days exist
+    if vnoon[0] == -2 or vother[0] == -2:
+      return (-2, np.nan, np.nan)
+    else:
+      # now both days available
+      Mg = vnoon[1] * (1-fother) + vother[1] * fother
+      SB = vnoon[2] * (1-fother) + vother[2] * fother
+      rc = -1 if vnoon[0] == -1 or vother[0] == -1 else 0
+      return (rc, Mg, SB)
+
+  def _plot_all_final(self):
+    spyear = 365.25 * 86400.
+    d0 = self.date_min_final; d = d0
+    dy = []; Mg = []; SB = []
+    while (d <= self.date_max_final):
+      v = self.getday(_yyyymmdd(d))
+      if v[0] == 0:
+        # only plot finals
+        dy.append((d-d0).total_seconds()/spyear)
+        Mg.append(v[1])
+        SB.append(v[2])
+      d += timedelta(days=1)
+    plt.subplot(211); plt.plot(dy,Mg,'b-'); plt.ylabel('Mg')
+    plt.subplot(212); plt.plot(dy,SB,'b-'); plt.ylabel('SB')
+    plt.xlabel('years since %s' % self.date_min_final)
+
+if __name__ == '__main__':
+
+  MgSB = Mg_SB_Daily()
+# print('16000101', MgSB.getday('16000101'))
+# print('20161130', MgSB.getday('20161130'))
+# print('20161201', MgSB.getday('20161201'))
+# print('20161202', MgSB.getday('20161202'))
+# t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S')
+# print(t, MgSB.gettime(t))
+
+  plt.figure()
+  MgSB._plot_all_final()
+# plt.savefig(os.sep.join(('gx','MgSB_plot_all_final.png')),
+#   pad_inches=0.33,bbox_inches='tight',dpi=100)
+  plt.show()
+

GEOS_RadiationShared/NRLSSI2/README.md

@@ -0,0 +1,70 @@
+=========================================================
+Instructions on how to pre-process NRLSSI2 data for GEOS5
+=========================================================
+Peter Norris, April 2024
+
+A. Introduction.
+In the GEOS5 RRTMG and RRTMGP solar codes, the Solar Spectral Irradiance (SSI) at TOA is calculated from
+NRLSSI2 Total Solar Irradiance (TSI) and Mg (facular) and SB (sunspot) indices. This data is currently
+updated yearly, when the finalized yearly products become available, which is generally in late January
+to late February of the following year. For example, the finalized 2023 data came out on 2024-02-23. This
+README describes how to get the NRLSSI2 source data and how to pre-process it into a text file read by
+the SolarGC.
+
+B. Get the data for the NEW year.
+The source data is obtained from https://www.ncei.noaa.gov/data/total-solar-irradiance/access/
+Note: you can use wget https://... from discover.
+There are two data files required for each new year:
+1. from ancilliary-data/, e.g., tsi-ssi_v02r01_model-input-time-series_s18820101_e20180331_c20180418.txt
+  !!! NOTE: Please read comments in docstring of Mg_SB_from_daily_file.py BEFORE downloading !!!
+  In particular,
+    -- the <filenames> list in that python code should ONLY be APPENDED to,
+         which will ensure that any existing final data is not overwritten,
+         thus ensuring historical reproducibility.
+2. from daily/, e.g., tsi_v02r01_daily_s20170101_e20171231_c20180227.nc
+  !!! NOTE: Please read comments in docstring of TSI_from_daily_files.py BEFORE downloading !!!
+  In particular,
+    -- there should be no files with overlapping time periods, and
+    -- you should NEVER overwrite any existing non-preliminary file,
+         since this may cause historical non-reproducibility.
+NOTE: Only two files are needed for the new year. Existing data has already been obtained
+for previous years. There is no need to get historical data for past years. Even though the
+file in part 1 above contains historical data, only the data for the new year inside it will
+be used if you follow these instructions, in order to maintain historical reproducibility.
+
+C. Decide where to put the source data.
+The source data is put in the directory DATADIR, where DATADIR is currently
+   /discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar/NRLSSI2/data
+This directory is managed by the GMAO Operations Team, so if you is not they, you will need to download
+the data somewhere else first and then work with "Ops" to add it to that directory. DATADIR here is NOT
+an environmental variable, just a shorthand. But whatever DATADIR is used must also be explicity set in
+the "main" of the preprocessing driver: TSI_Mg_SB_merged_from_daily.py.
+
+D. Preparing the preprocessor.
+1. Add the filename from B.1. to the filenames list in Mg_SB_from_daily_file.py. Specifically APPEND
+it to the END of the <filenames> list in the default argument of the __init__ of class Mg_SB_Daily
+in Mg_SB_from_daily_file.py.
+2. Set the OUTDIR in the "main" of TSI_Mg_SB_merged_from_daily.py to where you want to store the 
+pre-processed file for the neaw year.
+
+E. Running the preprocessor.
+python TSI_Mg_SB_merged_from_daily.py
+After you exit the program, rename the NRLSSI2.VYYYY.txt in OUTDIR to the correct year.
+
+F. Storing the Output file
+Get "Ops" to put the new NRLSSI2.VYYYY.txt from your OUTDIR into their directory:
+  /discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar
+
+G. Symlinks
+There is a symlink in later "Ops solar" directory named NRLSSI2.txt.
+It will typically be pointed to the file NRLSSI2.vYYYY.txt for the new year, but that decision must
+be approved by Ops and by the Modelling Team and Data Assimilation Team heads.
+
+H. Using the data.
+The GEOS5 run should have the following lines in its ASGCM.rc:
+USE_NRLSSI2: .TRUE.
+SOLAR_CYCLE_FILE_NAME: ExtData/g5gcm/solar/NRLSSI2.txt
+(Note: These will only be used if the solar code being run is RRTMG or RRTMGP. Chou-Suarez runs do
+not use NRLSSI2 data).
+PS: You are free to change the use of the symlink NRLSSI2.txt to any specific NRLSSI2.vYYYY.txt 
+if you want reproducibility to a specific historical run.

GEOS_RadiationShared/NRLSSI2/TSI_Mg_SB_merged_from_daily.py

@@ -0,0 +1,159 @@
+'''
+Get TSI and Solar Mg and SB indices from NRLSSI2 daily data
+'''
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+from calendar import monthrange
+from datetime import datetime, timedelta
+from TSI_from_daily_files import TSI_Daily
+from Mg_SB_from_daily_file import Mg_SB_Daily
+
+# because datetime.strftime doesnt work before 1900
+def _yyyymmdd(t):
+  return('%04d%02d%02d' % (t.year, t.month, t.day))
+
+class TSI_Mg_SB_Merged_Daily:
+  ''' Daily time series for TSI and Solar Mg and SB indices from NRLSSI2 '''
+
+  def __init__(self, DATADIR, verbose=False):
+
+    self.verbose  = verbose
+
+    # get each data set
+    zM = Mg_SB_Daily(DATADIR, verbose=verbose)
+    zT =   TSI_Daily(DATADIR, verbose=verbose)
+
+    # form the INTERSECTION:
+    # (each set is indexed by yyyymmdd, unique in each set)
+    self._data = {}
+    self.nfinal = 0
+    first = True
+    for yyyymmdd in (zM.keys() & zT.keys()):
+
+        # load data dictionary
+        rc1, Mg, SB       = zM.getday(yyyymmdd)
+        rc2, TSI, TSI_UNC = zT.getday(yyyymmdd)
+        rc = -1 if rc1 == -1 or rc2 == -1 else 0
+        self._data[yyyymmdd] = (rc, TSI, Mg, SB, TSI_UNC)
+
+        # keep track of min and max final times
+        if rc == 0:
+          self.nfinal += 1
+          d = datetime.strptime(yyyymmdd,'%Y%m%d').date()
+          if first:
+            self.date_min_final = d
+            self.date_max_final = d
+            first = False
+          else:
+            if d < self.date_min_final: self.date_min_final = d
+            if d > self.date_max_final: self.date_max_final = d
+
+  def getday(self, yyyymmdd):
+    '''
+    get daily TSI, Mg and SB values for daily string yyyymmdd
+    returns (rc, TSI, Mg, SB, TSI_UNC)
+    rc = 0 (final), -1 (prelim), -2 (unavailable)
+    '''
+    if yyyymmdd not in self._data:
+      return -2, np.nan, np.nan, np.nan
+    else:
+      return self._data[yyyymmdd]
+
+  def gettime(self, t):
+    '''
+    get time-interpolated TSI, Mg, and SB values for datetime t
+    returns (rc, TSI, Mg, SB, TSI_UNC)
+    rc = 0 (final), -1 (prelim), -2 (unavailable)
+    '''
+
+    # assume daily average valid at noon
+    tnoon = datetime(t.year,t.month,t.day,12)
+    vnoon = self.getday(_yyyymmdd(tnoon))
+    if t == tnoon: return vnoon
+
+    # other noon bracketing t
+    tother = tnoon + timedelta(days=(-1 if t < tnoon else +1))
+    vother = self.getday(_yyyymmdd(tother))
+
+    # fraction that the other daily average contributes
+    fother = abs((t-tnoon).total_seconds()) / 86400.
+
+    # only interpolate if both days exist
+    if vnoon[0] == -2 or vother[0] == -2:
+      return (-2, np.nan, np.nan, np.nan)
+    else:
+      # now both days available
+      TSI     = vnoon[1] * (1-fother) + vother[1] * fother
+      Mg      = vnoon[2] * (1-fother) + vother[2] * fother
+      SB      = vnoon[3] * (1-fother) + vother[3] * fother
+      TSI_UNC = vnoon[4] * (1-fother) + vother[4] * fother
+      rc = -1 if vnoon[0] == -1 or vother[0] == -1 else 0
+      return (rc, TSI, Mg, SB, TSI_UNC)
+
+  def final_date_range(self,
+      yyyymmdd0, yyyymmdd1,
+      no_gaps=True):
+    ''' final only lists over date range '''
+
+    d0 = datetime.strptime(yyyymmdd0,'%Y%m%d').date()
+    d1 = datetime.strptime(yyyymmdd1,'%Y%m%d').date()
+    d = d0; dates = []; TSI = []; Mg = []; SB = []
+    while (d <= d1):
+      v = self.getday(_yyyymmdd(d))
+      if v[0] == 0:
+        # only include finals
+        dates.append(d)
+        TSI.append(v[1])
+        Mg.append( v[2])
+        SB.append( v[3])
+      elif no_gaps:
+        raise RuntimeError('no_gaps: intervening non-final date found')
+      d += timedelta(days=1)
+    return dates, TSI, Mg, SB
+
+  def _plot_all_final(self):
+    dates, TSI, Mg, SB = self.final_date_range(
+      _yyyymmdd(self.date_min_final),
+      _yyyymmdd(self.date_max_final),
+      no_gaps=False)
+    spyear = 365.25 * 86400.
+    dy = [(d-self.date_min_final).total_seconds()/spyear for d in dates]
+    plt.subplot(311); plt.plot(dy,TSI,'b-'); plt.ylabel('TSI')
+    plt.subplot(312); plt.plot(dy, Mg,'b-'); plt.ylabel('Mg')
+    plt.subplot(313); plt.plot(dy, SB,'b-'); plt.ylabel('SB')
+    plt.xlabel('years since %s' % self.date_min_final)
+
+  def output_final_textfile(self, filename):
+    f = open(filename,'w')
+    f.write('# NRLSSI2 daily input\n')
+    f.write('# treat daily values as valid at 12:00 GMT\n')
+    f.write('# yyyy doy TSI:W/m2 MgIndex   SBindex\n')
+    for d, TSI, Mg, SB in zip(*self.final_date_range(
+        _yyyymmdd(self.date_min_final), _yyyymmdd(self.date_max_final))):
+      f.write('  %04d %03d %8.3f %8.6f %9.4f\n' % (
+        d.year, d.timetuple().tm_yday, TSI, Mg, SB))
+    f.close()
+
+if __name__ == '__main__':
+
+  DATADIR = '/discover/nobackup/projects/gmao/share/gmao_ops/fvInput/g5gcm/solar/NRLSSI2/data'
+  OUTDIR = os.sep.join((os.environ['NOBACKUP'],'NRLSSI2','output'))
+
+  z = TSI_Mg_SB_Merged_Daily(DATADIR)
+# print('16000101', z.getday('16000101'))
+# print('20161130', z.getday('20161130'))
+# print('20161201', z.getday('20161201'))
+# print('20161202', z.getday('20161202'))
+# t = datetime.strptime('2016-12-01 10:00:00','%Y-%m-%d %H:%M:%S')
+# print(t, z.gettime(t))
+
+  z.output_final_textfile(os.sep.join((OUTDIR,'NRLSSI2.vYYYY.txt')))
+
+  plt.figure(figsize=(6,12))
+  z._plot_all_final()
+# plt.savefig(os.sep.join(('gx','TSInMgSB_plot_all_final.png')),
+#   pad_inches=0.33,bbox_inches='tight',dpi=100)
+  plt.show()
+