additional changes, mostly to the font end

pcmdi_metrics/cloud_feedback/cloud_feedback_driver.py

@@ -25,7 +25,6 @@
 version = param.version
 input_files_json = param.input_files_json
 path = param.path
-xml_path = param.xml_path
 data_path = param.data_path
 figure_path = param.figure_path
 output_path = param.output_path
@@ -45,7 +44,6 @@
 print("version:", version)
 print("path:", path)
 print("input_files_json:", input_files_json)
-print("xml_path:", xml_path)
 print("figure_path:", figure_path)
 print("output_path:", output_path)
 print("output_json_filename:", output_json_filename)
@@ -57,9 +55,6 @@
 else:
  exps = ["amip", "amip-p4K"]
 
-# generate xmls pointing to the cmorized netcdf files
-os.makedirs(xml_path, exist_ok=True)
-
 filenames = dict()
 
 if input_files_json is not None:
@@ -108,9 +103,7 @@
  searchstring = os.path.join(
  ncfiles[exp][field]["path"], ncfiles[exp][field]["file"]
  )
- xmlname = os.path.join(xml_path, ".".join([exp, model, variant, field, "xml"]))
- os.system("cdscan -x " + xmlname + " " + searchstring)
- filenames[exp][field] = xmlname
+ filenames[exp][field] = searchstring
 
 if debug:
  with open(os.path.join(output_path, "filenames.json"), "w") as f:
@@ -123,17 +116,17 @@
 
 # add this model's results to the pre-existing json file containing other models' results:
 updated_fbk_dict, updated_obsc_fbk_dict = organize_fbk_jsons(
- fbk_dict, obsc_fbk_dict, model, variant, datadir=data_path
+ fbk_dict, obsc_fbk_dict, model, variant
 )
-updated_err_dict = organize_err_jsons(err_dict, model, variant, datadir=data_path)
+updated_err_dict = organize_err_jsons(err_dict, model, variant)
 
 ecs = None
 if get_ecs:
  # calculate ECS and add it to the pre-existing json file containing other models' results:
  ecs = compute_ECS(filenames)
  print("calc ECS done")
  print("ecs: ", ecs)
-updated_ecs_dict = organize_ecs_jsons(ecs, model, variant, datadir=data_path)
+updated_ecs_dict = organize_ecs_jsons(ecs, model, variant)
 
 os.makedirs(output_path, exist_ok=True)
 if debug:

pcmdi_metrics/cloud_feedback/lib/__init__.py

@@ -1,6 +1,6 @@
 from .argparse_functions import AddParserArgument # noqa
-from .cal_CloudRadKernel import CloudRadKernel # noqa
-from .compute_ECS import compute_ECS # noqa
+from .cal_CloudRadKernel_xr import CloudRadKernel # noqa
+from .compute_ECS_xr import compute_ECS # noqa
 from .lib_cloud_feedback import cloud_feedback_metrics_to_json # noqa
 from .organize_jsons import ( # noqa
  organize_ecs_jsons,

pcmdi_metrics/cloud_feedback/lib/cal_CloudRadKernel_xr.py

@@ -19,7 +19,7 @@
 # define necessary information
 # =============================================
 
-datadir = "../data/"
+datadir = "./data/"
 
 # Define a python dictionary containing the sections of the histogram to consider
 # These are the same as in Zelinka et al, GRL, 2016
@@ -190,7 +190,7 @@ def get_CRK_data(filepath):
 ###########################################################################
 def get_kernel_regrid(ctl):
  # Read in data and map kernels to lat/lon
-
+ 
  f = xc.open_mfdataset(datadir + "cloud_kernels2.nc", decode_times=False)
  f = f.rename({"mo": "time", "tau_midpt": "tau", "p_midpt": "plev"})
  f["time"] = ctl["time"].copy()

pcmdi_metrics/cloud_feedback/lib/compute_ECS_xr.py

@@ -6,132 +6,118 @@
 from scipy import stats
 import cftime
 
-
 def dict_to_dataset(DICT):
  # Convert a dictionary of Datasets to a single Dataset
- data_vars = {}
+ data_vars={}
  for var in list(DICT.keys()):
  data_vars[var] = (["year"], DICT[var][var].data)
  ds = xr.Dataset(
  data_vars,
  coords=dict(year=DICT[var].year),
  )
- return ds
-
+ return(ds)
 
-def get_branch_time(pi, ab):
+def get_branch_time(pi,ab):
  # Determine when abrupt overlaps with piControl
  # returns the following:
  # 1) st: the branch time in cftime.datetime format
  # 2) fn: the last abrupt year in cftime.datetime format
  # 3) extensions beyond the ~150-year overlap period to allow for 21-year rolling mean:
  # 3a) extend_st -- the date 10 years prior to st
  # 3b) extend_fn -- the date 10 years after fn
-
+ 
  var = ab.variable_id
  source = ab.source_id
  experiment = ab.experiment_id
- absub = ab.isel(time=slice(0, 12 * 150))
- lenab = int(absub[var].shape[0] / 12) - 1
- btip = int(absub.attrs["branch_time_in_parent"])  # timestep, not a date
+ absub = ab.isel(time=slice(0,12*150))
+ lenab = int(absub[var].shape[0]/12)-1
+ btip = int(absub.attrs['branch_time_in_parent']) # timestep, not a date
 
- if source == "GFDL-CM4" and (
- experiment == "1pctCO2"
- or experiment == "abrupt-4xCO2"
- or experiment == "historical"
- ):
- btip = 91250 # https://errata.es-doc.org/static/view.html?uid=2f6b5963-f87e-b2df-a5b0-2f12b6b68d32
- ptu = absub.attrs["parent_time_units"]
+ if source=='GFDL-CM4' and (experiment=='1pctCO2' or experiment=='abrupt-4xCO2' or experiment=='historical'): 
+ btip = 91250 # https://errata.es-doc.org/static/view.html?uid=2f6b5963-f87e-b2df-a5b0-2f12b6b68d32 
+ ptu = absub.attrs['parent_time_units']
  # print(' branch time in parent: '+str(btip)+' '+ptu)
  # print(' parent times: '+str(pi.time.dt.year[0].values)+' - '+str(pi.time.dt.year[-1].values))
- start = int(ptu.split(" ")[2][:4])
+ start = int(ptu.split(' ')[2][:4])
  cal = pi.time.dt.calendar
  # see http:https://cfconventions.org/cf-conventions/cf-conventions.html#calendar
- if cal == "360_day":
- dpy = 360
- elif cal == "noleap":
- dpy = 365
- elif cal == "proleptic_gregorian" or cal == "standard" or cal == "julian":
- dpy = 365.25
+ if cal=='360_day':
+ dpy=360
+ elif cal=='noleap':
+ dpy=365
+ elif cal=='proleptic_gregorian' or cal=='standard' or cal=='julian':
+ dpy=365.25
  else:
- raise (Exception("Not sure how to handle " + cal + " calendar"))
- if source == "CIESM" and (experiment == "abrupt-4xCO2" or experiment == "1pctCO2"):
+ raise(Exception('Not sure how to handle '+cal+' calendar'))
+ if source=='CIESM' and (experiment=='abrupt-4xCO2' or experiment=='1pctCO2'):
  # info from Yi Qin 10/1/20:
- # "The branch_time for 4xCO2 should be the 101-yr of these published 500-yr piControl data,
+ # "The branch_time for 4xCO2 should be the 101-yr of these published 500-yr piControl data, 
  # but I wrongly regarded the “raw” piControl with piControl-spinup as the parent_case, which is a 1000-yr long run.
- year0 = 101
- elif source == "KACE-1-0-G" and (
- experiment == "abrupt-4xCO2" or experiment == "1pctCO2"
- ):
+ year0 = 101 
+ elif source=='KACE-1-0-G' and (experiment=='abrupt-4xCO2' or experiment=='1pctCO2'):
  # info from Jae-Hee Lee (via Gavin Schmidt on 1/26/22):
- # "We used year 2150 of piControl as the initial condition for the KACE-1-0-G 1pctCO2 runs.
+ # "We used year 2150 of piControl as the initial condition for the KACE-1-0-G 1pctCO2 runs. 
  # so, you can assume that 1850 in the 1pctCO2 and 2150 in the piControl are the same
- year0 = 2150
+ year0 = 2150 
  else:
- year0 = int((btip / dpy) + start)
- year150 = int(year0 + lenab)
- st = cftime.datetime(year0, 1, 1).strftime("%Y-%m-%dT%H:%M:%S")
- fn = cftime.datetime(year150, 12, 31, 23, 59, 59).strftime("%Y-%m-%dT%H:%M:%S")
+ year0 = int((btip/dpy)+start)
+ year150 = int(year0+lenab)
+ st=cftime.datetime(year0, 1, 1).strftime("%Y-%m-%dT%H:%M:%S")
+ fn=cftime.datetime(year150, 12, 31,23,59,59).strftime("%Y-%m-%dT%H:%M:%S")
  # Add on 10 years before st and after fn to assist in computing 21-year rolling means
- extend_st = cftime.datetime(np.max((year0 - 10, 1)), 1, 1).strftime(
- "%Y-%m-%dT%H:%M:%S"
- )
- extend_fn = cftime.datetime(
- year150 + 10, 12, pi.time.dt.day[-1], 23, 59, 59
- ).strftime("%Y-%m-%dT%H:%M:%S")
+ extend_st=cftime.datetime(np.max((year0-10,1)), 1, 1).strftime("%Y-%m-%dT%H:%M:%S")
+ extend_fn=cftime.datetime(year150+10, 12, pi.time.dt.day[-1],23,59,59).strftime("%Y-%m-%dT%H:%M:%S") 
 
- return (extend_st, extend_fn, st, fn) # ,year0,year150)
+ return(extend_st,extend_fn,st,fn)#,year0,year150)
 
-
-def get_data(pi, ab):
+def get_data(pi,ab): 
  # Return the piControl running mean and the abrupt deviation from this
 
- var = ab.variable_id
- extend_st, extend_fn, st, fn = get_branch_time(pi, ab)
- absub = ab.isel(time=slice(0, 12 * 150))
- pisub = pi.sel(time=slice(extend_st, extend_fn))
- if len(pisub.time) < 1200:
- print("len(time)<1200...skipping")
+ var=ab.variable_id
+ extend_st,extend_fn,st,fn = get_branch_time(pi,ab)
+ absub = ab.isel(time=slice(0,12*150))
+ pisub = pi.sel(time=slice(extend_st,extend_fn))
+ if len(pisub.time)<1200:
+ print('len(time)<1200...skipping')
  return None
 
- bdate = pisub.time.dt.year.values
+ bdate = pisub.time.dt.year.values 
  # compute global means:
  GLcntl = pisub.spatial.average(var)
  GLpert = absub.spatial.average(var)
  # compute annual means:
- annpert = GLpert.groupby("time.year").mean("time")
- anncntl = GLcntl.groupby("time.year").mean("time")
+ annpert = GLpert.groupby('time.year').mean('time')
+ anncntl = GLcntl.groupby('time.year').mean('time')
  # compute 21-year centered rolling mean:
- runcntl = anncntl.rolling(year=21, min_periods=1, center=True).mean()
+ runcntl = anncntl.rolling(year = 21, min_periods = 1, center = True).mean()
  # subset the rolling mean for the actual overlapping time:
- runcntlsub = runcntl.sel(year=slice(int(st[:4]), int(fn[:4])))
+ runcntlsub = runcntl.sel(year=slice(int(st[:4]),int(fn[:4])))
  # set the abrupt run's year coord to match the overlapping piCon's year coord
- annpert = annpert.assign_coords(year=runcntlsub.coords["year"])
+ annpert=annpert.assign_coords(year=runcntlsub.coords['year'])
  # compute abrupt anomalies from piControl running mean:
  annanom = annpert - runcntlsub
+
+ return(runcntl,annanom)
 
- return (runcntl, annanom)
-
-
-def compute_abrupt_anoms(pifilepath, abfilepath):
+def compute_ECS(filepath):
  # compute annual- anad global-mean tas and EEI anomalies in abrupt4xCO2 w.r.t. piControl climo, then compute ERF, lambda, and ECS via Gregory
- cntl = {}
- anom = {}
- skip = False
- variables = list(pifilepath.keys())
+ cntl={}
+ anom={}
+ skip=False
+ variables = list(filepath["piControl"].keys())
  for var in variables:
- pi = xcdat.open_mfdataset(pifilepath[var], use_cftime=True)
- ab = xcdat.open_mfdataset(abfilepath[var], use_cftime=True)
-
- if len(ab.time) < 140 * 12:
- print(" len(ab.time)<140*12")
- skip = True
+ pi=xcdat.open_mfdataset(filepath["piControl"][var], use_cftime = True)
+ ab=xcdat.open_mfdataset(filepath["abrupt-4xCO2"][var], use_cftime = True)
+ 
+ if len(ab.time)<140*12:
+ print(' len(ab.time)<140*12')
+ skip=True
  break
- output = get_data(pi, ab)
+ output = get_data(pi,ab)
  if output:
- cntl[var], anom[var] = output
+ cntl[var],anom[var] = output
  else:
- skip = True
+ skip=True
  break
  if skip:
  # print('Skipping this model...')
@@ -140,18 +126,16 @@ def compute_abrupt_anoms(pifilepath, abfilepath):
  CNTL = dict_to_dataset(cntl)
  ANOM = dict_to_dataset(anom)
 
- CNTL["eei"] = CNTL["rsdt"] - CNTL["rsut"] - CNTL["rlut"]
- ANOM["eei"] = ANOM["rsdt"] - ANOM["rsut"] - ANOM["rlut"]
-
- x = ANOM["tas"].load()
- y = ANOM["eei"].load()
+ CNTL['eei'] = CNTL['rsdt']-CNTL['rsut']-CNTL['rlut']
+ ANOM['eei'] = ANOM['rsdt']-ANOM['rsut']-ANOM['rlut']
 
- return (x, y)
-
-
-def gregory_calcs(x, y):
- m, b, r, p, std_err = stats.linregress(x, y)
- ERF = b / 2
+ x = ANOM['tas'].load()
+ y = ANOM['eei'].load()
+
+ m, b, r, p, std_err = stats.linregress(x,y)
+ ERF = b/2
  lam = m
- ECS = -ERF / lam
- return (ERF, lam, ECS)
+ ecs = -ERF/lam
+
+ return ecs
+

pcmdi_metrics/cloud_feedback/param/my_param.py

@@ -12,7 +12,6 @@
 get_ecs = False
 
 # Output directory path (directory will be generated if it does not exist yet.)
-xml_path = "./xmls/"
 figure_path = "./figures/"
 output_path = "./output"
 output_json_filename = "_".join(["cloud_feedback", model, variant]) + ".json"