new functions added:

- remove_incomplete_timeseries(amoc_index): removes members that are incomplete
- get_amoc_strengths(amoc_index): returns dataset with AMOC strength at 26.5N and at 35N

get_subpolar_gyre_functions.py

@@ -10,6 +10,7 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable
 import cartopy.crs as ccrs
 import os
 import glob
+import json
 import xarray as xr
 
 
@@ -111,6 +112,9 @@ def get_subpolar_gyre(model, emember, print_info=False, paths_all_models=None):
         return None
     else:
         data = xr.open_mfdataset(paths[0], concat_dim="time").tos
+        if model == 'CNRM-CM6-1-HR':
+            # this crashes with a oom-kill, so re-chunking with 2-timestep-chunks seems to help...
+            data = data.chunk(dict(time=2))
 
         str_lon = get_lon_lat_keys(data, lon_keys)
         str_lat = get_lon_lat_keys(data, lat_keys)
@@ -127,6 +131,72 @@ def get_subpolar_gyre(model, emember, print_info=False, paths_all_models=None):
         return spg
 
 
+# set everything to nan where there are no full timeseries for the AMOC strength at 26.5N and 35N
+def remove_incomplete_timeseries(amoc_index):
+    amoc = get_amoc_strengths(amoc_index)
+    amoc_strength_265 = amoc.amoc_strength_265.astype(np.float)
+    amoc_strength_35  = amoc.amoc_strength_35.astype(np.float)
+    
+    incomplete = (np.isnan(amoc_strength_265).sum(dim="year") > 0)  & (np.isnan(amoc_strength_35).sum(dim="year") > 0 )
+
+    for m in range(len(amoc_strength_265.models.values)):
+        for em in range(len(amoc_strength_265.ensemble_members.values)):
+            if incomplete.isel(models=m, ensemble_members=em):
+                amoc_strength_265.isel(models=m, ensemble_members=em).values.fill(np.nan)
+                amoc_strength_35.isel(models=m, ensemble_members=em).values.fill(np.nan)
+                amoc_index.isel(models=m, ensemble_members=em).values.fill(np.nan)
+
+    # remove the UKESM1-0-LL r1i1p1f1 in the amoc-strength data
+    model   = 'UKESM1-0-LL'
+    emember = 'r1i1p1f1'
+    amoc_strength_265.sel(models = model, ensemble_members = emember).values.fill(np.nan)
+    amoc_strength_35.sel(models = model, ensemble_members = emember).values.fill(np.nan)
+
+    # save all CMIP6 amoc data in one file
+    xr.Dataset(dict(
+        strength_265N = amoc_strength_265,
+        strength_35N  = amoc_strength_35,
+        index         = amoc_index.groupby('time.year').mean(dim="time")
+    )).to_netcdf("CMIP6_amoc.nc", mode='w')
+    return amoc_index
+
+
+# getMatthew's AMOC strength
+def get_amoc_strengths(amoc_index):
+    with open('matthew/JSON_data/Figure_AR6_CMIP5-6_AMOC_35N_1000m.json', 'r') as handle:
+        json_load = json.load(handle)
+
+    amoc_c6_ts = np.ma.asarray(json_load["amoc_c6_ts"])
+    cmip6_models = json_load["cmip6_models"]
+    year = np.asarray(json_load["year"])
+    experiments = ['ssp119', 'ssp126', 'ssp245', 'ssp370', 'ssp585']
+    latitudes   = ['26.5N', '35N']
+    ''' 
+    Ensemble members:
+    These are the first 10 ensemble members r${ens_num}i1p1f1 in the respective experiments. 
+    Where "f1" was not available we have used "f2" or "f3" and so on.
+    '''
+    amoc_strength = xr.DataArray(
+        amoc_c6_ts,
+        dims = ('models', 'experiments', 'ensemble_members', 'latitudes', 'year'),
+        coords = dict(
+            models      = xr.DataArray(cmip6_models, dims="models", coords=dict(models=("models", cmip6_models))),
+            experiments = xr.DataArray(experiments, dims="experiments", coords=dict(experiments=("experiments", experiments))),
+            latitudes   = xr.DataArray(latitudes, dims="latitudes", coords=dict(latitudes=("latitudes", latitudes))),
+            year        = xr.DataArray(year, dims="year", coords=dict(year=("year", year))),
+            ensemble_members = amoc_index.ensemble_members
+            )
+        )
+    amoc_index = amoc_index.groupby('time.year').mean(dim="time")
+
+    # the historical period is same over all experiments, so just choose any
+    return xr.Dataset(
+        dict(
+            amoc_strength_265 = amoc_strength.sel(latitudes='26.5N', experiments='ssp119', drop=True).sel(year=amoc_index.year),
+            amoc_strength_35  = amoc_strength.sel(latitudes='35N',   experiments='ssp119', drop=True).sel(year=amoc_index.year)
+        )
+    )
+
 
 ## CREATE NEW OUTPUT FILE
 def create_new_amoc_index_file(output_file):