diff --git a/FACTS.py b/FACTS.py index d46c9f98..18e7cb64 100644 --- a/FACTS.py +++ b/FACTS.py @@ -395,7 +395,7 @@ def IdentifyClimateOutputFiles(pcfg,pipe_name): tcfg = pcfg[this_stage][this_task] if "climate_output_data" in tcfg.keys(): for this_file in tcfg['climate_output_data']: - if this_file.__contains__('climate.nc'): + if (this_file.__contains__('climate.nc') or this_file.__contains__('climate.h5')): pd['climate'] = '$SHARED/climate/' + mvar_replace_dict(mvar_dict, this_file) elif this_file.__contains__('gsat.nc'): pd['gsat'] = '$SHARED/climate/' + mvar_replace_dict(mvar_dict, this_file) diff --git a/experiments/ebm3.global/config.yml b/experiments/ebm3.global/config.yml index 044f60fa..9670377e 100644 --- a/experiments/ebm3.global/config.yml +++ b/experiments/ebm3.global/config.yml @@ -11,11 +11,7 @@ climate_step: fair2: module_set: "fair2" module: "climate" - pipeline_file: 'pipeline.yml' - # input_data_file: - # - "ebm3.global.fair2.fair2.climate_ohc.nc" - # - "ebm3.global.fair2.fair2.climate_gsat.nc" - # - "ebm3.global.fair2.fair2.climate_climate.nc" + pipeline_file: 'pipeline.global.yml' climate_output_data: - "ebm3.global.fair2.fair2.climate_ohc.nc" - "ebm3.global.fair2.fair2.climate_gsat.nc" @@ -24,6 +20,6 @@ climate_step: generates_climate_output: true sealevel_step: - ebm3: + ocean: module_set: "ebm3" - module: "thermalexpansion" \ No newline at end of file + module: "sterodynamics" \ No newline at end of file diff --git a/experiments/ebm3.local/config.yml b/experiments/ebm3.local/config.yml new file mode 100644 index 00000000..23d68c88 --- /dev/null +++ b/experiments/ebm3.local/config.yml @@ -0,0 +1,25 @@ +global-options: + nsamps: 50 + scenario: ssp585 + pyear_start: 2020 + pyear_end: 2150 + pyear_step: 10 + baseyear: 2005 + pipeline_file: "pipeline.yml" + +climate_step: + fair2: + module_set: "fair2" + module: "climate" + pipeline_file: 'pipeline.yml' + climate_output_data: + - "ebm3.local.fair2.fair2.climate_ohc.nc" + - "ebm3.local.fair2.fair2.climate_gsat.nc" + - "ebm3.local.fair2.fair2.climate_climate.nc" + - "ebm3.local.fair2.fair2.climate_oceantemp.nc" + generates_climate_output: true + +sealevel_step: + ocean: + module_set: "ebm3" + module: "sterodynamics" \ No newline at end of file diff --git a/experiments/emulandice2.protect.ssp585/config.yml b/experiments/emulandice2.protect.ssp585/config.yml new file mode 100644 index 00000000..1f1c110f --- /dev/null +++ b/experiments/emulandice2.protect.ssp585/config.yml @@ -0,0 +1,26 @@ +global-options: + nsamps: 2237 + scenario: ssp585 + pyear_start: 2020 + pyear_end: 2300 #2150 + pyear_step: 10 + baseyear: 2005 + pipeline_file: "pipeline.2300.yml" # "pipeline.yml" + +climate_step: + temperature: + module_set: "fair" + module: "temperature" + pipeline_file: "pipeline.global.yml" + generates_climate_output: true + +sealevel_step: + emuAIS: + module_set: "emulandice2" + module: "AIS" + emuGLA: + module_set: "emulandice2" + module: "glaciers" + emuGIS: + module_set: "emulandice2" + module: "GrIS" diff --git a/modules-data/modules-data.global_only.urls.txt b/modules-data/modules-data.global_only.urls.txt index 48956186..58ee43b6 100644 --- a/modules-data/modules-data.global_only.urls.txt +++ b/modules-data/modules-data.global_only.urls.txt @@ -10,4 +10,4 @@ https://zenodo.org/record/7478192/files/larmip_icesheet_project_data.tgz https://zenodo.org/record/7478192/files/ssp_landwaterstorage_preprocess_data.tgz https://zenodo.org/record/7478192/files/tlm_sterodynamics_preprocess_data.tgz https://zenodo.org/records/11506798/files/fair2_climate_project_data.tgz -https://zenodo.org/records/11506798/files/ebm3_thermal_expansion_data.tgz \ No newline at end of file +https://zenodo.org/records/16024082/files/ebm3_project_data.tgz \ No newline at end of file diff --git a/modules-data/modules-data.urls.txt b/modules-data/modules-data.urls.txt index 38667d38..ea2ecb56 100644 --- a/modules-data/modules-data.urls.txt +++ b/modules-data/modules-data.urls.txt @@ -30,3 +30,5 @@ https://zenodo.org/records/11506798/files/fair2_climate_project_data.tgz https://zenodo.org/records/11506798/files/ebm3_thermal_expansion_data.tgz https://zenodo.org/records/18199757/files/oelsmann24_vlm_data.tar.gz https://zenodo.org/records/15975768/files/caron18_GIA_preprocess_data.tgz +https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_ebm3_data.tgz +https://zenodo.org/records/16024082/files/ebm3_project_data.tgz \ No newline at end of file diff --git a/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py b/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py new file mode 100644 index 00000000..e966d8c3 --- /dev/null +++ b/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py @@ -0,0 +1,674 @@ +# Sterodynamics projections for PROTECT +# Victor Malagon Santos, 17.06.2024 + +import numpy as np +import os +import pandas as pd +import sys +import argparse +import energy_balance_model as ebm3 +import csv +import xarray as xr +from sklearn import linear_model +import re +import time +import netCDF4 +import h5py +import scipy +import pickle + +""" +Created on Fri Jun 7 11:24:19 2024 + +For any model you have up tp 2300, you want a set of scaling coefficients +for ssp126 and ssp585, as those thwo a typically the ones available up to 2300. +Those scenarios are used to get a pattern of regression coefficient used for +SSP interpolation and probabilistic projections using FAIR. + +@author: vmalagonsantos +""" + +# functions +def fetch_erfs_from_rcmip(path,scenarios): + ''' + function to read scenario-dependent ERF from RCMIP csv + + Input parameters: + path = path to RCMIP csv + scenarios = desired ssp's + + Output: Adds results of extreme sea-level analysis for all stations to the + station data dictionary. + + erfs = effective radiative forcing per scenario and year + erfyears = corresponding years (1750-2500) + ''' + + erfyears=np.arange(1750,2501) + ssp_idx = {'ssp119':212,'ssp126':231,'ssp245':308,'ssp370':59,'ssp585':404} #location in table + + with open(path) as csv_file: + csv_reader = csv.reader(csv_file) + rows = list(csv_reader) + + erfs = np.empty((len(scenarios),len(erfyears))) + + for s,scen in enumerate(scenarios): + try: + erfs[s,:] = rows[ssp_idx[scen]][7:] + except: + continue + + return erfs, erfyears + +def angd(lat0, lon0, qlat, qlon): + + # Convert the input from degrees to radians + (lat0, lon0) = np.radians((lat0, lon0)) + (qlat, qlon) = np.radians((qlat, qlon)) + + # Calculate the angle between the vectors + temp = np.arctan2(np.sqrt((np.cos(qlat)*np.sin(qlon-lon0))**2 + \ + (np.cos(lat0)*np.sin(qlat) - np.sin(lat0)*np.cos(qlat) * np.cos(qlon-lon0))**2),\ + (np.sin(lat0)*np.sin(qlat) + np.cos(lat0)*np.cos(qlat)*np.cos(qlon-lon0))) + + # Convert the results from radians to degrees and return + return(np.degrees(temp)) + +def NearestPoint(qlat, qlon, lats, lons, tol = None): + + # Get the distance between the query point and all the possible points + dist = angd(lats, lons, qlat, qlon) + + # Which is the closest point + nearest_idx = np.argmin(dist) + + # Is the point within the tolerance? + if isinstance(tol, (int, float)): + if dist[nearest_idx] > tol: + return(None) + + return(nearest_idx) + + +def NearestPoints(qlats, qlons, lats, lons, tol): + + if len(qlats) != len(qlons): + raise Exception("Query lats ({}) and lons ({}) differ in length".format(len(qlats), len(qlons))) + + idx = map(lambda qlat,qlon: NearestPoint(qlat, qlon, lats, lons, tol), qlats, qlons) + + return(list(idx)) + +def ReadLocationFile(location_file): + + # Initialize variables to hold data and site information + names = [] + ids = [] + lats = [] + lons = [] + + # Compile the regex for finding commented lines + comment_regex = re.compile(r'^#') + + # Open the rate file + with open(location_file, 'r') as f: + + # Loop over the lines of the file + for line in f: + + # Skip commented lines + if re.search(comment_regex, line): + continue + + # Split the line into components + (this_name, this_id, this_lat, this_lon) = line.split("\t") + + # Store the information + names.append(this_name) + ids.append(int(this_id)) + lats.append(float(this_lat)) + lons.append(float(this_lon)) + + # Cast everything as numpy arrays + names = np.array(names) + ids = np.array(ids) + lats = np.array(lats) + lons = np.array(lons) + + # Return variables + return(names, ids, lats, lons) + + +# Smooth ZOS and ZOSTOGA over 19 year smoothing window +def nanSmooth(x, w=19): + idx = np.flatnonzero(~np.isnan(x)) + temp = x + if len(idx) > 0: + temp[idx] = Smooth(x[idx], w) + return(temp) + +def Smooth(x, w=19): + out0 = np.convolve(x, np.ones(w,dtype='double'), 'valid')/w + r = np.arange(1,w-1,2, dtype="double") + start = np.cumsum(x[:w-1])[::2]/r + stop = (np.cumsum(x[:-w:-1])[::2]/r)[::-1] + y = np.concatenate((start, out0, stop)) + return(y) + +def emb3_thermalexpansion_postprocess(scenario, pipeline_id, nsamps, seed, pyear_start, pyear_end, pyear_step, locationfile, baseyear, climate_data_file, rfmip, params, zosdir): + if pyear_end < 2151: + ryear_end = 2100 + elif pyear_end > 2300: + ryear_end = 2300 + else: + ryear_end = pyear_end + + targyears = np.arange(pyear_start, ryear_end) # for regression + projyears = np.arange(pyear_start, pyear_end+1, pyear_step) + + (_, site_ids, site_lats, site_lons) = ReadLocationFile(locationfile) + + # get temperature from FaIR simulations: INPUT FROM CLIMATE STEP. We need both gsat and oceantemp + cf = xr.open_dataset(climate_data_file, group=scenario, engine='netcdf4', ) + gsat = cf['surface_temperature'].sel(years=projyears) - cf['surface_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + otemp = cf['deep_ocean_temperature'].sel(years=projyears) - cf['deep_ocean_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + # INPUT temperature file from ebm3 global + gte_file = f'{pipeline_id}_globalsl.nc' + + # INPUT get models and parameters + paramdir = params + + #forcing for ebm + scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections + rfmipfile = rfmip #path # RFMIP FILE IS INPUT, NEW DATA + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries + + ebm_param = pd.read_csv(paramdir) + + # find available zos simulations - ignoring hidden files + model_has_zos = [f for f in os.listdir(zosdir) if not f.startswith('.')] + + ## FITTING + slopes = [] + intercept = [] + + # if projections are below 2150, look for 2100 zos simulations available for requested SSP, run ebm adn estimate slopes + if pyear_end < 2151: # # keep those up to 2300 or 2500 + for m, model in enumerate(model_has_zos): + ssp1_file = [] + hist_file = [] + # model = 'ACCESS-ESM1-5' + + # find matching ebm parameters to model + # df['company_name'].eq('ABC').any() + model_param = ebm_param[ebm_param["model"].eq(model)] + modeldir = zosdir + model + + zos_runs2_2100 = [f for f in os.listdir(modeldir) if scenario in f] # regex would be better here + if not zos_runs2_2100: + print(' ') + print(scenario + ' scenario not found in model ' + model + '. Moving on.') + continue + print(' ') + print('Estimating regression parameters for model: ' + model) + + + # get variants for 2100 simulation + + for v, filename in enumerate(zos_runs2_2100): + variant = [f for f in re.split('_', filename) if f.startswith('r')][0] + # get EBM parameters for both model and variant + model_var_param = model_param[model_param["run"].str.contains(variant)] + if not model_var_param.empty: #if not empty, check if that variant is availble for both hist and ssps + # get zos, historical + ssp + hist_file = [f for f in os.listdir(modeldir) if 'hist' in f and variant in f] + ssp_file = [f for f in os.listdir(modeldir) if scenario in f and variant in f] + if (hist_file and ssp_file): # if there is zos sims for hist and ssps and ebm parameters, continue with emulation + break + if model_var_param.empty: + print('No EBM parameters found for model ' + model) + continue + + + # get model parameters + c1 = np.array(model_var_param.get('C1')) + c2 = np.array(model_var_param.get('C2')) + c3 = np.array(model_var_param.get('C3')) + k1 = np.array(model_var_param.get('kappa1')) + k2 = np.array(model_var_param.get('kappa2')) + k3 = np.array(model_var_param.get('kappa3')) + e = np.array(model_var_param.get('epsilon')) + f4 = np.array(model_var_param.get('F_4xCO2')) + se = np.array(model_var_param.get('sigma_eta')) + sx = np.array(model_var_param.get('sigma_xi')) + gamma = np.array(model_var_param.get('gamma')) + nit = np.array(model_var_param.get('nit')) + + + # get forcing + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, [scenario]) #get ERF timeseries + + + #Model setup 3lm + ebm_obj = ebm3.EnergyBalanceModel( + ocean_heat_capacity = [c1[0], c2[0], c3[0]], + ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value + deep_ocean_efficacy = e[0], + forcing_4co2 = f4[0], + stochastic_run=False, + sigma_eta = se[0], + sigma_xi = sx[0], + gamma_autocorrelation = gamma[0], + seed = None, + timestep = 1, + n_timesteps = 2500-1750, + ) + + # Run model for different scenarios + ebm_obj.add_forcing( erfs[0,:], erfyears) + ebm_obj.run() + temp3 = ebm_obj.temperature + + attrs={'description': 'surface and deep temperature obtained by EBM-3LM', + 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} + + temp3_xr = xr.Dataset({'temperature':(['years','layers'], temp3)}, + coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears}, attrs=attrs) + + # reference to base year + temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) + # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + try: + hist_xr = xr.open_dataset(modeldir + '/' + hist_file[0],decode_times=False) + ssp_xr = xr.open_dataset(modeldir + '/' + ssp_file[0],decode_times=False) + except: + print('Model misses either historical or ssp simulation. Moving on') + continue + + if ssp_xr['time'][0].values != 60280.5: + ssp_xr.coords['time'] = ((ssp_xr['time'] + 60280.5).astype('float')) + + # cocatenate + zos = xr.concat([hist_xr, ssp_xr], dim='time') + + #% create new coordinates for year + zos.coords['years'] = ((zos['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + + # Calculate annual means + zos_am = zos.groupby('years').mean(dim='time') + zos_am = zos_am.where(zos_am['zos'] < 99999) # replace land values with nans + + # reference zos + # baseyear_idx = np.flatnonzero(datayears == baseyear) + # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) + try: + zos_am = zos_am.sel(years=targyears) - zos_am.sel(years=baseyear) + except: + print('Model has inconsistent time. Moving on') + continue + + #% fitting stage + nlon = len(zos_am['lon']) + nlat = len(zos_am['lat']) + + # get temperature for desired scenario + Ts = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='surface')) + Ti = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='inter')) + Td = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='deep')) + + + zos = np.array(zos_am['zos'].sel(years=targyears)) + + slope = np.zeros((3, nlat, nlon)) + slope[:] = np.nan + smoothwin = 19 + + reg = linear_model.LinearRegression() + # for i in tqdm(range(nlat)): + # sleep(3) + for i in range(nlat): + for j in range(nlon): + + y = np.array(Smooth(zos[:,i,j],w=smoothwin)) # reducing varability + x1 = Ts.flatten() + x2 = Ti.flatten() + x3 = Td.flatten() + x = np.transpose((x1,x2,x3)) + + try: + + reg.fit(x, y) + slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface + slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter + slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep + # intercept[m,i,j] = reg.intercept_ + + except ValueError: + slope[0,i,j] = np.nan + slope[1,i,j] = np.nan + slope[2,i,j] = np.nan + # intercept[m,i,j] = np.nan + + slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? + + + # if projections are above 2150, look for 2300 zos simulations available for SSP126 and SSP585, run ebm, estimate slopes, and emulate requested scenario + if pyear_end >= 2151: + for m, model in enumerate(model_has_zos): + + # find matching ebm parameters to model + model_param = ebm_param[ebm_param["model"].str.contains(model)] + modeldir = zosdir + model + + # keep those up to 2300 or 2500 + zos_runs2_2300 = [f for f in os.listdir(modeldir) if '2300' in f or '2500' in f] # regex would be better here + if not zos_runs2_2300: + print(' ') + print('Model ' + model + ' only runs to 2100. Skipping and looking for longer simulations' ) + continue + print(' ') + print('Estimating parameters for model: ' + model) + + for v, filename in enumerate(zos_runs2_2300): + ssp1_file = [] + ssp5_file = [] + hist_file = [] + + variant = [f for f in re.split('_', filename) if f.startswith('r')][0] + + # get EBM parameters for both model and variant + model_var_param = model_param[model_param["run"].str.contains(variant)] + if not model_var_param.empty: #if not empty, check if that variant is availble for both hist and ssps + # get zos, historical + ssp + hist_file = [f for f in os.listdir(modeldir) if 'hist' in f and variant in f] + ssp1_file = [f for f in os.listdir(modeldir) if 'ssp126' in f and variant in f] + ssp5_file = [f for f in os.listdir(modeldir) if 'ssp585' in f and variant in f] + if (hist_file and ssp1_file and ssp5_file): # if there is zos sims for hist and ssps and ebm parameters, continue with emulation + break + + + # define EBM parameters + c1 = np.array(model_var_param.get('C1')) + c2 = np.array(model_var_param.get('C2')) + c3 = np.array(model_var_param.get('C3')) + k1 = np.array(model_var_param.get('kappa1')) + k2 = np.array(model_var_param.get('kappa2')) + k3 = np.array(model_var_param.get('kappa3')) + e = np.array(model_var_param.get('epsilon')) + f4 = np.array(model_var_param.get('F_4xCO2')) + se = np.array(model_var_param.get('sigma_eta')) + sx = np.array(model_var_param.get('sigma_xi')) + gamma = np.array(model_var_param.get('gamma')) + nit = np.array(model_var_param.get('nit')) + + scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries + + + temp3 = np.empty((len(erfyears), 3, len(scenarios))) + + for s, scen in enumerate(scenarios): + + #Model setup 3lm + ebm_obj = ebm3.EnergyBalanceModel( + ocean_heat_capacity = [c1[0], c2[0], c3[0]], + ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value + deep_ocean_efficacy = e[0], + forcing_4co2 = f4[0], + stochastic_run=False, + sigma_eta = se[0], + sigma_xi = sx[0], + gamma_autocorrelation = gamma[0], + seed = None, + timestep = 1, + n_timesteps = 2500-1750, + ) + + # Run model for different scenarios + ebm_obj.add_forcing( erfs[s,:], erfyears) + ebm_obj.run() + temp3[:,:,s] = ebm_obj.temperature + + attrs={'description': 'surface and deep temperature obtained by EBM-3LM', + 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} + + temp3_xr = xr.Dataset({'temperature':(['years','layers', 'scenarios'], temp3)}, + coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears, 'scenarios': scenarios}, attrs=attrs) + + # reference to base year + temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) + # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + try: + hist_xr = xr.open_dataset(modeldir + '/' + hist_file[0],decode_times=False) + ssp1_xr = xr.open_dataset(modeldir + '/' + ssp1_file[0],decode_times=False) + ssp5_xr = xr.open_dataset(modeldir + '/' + ssp5_file[0],decode_times=False) + except: + print('Model misses one of the scenarios used for emulation. Moving on.') + continue + + # some models' time is references to 2015 instead of 1850 + if ssp1_xr['time'][0].values != 60280.5: + ssp1_xr.coords['time'] = ((ssp1_xr['time'] + 60280.5).astype('float')) + ssp5_xr.coords['time'] = ((ssp5_xr['time'] + 60280.5).astype('float')) + + # cocatenate + zos1 = xr.concat([hist_xr, ssp1_xr], dim='time') + zos5 = xr.concat([hist_xr, ssp5_xr], dim='time') + + #% create new coordinates for year + zos1.coords['years'] = ((zos1['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + zos5.coords['years'] = ((zos5['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + + # Calculate annual means + zos1_am = zos1.groupby('years').mean(dim='time') + zos1_am = zos1_am.where(zos1_am['zos'] < 99999) # replace land values with nans + zos5_am = zos5.groupby('years').mean(dim='time') + zos5_am = zos5_am.where(zos5_am['zos'] < 99999) # replace land values with nans + + # reference zos + # baseyear_idx = np.flatnonzero(datayears == baseyear) + # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) + zos1_am = zos1_am.sel(years=targyears) - zos1_am.sel(years=baseyear) + zos5_am = zos5_am.sel(years=targyears) - zos5_am.sel(years=baseyear) + + #% fitting stage + nlon = len(zos1_am['lon']) + nlat = len(zos1_am['lat']) + + + Ts = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='surface').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='surface').drop_vars('scenarios')),dim='years') + + Ti = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='inter').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='inter').drop_vars('scenarios')),dim='years') + + Td = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='deep').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='deep').drop_vars('scenarios')),dim='years') + + + zos1 = np.array(zos1_am['zos'].sel(years=targyears)) + zos5 = np.array(zos5_am['zos'].sel(years=targyears)) + + slope = np.zeros((3, nlat, nlon)) + slope[:] = np.nan + smoothwin = 19 + + reg = linear_model.LinearRegression() + # for i in tqdm(range(nlat)): + # sleep(3) + for i in range(nlat): + for j in range(nlon): + + y = np.concatenate((Smooth(zos1[:,i,j],w=smoothwin), Smooth(zos5[:,i,j],w=smoothwin))) # reducing varability + x1 = np.array(Ts['temperature'][:]).flatten() + x2 = np.array(Ti['temperature'][:]).flatten() + x3 = np.array(Td['temperature'][:]).flatten() + x = np.transpose(np.stack((x1,x2,x3))) + + try: + + reg.fit(x, y) + slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface + slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter + slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep + # intercept[m,i,j] = reg.intercept_ + + except ValueError: + slope[0,i,j] = np.nan + slope[1,i,j] = np.nan + slope[2,i,j] = np.nan + # intercept[m,i,j] = np.nan + + slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? + + lons, lats = np.meshgrid(hist_xr['lon'], hist_xr['lat']) + lats = lats.flatten() + lons = lons.flatten() + + site_ids_map = np.array(NearestPoints(site_lats, site_lons, lats, lons, tol=None)) + + # # remove outliers + # slopes_sum = np.empty((len(slopes), len(lats))) + # # slopes_sum = np.array([np.sum(s,0).flatten() for s in slopes]) + # slopes_sum = np.array([s[0,:,:].flatten() for s in slopes]) + # slopes_sum_std = np.std(slopes_sum,0) + # [ev, ef] = np.unique(np.where(slopes_sum > slopes_sum_std * 3 )[0], return_counts = True) + # mi = [i for i,v in enumerate(ef) if v > 1000] + # slopes = [sl for s, sl in enumerate(slopes) if not s in mi ] + + #nsamps + samples = np.array(gsat['samples']) + + # RESAMPLE SLOPE PARAMETERS + nsims=len(slopes) + + rng = np.random.default_rng(seed) + if nsamps > nsims: + run_idx = np.arange(nsims) + sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) + else: + run_idx = rng.choice(nsims, nsamps, nsamps>nsims) + sample_idx = np.arange(nsamps) + + + ## resample slopes + slopes_resampled = [slopes[i] for i in sample_idx] + + dsl = np.empty((nsamps, len(projyears), len(site_ids))) + dsl[:] = np.nan + + for sample in samples: + slope = slopes_resampled[sample] + print(sample) + + # prepare slopes + slope_s = slope[0,:,:].flatten() + slope_s = np.array([slope_s[x] for x in site_ids_map]) + + slope_i = slope[1,:,:].flatten() + slope_i = np.array([slope_i[x] for x in site_ids_map]) + + slope_d = slope[2,:,:].flatten() + slope_d = np.array([slope_d[x] for x in site_ids_map]) + + # obtain fair temepratures + Tfs = np.array(gsat.sel(samples=sample)) # surface temperature from fair + Tfi = np.array(otemp.sel(layers=1).sel(samples=sample)) # intermediate temperature from fair + Tfd = np.array(otemp.sel(layers=2).sel(samples=sample)) # deep temperature from fair + + + #% project + dsl[sample,:,:] = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) + + # dsl = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) + + ncvar_attributes = {"description": "Dynamic Sea Level"} + + nc_missing_value = np.nan + # Generate the output xarray + dsl_xr = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), dsl, {"units":"mm", "missing_value":nc_missing_value}), + "lat": (("locations"), site_lats), + "lon": (("locations"), site_lons)}, + coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) + + + # make sure projections have 0 mean. This weighted average only works for regular grids. + # Must be edited to accomodate other grids if zos is not in 1x1. + + weights = np.cos(np.deg2rad(dsl_xr.lat)) + weights.name = "weights" + + dsl_xr_weighted = dsl_xr.weighted(weights) + dsl_xr_weighted + weighted_mean = dsl_xr_weighted.mean(("locations")) + + + # Appply 0mean correction to dsl, and add GTE to get sterodynamics + + gte = xr.open_dataset(gte_file) + + + gte = gte['sea_level_change'].values[:,:,0] + sdsl = dsl*1000 - np.repeat(np.array(weighted_mean['sea_level_change'])[:, :, np.newaxis], len(site_ids), axis=2) + np.repeat(gte[0:nsamps, :, np.newaxis], len(site_ids), axis=2) + + ncvar_attributes = {"description": "Local SLR contributions from thermal expansion and dynamic sea-level using EBM3", + "history": "Created " + time.ctime(time.time()), + "source": "SLR Framework: PROTECT 2300", + "scenario": scenario, # change to scenario later + "baseyear": baseyear} + + nc_missing_value = np.nan + # Generate the output xarray + local_out = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), sdsl, {"units":"mm", "missing_value":nc_missing_value}), + "lat": (("locations"), site_lats), + "lon": (("locations"), site_lons)}, + coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) + # Write these samples to a temporary netcdf file + local_out.to_netcdf("{0}_localsl.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) + local_outq = local_out.quantile([0.01,0.05,0.17,0.50,0.83,0.95,0.99], dim='samples') + local_outq.to_netcdf("{0}_quantiles.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) + +if __name__ == '__main__': + + # Initialize the command-line argument parser + parser = argparse.ArgumentParser(description="Compute dynamic sea level for EBM3 workflow",\ + epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") + + # Define the command line arguments to be expected + parser.add_argument('--scenario', help="SSP scenario (i.e ssp585)", default='ssp585') + parser.add_argument('--nsamps', help="Number of samples to generate [default=20000]", default=20000, type=int) + parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) + parser.add_argument('--pyear_start', help="Year for which projections start [default=2000]", default=2000, type=int) + parser.add_argument('--pyear_end', help="Year for which projections end [default=2300]", default=2300, type=int) + parser.add_argument('--pyear_step', help="Step size in years between pyear_start and pyear_end at which projections are produced [default=5]", default=5, type=int) + parser.add_argument('--locationfile', help="File that contains name, id, lat, and lon of points for localization", default="location.lst") + parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) + parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") + parser.add_argument('--climate_data_file',type=str) + parser.add_argument('--rfmip', help='rfmip file',default='rfmip-radiative-forcing-annual-means-v4-0-0.csv') + parser.add_argument('--params', help='CMIP6 Params cvs', default='4xCO2_cummins_ebm3_cmip6.csv') + parser.add_argument('--zosdir',help='Path to CMIP6 ZOS directory', default='cmip6/zos/') + + # Parse the arguments + args = parser.parse_args() + + emb3_thermalexpansion_postprocess(args.scenario, + args.pipeline_id, + args.nsamps, + args.seed, + args.pyear_start, + args.pyear_end, + args.pyear_step, + args.locationfile, + args.baseyear, + args.climate_data_file, + args.rfmip, + args.params, + args.zosdir) + + # Done + sys.exit() \ No newline at end of file diff --git a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py b/modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py similarity index 96% rename from modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py rename to modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py index 6b882dbd..564ce43a 100644 --- a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py +++ b/modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py @@ -1,5 +1,5 @@ -# GMSLR projection program used for IPCC WG1 AR5 -# Translated from IDL to Python 2.7 by Jonathan Gregory 23.10.19 +# Thermal expansion projections for PROTECT +# Victor Malagon Santos, 17.06.2024 import os import numpy as np @@ -138,7 +138,7 @@ def ebm3_thermalexpansion_project(scenario, climate_data_file, coef_file, params if __name__ == '__main__': # Initialize the command-line argument parser - parser = argparse.ArgumentParser(description="Run the glacier projection stage for the AR5 SLR projection workflow",\ + parser = argparse.ArgumentParser(description="Run the global mean thermosteric sea-level rise projections",\ epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") # Define the command line arguments to be expected diff --git a/modules/ebm3/sterodynamics/energy_balance_model.py b/modules/ebm3/sterodynamics/energy_balance_model.py new file mode 100644 index 00000000..b110180b --- /dev/null +++ b/modules/ebm3/sterodynamics/energy_balance_model.py @@ -0,0 +1,519 @@ +"""n-layer energy balance representation of Earth's climate.""" + +import numpy as np +import scipy.linalg +import scipy.sparse.linalg +import scipy.stats +import xarray as xr + +DOUBLING_TIME_1PCT = np.log(2) / np.log(1.01) + +#: Radius of Earth, m +earth_radius = 6371000 + +#: Mass of Earth's atmosphere, kg. +mass_atmosphere = 5.1352e18 + +#: Length of a tropical year, s. +seconds_per_year = 60 * 60 * 24 * 365.24219 + + +class EnergyBalanceModel: + """Energy balance model that converts forcing to temperature. + + The energy balance model is converted to an impulse-response formulation + (hence the IR part of FaIR) to allow efficient evaluation. The benefits of + this are increased as once derived, the "layers" of the energy balance + model do not communicate with each other. The model description can be + found in [Leach2021]_, [Cummins2020]_, [Tsutsui2017]_ and [Geoffroy2013]_. + + Parameters + ---------- + ocean_heat_capacity : ``np.ndarray`` + Ocean heat capacity of each layer (top first), W m-2 yr K-1 + ocean_heat_transfer : ``np.ndarray`` + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : float + efficacy of deepest ocean layer. See e.g. [Geoffroy2013]_. + forcing_4co2 : float + effective radiative forcing from a quadrupling of atmospheric + CO2 concentrations above pre-industrial. + stochastic_run : bool + Activate the stochastic variability component from [Cummins2020]_. + sigma_eta : float + Standard deviation of stochastic forcing component from [Cummins2020]_. + sigma_xi : float + Standard deviation of stochastic disturbance applied to surface + layer. See [Cummins2020]_. + gamma_autocorrelation : float + Stochastic forcing continuous-time autocorrelation parameter. + See [Cummins2020]_. + seed : int or None + Random seed to use for stochastic variability. + timestep : float + Time interval of the model (yr) + + Raises + ------ + ValueError + if the shapes of ``ocean_heat_capacity`` and ``ocean_heat_transfer`` differ. + ValueError + if there are not at least two layers in the energy balance model. + """ + + def __init__( + self, + ocean_heat_capacity, + ocean_heat_transfer, + deep_ocean_efficacy=1, + forcing_4co2=8, + stochastic_run=False, + sigma_eta=0.5, + sigma_xi=0.5, + gamma_autocorrelation=2, + seed=None, + timestep=1, + n_timesteps=1, + ): + """Initialise the EnergyBalanceModel.""" + # adjust ocean heat capacity to be a rate: units W m-2 K-1 + self.ocean_heat_transfer = np.asarray(ocean_heat_transfer) + self.deep_ocean_efficacy = deep_ocean_efficacy + self.forcing_4co2 = forcing_4co2 + self.stochastic_run = stochastic_run + self.sigma_eta = sigma_eta + self.sigma_xi = sigma_xi + self.gamma_autocorrelation = gamma_autocorrelation + self.ocean_heat_capacity = np.asarray(ocean_heat_capacity) / timestep + self.n_temperature_boxes = len(self.ocean_heat_capacity) + if len(self.ocean_heat_transfer) != self.n_temperature_boxes: + raise ValueError( + "ocean_heat_capacity and ocean_heat_transfer must be arrays of the " + "same shape." + ) + if self.n_temperature_boxes < 2: + raise ValueError( + "At least two ocean layers must be specified in the energy balance " + "model." + ) + self.temperature = np.zeros((1, self.n_temperature_boxes + 1)) + self.n_timesteps = n_timesteps + self.n_matrix = self.n_temperature_boxes + 1 + self.seed = seed + self.timestep = timestep + + def _eb_matrix(self): + """Define the matrix of differential equations. + + Returns + ------- + eb_matrix_eigenvalues : `np.ndarray` + 1D array of eigenvalues of the energy balance matrix. + eb_matrix_eigenvectors : `np.ndarray` + 2D array of eigenvectors (an array of 1D eigenvectors) of the + energy balance matrix. + """ + # two box model + # [x x] + # [x x] + + # three box model + # [x x 0] + # [x x ex] + # [0 x x] + + # four box model + # [x x 0 0] + # [x x x 0] + # [0 x x ex] + # [0 0 x x] + + # put the efficacy of deep ocean in the right place + # making a vector avoids if statements + n_box = self.n_temperature_boxes + eb_matrix = np.zeros((n_box, n_box)) + epsilon_array = np.ones(n_box) + epsilon_array[n_box - 2] = self.deep_ocean_efficacy + + # First row + eb_matrix[0, :2] = [ + -( + self.ocean_heat_transfer[0] + + epsilon_array[0] * self.ocean_heat_transfer[1] + ) + / self.ocean_heat_capacity[0], + epsilon_array[0] + * self.ocean_heat_transfer[1] + / self.ocean_heat_capacity[0], + ] + + # Last row + eb_matrix[-1, -2:] = [ + self.ocean_heat_transfer[-1] / self.ocean_heat_capacity[-1], + -self.ocean_heat_transfer[-1] / self.ocean_heat_capacity[-1], + ] + + # Intermediate rows where n>2 + for row in range(1, n_box - 1): + eb_matrix[row, row - 1 : row + 2] = [ + self.ocean_heat_transfer[row] / self.ocean_heat_capacity[row], + -( + self.ocean_heat_transfer[row] + + epsilon_array[row] * self.ocean_heat_transfer[row + 1] + ) + / self.ocean_heat_capacity[row], + epsilon_array[row] + * self.ocean_heat_transfer[row + 1] + / self.ocean_heat_capacity[row], + ] + + # Prepend eb_matrix with stochastic terms if this is a stochastic run: + # Cummins et al. (2020) eqs. 13 and 14 + eb_matrix = np.insert(eb_matrix, 0, np.zeros(n_box), axis=0) + prepend_col = np.zeros(n_box + 1) + prepend_col[0] = -self.gamma_autocorrelation + prepend_col[1] = 1 / self.ocean_heat_capacity[0] + eb_matrix = np.insert(eb_matrix, 0, prepend_col, axis=1) + return eb_matrix + + @property + def eb_matrix_d(self): + """Return the discretised matrix exponential.""" + _eb_matrix_d = scipy.linalg.expm(self._eb_matrix()) + return _eb_matrix_d + + def _forcing_vector(self): + forcing_vector = np.zeros(self.n_temperature_boxes + 1) + forcing_vector[0] = self.gamma_autocorrelation + return forcing_vector + + @property + def forcing_vector_d(self): + """Return the discretised forcing vector.""" + return scipy.linalg.solve( + self._eb_matrix(), + (self.eb_matrix_d - np.identity(self.n_temperature_boxes + 1)) + @ self._forcing_vector(), + ) + + @property + def stochastic_d(self): + """Return the stochastic matrix.""" + # define stochastic matrix + _stochastic_d = np.zeros((self.n_timesteps, self.n_temperature_boxes + 1)) + + # stochastic stuff + if self.stochastic_run: + eb_matrix = self._eb_matrix() + q_mat = np.zeros((self.n_matrix, self.n_matrix)) + q_mat[0, 0] = self.sigma_eta**2 + q_mat[1, 1] = (self.sigma_xi / self.ocean_heat_capacity[0]) ** 2 + # use Van Loan (1978) to compute the matrix exponential + h_mat = np.zeros((self.n_matrix * 2, self.n_matrix * 2)) + h_mat[: self.n_matrix, : self.n_matrix] = -eb_matrix + h_mat[: self.n_matrix, self.n_matrix :] = q_mat + h_mat[self.n_matrix :, self.n_matrix :] = eb_matrix.T + g_mat = scipy.sparse.linalg.expm(h_mat) + q_mat_d = ( + g_mat[self.n_matrix :, self.n_matrix :].T + @ g_mat[: self.n_matrix, self.n_matrix :] + ) + q_mat_d = q_mat_d.astype(np.float64) + _stochastic_d = scipy.stats.multivariate_normal.rvs( + size=self.n_timesteps, + mean=np.zeros(self.n_matrix), + cov=q_mat_d, + random_state=self.seed, + ) + + return _stochastic_d + + def impulse_response(self): + """Convert the energy balance to impulse response representation.""" + eb_matrix = self._eb_matrix() + + # calculate the eigenvectors and eigenvalues on the energy balance + # (determininstic) part of the matrix, these are the timescales of responses + eb_matrix_eigenvalues, eb_matrix_eigenvectors = scipy.linalg.eig( + eb_matrix[1:, 1:] + ) + self.timescales = -self.timestep / (np.real(eb_matrix_eigenvalues)) + self.response_coefficients = ( + self.timescales + * ( + eb_matrix_eigenvectors[0, :] + * scipy.linalg.inv(eb_matrix_eigenvectors)[:, 0] + ) + / (self.ocean_heat_capacity[0] * self.timestep) + ) + + def emergent_parameters(self, forcing_2co2_4co2_ratio=0.5): + """Calculate emergent parameters from the energy balance parameters. + + Parameters + ---------- + forcing_2co2_4co2_ratio : float + ratio of (effective) radiative forcing converting a quadrupling of + CO2 to a doubling of CO2. + """ + # requires impulse response step + if not hasattr(self, "timescales"): + self.impulse_response() + self.ecs = ( + self.forcing_4co2 + * forcing_2co2_4co2_ratio + * np.sum(self.response_coefficients) + ) + self.tcr = ( + self.forcing_4co2 + * forcing_2co2_4co2_ratio + * np.sum( + self.response_coefficients + * ( + 1 + - self.timescales + / DOUBLING_TIME_1PCT + * (1 - np.exp(-DOUBLING_TIME_1PCT / self.timescales)) + ) + ) + ) + + def add_forcing(self, forcing, timestep): + """Add a forcing time series to EnergyBalanceModel. + + Parameters + ---------- + forcing : np.ndarray + time series of [effective] radiative forcing + timestep : float + Model timestep, in years + """ + self.forcing = forcing + self.timestep = timestep + self.n_timesteps = len(forcing) + + def run(self): + """Run the EnergyBalanceModel.""" + # internal variables + n_box = self.n_matrix - 1 + forcing_vector = self._forcing_vector() + + # Calculate the matrix exponential + eb_matrix = self._eb_matrix() + eb_matrix_d = scipy.linalg.expm(eb_matrix) + + # Solve for temperature + forcing_vector_d = scipy.linalg.solve( + eb_matrix, (eb_matrix_d - np.identity(self.n_matrix)) @ forcing_vector + ) + + solution = np.zeros((self.n_timesteps, self.n_matrix)) + solution[0, :] = self.temperature[0, :] + for i in range(1, self.n_timesteps): + solution[i, :] = ( + eb_matrix_d @ solution[i - 1, :] + + forcing_vector_d * self.forcing[i - 1] + + self.stochastic_d[i - 1, :] + ) + + self.temperature = solution[:, 1:] + self.stochastic_forcing = solution[:, 0] + self.toa_imbalance = ( + self.forcing + - self.ocean_heat_transfer[0] * self.temperature[:, 0] + + (1 - self.deep_ocean_efficacy) + * self.ocean_heat_transfer[n_box - 1] + * (self.temperature[:, n_box - 2] - self.temperature[:, n_box - 1]) + ) + self.ocean_heat_content_change = np.cumsum( + self.toa_imbalance + * self.timestep + * earth_radius**2 + * 4 + * np.pi + * seconds_per_year + ) + + +def multi_ebm( + configs, + ocean_heat_capacity, + ocean_heat_transfer, + deep_ocean_efficacy, + stochastic_run, + sigma_eta, + sigma_xi, + gamma_autocorrelation, + seed, + use_seed, + forcing_4co2, + timestep, + timebounds, +): + """Create several instances of the EnergyBalanceModel. + + This allows efficient parallel implementation in FaIR. + We have to use a for loop in this function as is does not look like the linear + algebra functions in scipy are naturally parallel. + + Parameters + ---------- + configs : list + A named list of climate configurations. + ocean_heat_capacity : `np.ndarray` + Ocean heat capacity of each layer (top first), W m-2 yr K-1 + ocean_heat_transfer : `np.ndarray` + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : float + efficacy of deepest ocean layer. See e.g. [Geoffroy2013]_. + stochastic_run : bool + Activate the stochastic variability component from [Cummins2020]_. + sigma_eta : float + Standard deviation of stochastic forcing component from [Cummins2020]_. + sigma_xi : float + Standard deviation of stochastic disturbance applied to surface + layer. See [Cummins2020]_. + gamma_autocorrelation : float + Stochastic forcing continuous-time autocorrelation parameter. + See [Cummins2020]_. + seed : int or None + Random seed to use for stochastic variability. + use_seed : bool + Whether or not to use the random seed. + forcing_4co2 : float + effective radiative forcing from a quadrupling of atmospheric + CO2 concentrations above pre-industrial. + timestep : float + Time interval of the model (yr) + timebounds : np.ndarray + Vector representing the time snapshots to calculate temperature on. + """ + n_configs = ocean_heat_capacity.shape[0] + n_layers = ocean_heat_capacity.shape[1] + n_timebounds = len(timebounds) + ebms = xr.Dataset( + { + "eb_matrix_d": ( + ["config", "eb_dim0", "eb_dim1"], + np.ones((n_configs, n_layers + 1, n_layers + 1)) * np.nan, + ), + "forcing_vector_d": ( + ["config", "eb_dim0"], + np.ones((n_configs, n_layers + 1)) * np.nan, + ), + "stochastic_d": ( + ["timebounds", "config", "eb_dim0"], + np.ones((n_timebounds, n_configs, n_layers + 1)) * np.nan, + ), + "ecs": (["config"], np.ones(n_configs) * np.nan), + "tcr": (["config"], np.ones(n_configs) * np.nan), + }, + coords={ + "timebounds": timebounds, + "config": configs, + "eb_dim0": np.arange(-1, n_layers), + "eb_dim1": np.arange(-1, n_layers), + }, + ) + + for i_conf, config in enumerate(configs): + ebm = EnergyBalanceModel( + ocean_heat_capacity=ocean_heat_capacity[i_conf, :], + ocean_heat_transfer=ocean_heat_transfer[i_conf, :], + deep_ocean_efficacy=deep_ocean_efficacy[i_conf], + stochastic_run=stochastic_run[i_conf], + sigma_eta=sigma_eta[i_conf], + sigma_xi=sigma_xi[i_conf], + gamma_autocorrelation=gamma_autocorrelation[i_conf], + seed=seed.data[i_conf] if use_seed[i_conf] else None, + forcing_4co2=forcing_4co2[i_conf], + timestep=timestep, + n_timesteps=n_timebounds, + ) + ebms["eb_matrix_d"].loc[dict(config=config)] = ebm.eb_matrix_d + ebms["forcing_vector_d"].loc[dict(config=config)] = ebm.forcing_vector_d + ebms["stochastic_d"].loc[dict(config=config)] = ebm.stochastic_d + ebm.emergent_parameters() + ebms["ecs"].loc[dict(config=config)] = ebm.ecs + ebms["tcr"].loc[dict(config=config)] = ebm.tcr + + return ebms + + +def step_temperature(state_old, eb_matrix_d, forcing_vector_d, stochastic_d, forcing): + """Advance parallel energy balance models forward one timestep. + + Parameters + ---------- + state_old : np.ndarray + stacked arrays of forcing and temperature of layers in previous timestep + eb_matrix_d : np.ndarray + stacked discretised energy balance matrices + forcing_vector_d : np.ndarray + stacked discretised forcing vectors + _stochastic_d : np.ndarray + stacked matrices of stochastic internal variability + forcing: np.ndarray + stacked vectors of [effective] radiative forcing + + Returns + ------- + state_new : np.ndarray + stacked arrays of forcing and temperature of layers in this timestep + """ + state_new = ( + (eb_matrix_d[0, ...] @ state_old[0, ..., None])[..., 0] + + forcing_vector_d[0, ...] * forcing[0, ..., 0, None] + + stochastic_d[0, ...] + ) + + return state_new + + +def calculate_toa_imbalance_postrun( + state, + forcing, + ocean_heat_transfer, + deep_ocean_efficacy, +): + """Calculate top of atmosphere energy imbalance. + + The calculation is performed after the scenario has been run to avoid + looping, since no dynamic state changes affect the calculation. + + Parameters + ---------- + state : np.ndarray + stacked arrays of forcing and temperature of layers across the run + forcing : np.ndarray + stacked arrays of [effective] radiative forcing across the run + ocean_heat_transfer : np.ndarray + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : np.ndarray + efficacy of deepest ocean layer. + + Returns + ------- + toa_imbalance : np.ndarray + Top of atmsophere energy imbalance. + """ + toa_imbalance = ( + forcing + - ocean_heat_transfer[..., 0] * state[..., 1] + + (1 - deep_ocean_efficacy) + * ocean_heat_transfer[..., -1] + * (state[..., -2] - state[..., -1]) + ) + return toa_imbalance diff --git a/modules/ebm3/thermalexpansion/pipeline.global.yml b/modules/ebm3/sterodynamics/pipeline.global.yml similarity index 91% rename from modules/ebm3/thermalexpansion/pipeline.global.yml rename to modules/ebm3/sterodynamics/pipeline.global.yml index 65518ae1..4872e770 100644 --- a/modules/ebm3/thermalexpansion/pipeline.global.yml +++ b/modules/ebm3/sterodynamics/pipeline.global.yml @@ -16,6 +16,6 @@ project: - "pyear_step" - "seed" upload_and_extract_input_data: - - "ebm3_thermal_expansion_data.tgz" + - "ebm3_project_data.tgz" global_total_files: - "%PIPELINE_ID%_globalsl.nc" diff --git a/modules/ebm3/thermalexpansion/pipeline.yml b/modules/ebm3/sterodynamics/pipeline.yml similarity index 77% rename from modules/ebm3/thermalexpansion/pipeline.yml rename to modules/ebm3/sterodynamics/pipeline.yml index bb97291b..06586c6d 100644 --- a/modules/ebm3/thermalexpansion/pipeline.yml +++ b/modules/ebm3/sterodynamics/pipeline.yml @@ -16,7 +16,7 @@ project: - "pyear_step" - "seed" upload_and_extract_input_data: - - "ebm3_thermal_expansion_data.tgz" + - "ebm3_project_data.tgz" global_total_files: - "%PIPELINE_ID%_globalsl.nc" @@ -24,7 +24,7 @@ postprocess: task1: executable: "python3" python_dependencies: "scikit-learn tqdm" - script: "ebm3_thermalexpansion_postprocess.py" + script: "ebm3_oceandynamics_project.py" options: - "scenario" - "pipeline_id" @@ -40,15 +40,18 @@ postprocess: - "rfmip" - "params" - "zosdir" + upload_input_data: + - "energy_balance_model.py" copy_input_data: project: task1: - "%PIPELINE_ID%_globalsl.nc" upload_and_extract_input_data: - - "ebm3_thermal_expansion_postprocessing_data.tgz" - - "tlm_sterodynamics_cmip6_data.tgz" + - "ebm3_project_data.tgz" + - "tlm_sterodynamics_cmip6_ebm3_data.tgz" copy_shared_data: - - '$SHARED/location.lst' + - "$SHARED/location.lst" + download_output_data: + - "%PIPELINE_ID%_quantiles.nc" local_total_files: - "%PIPELINE_ID%_localsl.nc" - diff --git a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py b/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py deleted file mode 100644 index bbb87b05..00000000 --- a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py +++ /dev/null @@ -1,490 +0,0 @@ -import numpy as np -import os -import pandas as pd -import sys -import argparse -from fair import energy_balance_model as ebm3 -import csv -import xarray as xr -from sklearn import linear_model -import re -import time -import netCDF4 -import h5py -import scipy -import pickle -from tqdm import tqdm - -""" -Created on Fri Jun 7 11:24:19 2024 - -For any model you have up tp 2300, you want a set of scaling coefficients -for ssp126 and ssp585, as those thwo a typically the ones available up to 2300. -Those scenarios are used to get a pattern of regression coefficient used for -SSP interpolation and probabilistic projections using FAIR. - -@author: vmalagonsantos -""" - -# functions -def fetch_erfs_from_rcmip(path,scenarios): - ''' - function to read scenario-dependent ERF from RCMIP csv - - Input parameters: - path = path to RCMIP csv - scenarios = desired ssp's - - Output: Adds results of extreme sea-level analysis for all stations to the - station data dictionary. - - erfs = effective radiative forcing per scenario and year - erfyears = corresponding years (1750-2500) - ''' - - erfyears=np.arange(1750,2501) - ssp_idx = {'ssp119':212,'ssp126':231,'ssp245':308,'ssp370':59,'ssp585':404} #location in table - - with open(path) as csv_file: - csv_reader = csv.reader(csv_file) - rows = list(csv_reader) - - erfs = np.empty((len(scenarios),len(erfyears))) - - for s,scen in enumerate(scenarios): - try: - erfs[s,:] = rows[ssp_idx[scen]][7:] - except: - continue - - return erfs, erfyears - -def angd(lat0, lon0, qlat, qlon): - - # Convert the input from degrees to radians - (lat0, lon0) = np.radians((lat0, lon0)) - (qlat, qlon) = np.radians((qlat, qlon)) - - # Calculate the angle between the vectors - temp = np.arctan2(np.sqrt((np.cos(qlat)*np.sin(qlon-lon0))**2 + \ - (np.cos(lat0)*np.sin(qlat) - np.sin(lat0)*np.cos(qlat) * np.cos(qlon-lon0))**2),\ - (np.sin(lat0)*np.sin(qlat) + np.cos(lat0)*np.cos(qlat)*np.cos(qlon-lon0))) - - # Convert the results from radians to degrees and return - return(np.degrees(temp)) - -def NearestPoint(qlat, qlon, lats, lons, tol = None): - - # Get the distance between the query point and all the possible points - dist = angd(lats, lons, qlat, qlon) - - # Which is the closest point - nearest_idx = np.argmin(dist) - - # Is the point within the tolerance? - if isinstance(tol, (int, float)): - if dist[nearest_idx] > tol: - return(None) - - return(nearest_idx) - - -def NearestPoints(qlats, qlons, lats, lons, tol): - - if len(qlats) != len(qlons): - raise Exception("Query lats ({}) and lons ({}) differ in length".format(len(qlats), len(qlons))) - - idx = map(lambda qlat,qlon: NearestPoint(qlat, qlon, lats, lons, tol), qlats, qlons) - - return(list(idx)) - -def ReadLocationFile(location_file): - - # Initialize variables to hold data and site information - names = [] - ids = [] - lats = [] - lons = [] - - # Compile the regex for finding commented lines - comment_regex = re.compile(r'^#') - - # Open the rate file - with open(location_file, 'r') as f: - - # Loop over the lines of the file - for line in f: - - # Skip commented lines - if re.search(comment_regex, line): - continue - - # Split the line into components - (this_name, this_id, this_lat, this_lon) = line.split("\t") - - # Store the information - names.append(this_name) - ids.append(int(this_id)) - lats.append(float(this_lat)) - lons.append(float(this_lon)) - - # Cast everything as numpy arrays - names = np.array(names) - ids = np.array(ids) - lats = np.array(lats) - lons = np.array(lons) - - # Return variables - return(names, ids, lats, lons) - - -# Smooth ZOS and ZOSTOGA over 19 year smoothing window -def nanSmooth(x, w=19): - idx = np.flatnonzero(~np.isnan(x)) - temp = x - if len(idx) > 0: - temp[idx] = Smooth(x[idx], w) - return(temp) - -def Smooth(x, w=19): - out0 = np.convolve(x, np.ones(w,dtype='double'), 'valid')/w - r = np.arange(1,w-1,2, dtype="double") - start = np.cumsum(x[:w-1])[::2]/r - stop = (np.cumsum(x[:-w:-1])[::2]/r)[::-1] - y = np.concatenate((start, out0, stop)) - return(y) - -def emb3_thermalexpansion_postprocess(scenario, pipeline_id, nsamps, seed, pyear_start, pyear_end, pyear_step, locationfile, baseyear, climate_data_file, rfmip, params, zosdir): - targyears = np.arange(baseyear,2301) # for regression - projyears = np.arange(pyear_start,pyear_end+1, pyear_step) - (_, site_ids, site_lats, site_lons) = ReadLocationFile(locationfile) - - # get temperature from FaIR simulations: INPUT FROM CLIMATE STEP. We need both gsat and oceantemp - cf = xr.open_dataset(climate_data_file, group=scenario, engine='netcdf4', ) - gsat = cf['surface_temperature'].sel(years=projyears) - cf['surface_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - otemp = cf['deep_ocean_temperature'].sel(years=projyears) - cf['deep_ocean_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - - # INPUT temperature file from ebm3 global - gte_file = f'{pipeline_id}_globalsl.nc' - - # INPUT get models and parameters - paramdir = params - - #forcing for ebm - scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections - rfmipfile = rfmip #path # RFMIP FILE IS INPUT, NEW DATA - erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries - - ebm_param = pd.read_csv(paramdir) - - # find available zos simulations - ignoring hidden files - model_has_zos = [f for f in os.listdir(zosdir) if not f.startswith('.')] - - ## FITTING - slopes = [] - intercept = [] - - for m, model in enumerate(model_has_zos): - # model = model_has_zos[0] - - # model = 'MRI-ESM2-0' - - # find matching ebm parameters to model - model_param = ebm_param[ebm_param["model"].str.contains(model)] - - modeldir = zosdir + model - # keep those up to 2300 or 2500 - zos_runs2_2300 = [f for f in os.listdir(modeldir) if '2300' in f or '2500' in f] # regex would be better here - if not zos_runs2_2300: - print('Model ' + model + ' only runs to 2100. Skipping and looking for longer simulations' ) - - continue - - print('Estimating parameters for model: ' + model) - - - # get variants with 2300 simulation - if len(zos_runs2_2300) == 1: - variant_2300 = [f for f in re.split('_', zos_runs2_2300[0]) if f.startswith('r')] - model_var_param = model_param[ebm_param["run"].str.contains(variant_2300[0])] - else: - for v, filename in enumerate(zos_runs2_2300): - variant_2300 = [f for f in re.split('_', zos_runs2_2300[v]) if f.startswith('r')] - # get EBM parameters for both model and variant - model_var_param = model_param[ebm_param["run"].str.contains(variant_2300[0])] - if not model_var_param.empty: - continue - - c1 = np.array(model_var_param.get('C1')) - c2 = np.array(model_var_param.get('C2')) - c3 = np.array(model_var_param.get('C3')) - k1 = np.array(model_var_param.get('kappa1')) - k2 = np.array(model_var_param.get('kappa2')) - k3 = np.array(model_var_param.get('kappa3')) - e = np.array(model_var_param.get('epsilon')) - f4 = np.array(model_var_param.get('F_4xCO2')) - se = np.array(model_var_param.get('sigma_eta')) - sx = np.array(model_var_param.get('sigma_xi')) - gamma = np.array(model_var_param.get('gamma')) - nit = np.array(model_var_param.get('nit')) - - - temp3 = np.empty((len(erfyears), 3, len(scenarios))) - - for s, scen in enumerate(scenarios): - - #Model setup 3lm - ebm_obj = ebm3.EnergyBalanceModel( - ocean_heat_capacity = [c1[0], c2[0], c3[0]], - ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value - deep_ocean_efficacy = e[0], - forcing_4co2 = f4[0], - stochastic_run=False, - sigma_eta = se[0], - sigma_xi = sx[0], - gamma_autocorrelation = gamma[0], - seed = None, - timestep = 1, - n_timesteps = 2500-1750, - ) - - # Run model for different scenarios - ebm_obj.add_forcing( erfs[s,:], erfyears) - ebm_obj.run() - temp3[:,:,s] = ebm_obj.temperature - - attrs={'description': 'surface and deep temperature obtained by EBM-3LM', - 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} - - temp3_xr = xr.Dataset({'temperature':(['years','layers', 'scenarios'], temp3)}, - coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears, 'scenarios': scenarios}, attrs=attrs) - - # reference to base year - temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) - # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - - # get zos, historical + ssp - model_dir = zosdir + model - - hist_file = [f for f in os.listdir(model_dir) if 'hist' in f and variant_2300[0] in f] - ssp1_file = [f for f in os.listdir(model_dir) if 'ssp1' in f and variant_2300[0] in f and 'ssp126' in f] - ssp5_file = [f for f in os.listdir(model_dir) if 'ssp5' in f and variant_2300[0] in f and 'ssp585' in f] - - if not ssp1_file or not ssp5_file: - continue - - hist_xr = xr.open_dataset(model_dir + '/' + hist_file[0],decode_times=False) - ssp1_xr = xr.open_dataset(model_dir + '/' + ssp1_file[0],decode_times=False) - ssp5_xr = xr.open_dataset(model_dir + '/' + ssp5_file[0],decode_times=False) - - if ssp1_xr['time'][0].values != 60280.5: - ssp1_xr.coords['time'] = ((ssp1_xr['time'] + 60280.5).astype('float')) - ssp5_xr.coords['time'] = ((ssp5_xr['time'] + 60280.5).astype('float')) - - # cocatenate - zos1 = xr.concat([hist_xr, ssp1_xr], dim='time') - zos5 = xr.concat([hist_xr, ssp5_xr], dim='time') - - #% create new coordinates for year - zos1.coords['years'] = ((zos1['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR - zos5.coords['years'] = ((zos5['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR - - # Calculate annual means - zos1_am = zos1.groupby('years').mean(dim='time') - zos1_am = zos1_am.where(zos1_am['zos'] < 99999) # replace land values with nans - zos5_am = zos5.groupby('years').mean(dim='time') - zos5_am = zos5_am.where(zos5_am['zos'] < 99999) # replace land values with nans - - # reference zos - # baseyear_idx = np.flatnonzero(datayears == baseyear) - # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) - zos1_am = zos1_am.sel(years=targyears) - zos1_am.sel(years=baseyear) - zos5_am = zos5_am.sel(years=targyears) - zos5_am.sel(years=baseyear) - - #% fitting stage - nlon = len(zos1_am['lon']) - nlat = len(zos1_am['lat']) - - - Ts = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='surface').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='surface').drop_vars('scenarios')),dim='years') - - Ti = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='inter').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='inter').drop_vars('scenarios')),dim='years') - - Td = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='deep').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='deep').drop_vars('scenarios')),dim='years') - - - zos1 = np.array(zos1_am['zos'].sel(years=targyears)) - zos5 = np.array(zos5_am['zos'].sel(years=targyears)) - - slope = np.zeros((3, nlat, nlon)) - slope[:] = np.nan - smoothwin = 19; - - reg = linear_model.LinearRegression() - # for i in tqdm(range(nlat)): - # sleep(3) - for i in tqdm(range(nlat)): - for j in range(nlon): - - y = np.concatenate((Smooth(zos1[:,i,j],w=smoothwin), Smooth(zos5[:,i,j],w=smoothwin))) # reducing varability - x1 = np.array(Ts['temperature'][:]).flatten() - x2 = np.array(Ti['temperature'][:]).flatten() - x3 = np.array(Td['temperature'][:]).flatten() - x = np.transpose(np.stack((x1,x2,x3))) - - try: - - reg.fit(x, y) - slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface - slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter - slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep - # intercept[m,i,j] = reg.intercept_ - - except ValueError: - slope[0,i,j] = np.nan - slope[1,i,j] = np.nan - slope[2,i,j] = np.nan - # intercept[m,i,j] = np.nan - - slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? - - lons, lats = np.meshgrid(zos1_am['lon'], zos1_am['lat']) - lats = lats.flatten() - lons = lons.flatten() - - site_ids_map = np.array(NearestPoints(site_lats, site_lons, lats, lons, tol=None)) - - #nsamps = 100 - samples = np.array(gsat['samples']) - - # RESAMPLE SLOPE PARAMETERS - nsims=len(slopes) - - rng = np.random.default_rng(seed) - if nsamps > nsims: - run_idx = np.arange(nsims) - sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) - else: - run_idx = rng.choice(nsims, nsamps, nsamps>nsims) - sample_idx = np.arange(nsamps) - - - ## need to figure out of to do this resamples properly, so from now on this is dummy data - slopes_resampled = [slopes[i] for i in sample_idx] - # aa = np.dstack(slopes_resampled) - - dsl = np.empty((nsamps, len(projyears), len(site_ids))) - dsl[:] = np.nan - - for sample in tqdm(samples): - slope = slopes_resampled[sample] - print(sample) - - # prepare slopes - slope_s = slope[0,:,:].flatten() - slope_s = np.array([slope_s[x] for x in site_ids_map]) - - slope_i = slope[1,:,:].flatten() - slope_i = np.array([slope_i[x] for x in site_ids_map]) - - slope_d = slope[2,:,:].flatten() - slope_d = np.array([slope_d[x] for x in site_ids_map]) - - # obtain fair temepratures - Tfs = np.array(gsat.sel(samples=sample)) # surface temperature from fair - Tfi = np.array(otemp.sel(layers=1).sel(samples=sample)) # intermediate temperature from fair - Tfd = np.array(otemp.sel(layers=2).sel(samples=sample)) # deep temperature from fair - - - #% project - dsl[sample,:,:] = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) - - # dsl = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) - - ncvar_attributes = {"description": "Dynamic Sea Level"} - - nc_missing_value = np.nan - # Generate the output xarray - dsl_xr = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), dsl, {"units":"mm", "missing_value":nc_missing_value}), - "lat": (("locations"), site_lats), - "lon": (("locations"), site_lons)}, - coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) - - - # make sure projections have 0 mean. This weighted average only works for regular grids. - # Must be edited to accomodate other grids if zos is not in 1x1. - - weights = np.cos(np.deg2rad(dsl_xr.lat)) - weights.name = "weights" - - dsl_xr_weighted = dsl_xr.weighted(weights) - dsl_xr_weighted - weighted_mean = dsl_xr_weighted.mean(("locations")) - - - # Appply 0mean correction to dsl, and add GTE to get sterodynamics - - gte = xr.open_dataset(gte_file) - - - gte = gte['sea_level_change'].values[:,:,0] - sdsl = dsl*1000 - np.repeat(np.array(weighted_mean['sea_level_change'])[:, :, np.newaxis], len(site_ids), axis=2) + np.repeat(gte[0:nsamps, :, np.newaxis], len(site_ids), axis=2) - - ncvar_attributes = {"description": "Local SLR contributions from thermal expansion and dynamic sea-level using EBM3", - "history": "Created " + time.ctime(time.time()), - "source": "SLR Framework: PROTECT 2300", - "scenario": scenario, # change to scenario later - "baseyear": baseyear} - - nc_missing_value = np.nan - # Generate the output xarray - local_out = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), sdsl, {"units":"mm", "missing_value":nc_missing_value}), - "lat": (("locations"), site_lats), - "lon": (("locations"), site_lons)}, - coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) - # Write these samples to a temporary netcdf file - local_out.to_netcdf("{0}_localsl.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) - -if __name__ == '__main__': - - # Initialize the command-line argument parser - parser = argparse.ArgumentParser(description="Run the pre-processing stage for the TLM ocean dynamics workflow",\ - epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") - - # Define the command line arguments to be expected - parser.add_argument('--scenario', help="SSP scenario (i.e ssp585) or temperature target (i.e. tlim2.0win0.25)", default='ssp585') - parser.add_argument('--nsamps', help="Number of samples to generate [default=20000]", default=20000, type=int) - parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) - parser.add_argument('--pyear_start', help="Year for which projections start [default=2000]", default=2000, type=int) - parser.add_argument('--pyear_end', help="Year for which projections end [default=2300]", default=2300, type=int) - parser.add_argument('--pyear_step', help="Step size in years between pyear_start and pyear_end at which projections are produced [default=5]", default=5, type=int) - parser.add_argument('--locationfile', help="File that contains name, id, lat, and lon of points for localization", default="location.lst") - parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) - parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") - parser.add_argument('--climate_data_file',type=str) - parser.add_argument('--rfmip', help='rfmip file',default='rfmip-radiative-forcing-annual-means-v4-0-0.csv') - parser.add_argument('--params', help='CMIP6 Params cvs', default='4xCO2_cummins_ebm3_cmip6.csv') - parser.add_argument('--zosdir',help='Path to CMIP6 ZOS directory', default='cmip6/zos/') - - # Parse the arguments - args = parser.parse_args() - - emb3_thermalexpansion_postprocess(args.scenario, - args.pipeline_id, - args.nsamps, - args.seed, - args.pyear_start, - args.pyear_end, - args.pyear_step, - args.locationfile, - args.baseyear, - args.climate_data_file, - args.rfmip, - args.params, - args.zosdir) - - # Done - sys.exit() \ No newline at end of file diff --git a/modules/emulandice2/AIS/pipeline.2300.yml b/modules/emulandice2/AIS/pipeline.2300.yml new file mode 100644 index 00000000..06ea1919 --- /dev/null +++ b/modules/emulandice2/AIS/pipeline.2300.yml @@ -0,0 +1,88 @@ + +project: + task1: + python_dependencies: "xarray dask[complete] numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - AIS + - --region + - ALL + - --emu_file + - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_AIS_project_2300_240528_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_EAIS1_globalsl.nc" + - "%PIPELINE_ID%_EAIS2_globalsl.nc" + - "%PIPELINE_ID%_EAIS3_globalsl.nc" + - "%PIPELINE_ID%_EAIS4_globalsl.nc" + - "%PIPELINE_ID%_EAIS5_globalsl.nc" + - "%PIPELINE_ID%_EAIS6_globalsl.nc" + - "%PIPELINE_ID%_EAIS7_globalsl.nc" + - "%PIPELINE_ID%_PEN_globalsl.nc" + - "%PIPELINE_ID%_WAIS1_globalsl.nc" + - "%PIPELINE_ID%_WAIS2_globalsl.nc" + - "%PIPELINE_ID%_WAIS3_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data2.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_EAIS1_globalsl.nc" + - "%PIPELINE_ID%_EAIS2_globalsl.nc" + - "%PIPELINE_ID%_EAIS3_globalsl.nc" + - "%PIPELINE_ID%_EAIS4_globalsl.nc" + - "%PIPELINE_ID%_EAIS5_globalsl.nc" + - "%PIPELINE_ID%_EAIS6_globalsl.nc" + - "%PIPELINE_ID%_EAIS7_globalsl.nc" + - "%PIPELINE_ID%_PEN_globalsl.nc" + - "%PIPELINE_ID%_WAIS1_globalsl.nc" + - "%PIPELINE_ID%_WAIS2_globalsl.nc" + - "%PIPELINE_ID%_WAIS3_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/AIS/pipeline.global.2150.yml b/modules/emulandice2/AIS/pipeline.global.2150.yml index 88b52e4d..61ec5188 100644 --- a/modules/emulandice2/AIS/pipeline.global.2150.yml +++ b/modules/emulandice2/AIS/pipeline.global.2150.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_2150_240214_data.tgz" + - "emulandice2_AIS_project_2150_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/AIS/pipeline.global.2300.yml b/modules/emulandice2/AIS/pipeline.global.2300.yml index 2752d2db..6e9364f8 100644 --- a/modules/emulandice2/AIS/pipeline.global.2300.yml +++ b/modules/emulandice2/AIS/pipeline.global.2300.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_2300_240214_data.tgz" + - "emulandice2_AIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/AIS/pipeline.global.yml b/modules/emulandice2/AIS/pipeline.global.yml index af3fa527..61ec5188 100644 --- a/modules/emulandice2/AIS/pipeline.global.yml +++ b/modules/emulandice2/AIS/pipeline.global.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + - emu_file/AIS_ALL_Kori_PISM_CISM_ElmerIce_pow_exp_10_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_data.tgz" + - "emulandice2_AIS_project_2150_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/AIS/pipeline.yml b/modules/emulandice2/AIS/pipeline.yml index 50654cb7..c3f6ec22 100644 --- a/modules/emulandice2/AIS/pipeline.yml +++ b/modules/emulandice2/AIS/pipeline.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + - emu_file/AIS_ALL_Kori_PISM_CISM_ElmerIce_pow_exp_10_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -27,7 +27,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_data.tgz" + - "emulandice2_AIS_project_2150_240214_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.2300.yml b/modules/emulandice2/GrIS/pipeline.2300.yml new file mode 100644 index 00000000..6a332691 --- /dev/null +++ b/modules/emulandice2/GrIS/pipeline.2300.yml @@ -0,0 +1,80 @@ + +project: + task1: + python_dependencies: "xarray dask[complete] numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - GIS + - --region + - ALL + - --emu_file + - emu_file/GIS_ALL_CISM_pow_exp_01_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + - "cyear_start" + - "cyear_end" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_GrIS_project_2300_240528_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_CW_globalsl.nc" + - "%PIPELINE_ID%_SW_globalsl.nc" + - "%PIPELINE_ID%_SE_globalsl.nc" + - "%PIPELINE_ID%_NE_globalsl.nc" + - "%PIPELINE_ID%_NO_globalsl.nc" + - "%PIPELINE_ID%_NW_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data2.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_CW_globalsl.nc" + - "%PIPELINE_ID%_SW_globalsl.nc" + - "%PIPELINE_ID%_SE_globalsl.nc" + - "%PIPELINE_ID%_NE_globalsl.nc" + - "%PIPELINE_ID%_NO_globalsl.nc" + - "%PIPELINE_ID%_NW_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/GrIS/pipeline.global.2100.yml b/modules/emulandice2/GrIS/pipeline.global.2100.yml index 1932fb9b..ef53561b 100644 --- a/modules/emulandice2/GrIS/pipeline.global.2100.yml +++ b/modules/emulandice2/GrIS/pipeline.global.2100.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_2100_240214_data.tgz" + - "emulandice2_GrIS_project_2100_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.global.2300.yml b/modules/emulandice2/GrIS/pipeline.global.2300.yml index 9eba19f9..3d265a37 100644 --- a/modules/emulandice2/GrIS/pipeline.global.2300.yml +++ b/modules/emulandice2/GrIS/pipeline.global.2300.yml @@ -1,7 +1,7 @@ project: task1: - python_dependencies: "xarray dask numpy scipy" + python_dependencies: "xarray dask[complete] numpy scipy" executable: "python3" script_noupload: "emulandice_project.py" arguments: @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_2300_240214_data.tgz" + - "emulandice2_GrIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.yml b/modules/emulandice2/GrIS/pipeline.yml index e5cf200a..aa9f0815 100644 --- a/modules/emulandice2/GrIS/pipeline.yml +++ b/modules/emulandice2/GrIS/pipeline.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/GIS_ALL_CISM_pow_exp_01_EMULATOR.RData + - emu_file/GIS_ALL_CISM_IMAUICE_ElmerIce_GISM_pow_exp_01_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -27,7 +27,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_data.tgz" + - "emulandice2_GrIS_project_2100_240214_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/emulandice_config.sh b/modules/emulandice2/emulandice_config.sh index 291309a5..0344473b 100755 --- a/modules/emulandice2/emulandice_config.sh +++ b/modules/emulandice2/emulandice_config.sh @@ -1,7 +1,8 @@ #!/bin/bash +set -e -ROOTDIR=`dirname $0` -MODE=$1 +ROOTDIR="$(dirname "$0")" +MODE=${1:-} case $MODE in "clean") diff --git a/modules/emulandice2/glaciers/pipeline.2300.yml b/modules/emulandice2/glaciers/pipeline.2300.yml new file mode 100644 index 00000000..01e7c4a3 --- /dev/null +++ b/modules/emulandice2/glaciers/pipeline.2300.yml @@ -0,0 +1,141 @@ + +project: + task1: + python_dependencies: "xarray dask numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - GLA + - --region + - RGI01 + - RGI02 + - RGI03 + - RGI04 + - RGI05 + - RGI06 + - RGI07 + - RGI08 + - RGI09 + - RGI10 + - RGI11 + - RGI12 + - RGI13 + - RGI14 + - RGI15 + - RGI16 + - RGI17 + - RGI18 + - RGI19 + - --emu_file + - emu_file/GLA_RGI01_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI02_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI03_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI04_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI05_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI06_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI07_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI08_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI09_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI10_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI11_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI12_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI13_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI14_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI15_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI16_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI17_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI18_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI19_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + - "cyear_start" + - "cyear_end" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_glaciers_project_2300_240320_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_RGI01_globalsl.nc" + - "%PIPELINE_ID%_RGI02_globalsl.nc" + - "%PIPELINE_ID%_RGI03_globalsl.nc" + - "%PIPELINE_ID%_RGI04_globalsl.nc" + - "%PIPELINE_ID%_RGI05_globalsl.nc" + - "%PIPELINE_ID%_RGI06_globalsl.nc" + - "%PIPELINE_ID%_RGI07_globalsl.nc" + - "%PIPELINE_ID%_RGI08_globalsl.nc" + - "%PIPELINE_ID%_RGI09_globalsl.nc" + - "%PIPELINE_ID%_RGI10_globalsl.nc" + - "%PIPELINE_ID%_RGI11_globalsl.nc" + - "%PIPELINE_ID%_RGI12_globalsl.nc" + - "%PIPELINE_ID%_RGI13_globalsl.nc" + - "%PIPELINE_ID%_RGI14_globalsl.nc" + - "%PIPELINE_ID%_RGI15_globalsl.nc" + - "%PIPELINE_ID%_RGI16_globalsl.nc" + - "%PIPELINE_ID%_RGI17_globalsl.nc" + - "%PIPELINE_ID%_RGI18_globalsl.nc" + - "%PIPELINE_ID%_RGI19_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_RGI01_globalsl.nc" + - "%PIPELINE_ID%_RGI02_globalsl.nc" + - "%PIPELINE_ID%_RGI03_globalsl.nc" + - "%PIPELINE_ID%_RGI04_globalsl.nc" + - "%PIPELINE_ID%_RGI05_globalsl.nc" + - "%PIPELINE_ID%_RGI06_globalsl.nc" + - "%PIPELINE_ID%_RGI07_globalsl.nc" + - "%PIPELINE_ID%_RGI08_globalsl.nc" + - "%PIPELINE_ID%_RGI09_globalsl.nc" + - "%PIPELINE_ID%_RGI10_globalsl.nc" + - "%PIPELINE_ID%_RGI11_globalsl.nc" + - "%PIPELINE_ID%_RGI12_globalsl.nc" + - "%PIPELINE_ID%_RGI13_globalsl.nc" + - "%PIPELINE_ID%_RGI14_globalsl.nc" + - "%PIPELINE_ID%_RGI15_globalsl.nc" + - "%PIPELINE_ID%_RGI16_globalsl.nc" + - "%PIPELINE_ID%_RGI17_globalsl.nc" + - "%PIPELINE_ID%_RGI18_globalsl.nc" + - "%PIPELINE_ID%_RGI19_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/glaciers/pipeline.global.2300.yml b/modules/emulandice2/glaciers/pipeline.global.2300.yml index 47b0fad1..b2b657b7 100644 --- a/modules/emulandice2/glaciers/pipeline.global.2300.yml +++ b/modules/emulandice2/glaciers/pipeline.global.2300.yml @@ -65,7 +65,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_glaciers_project_2300_240223_data.tgz" + - "emulandice2_glaciers_project_2300_240320_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/glaciers/pipeline.yml b/modules/emulandice2/glaciers/pipeline.yml index 4be08ffc..b8f21b30 100644 --- a/modules/emulandice2/glaciers/pipeline.yml +++ b/modules/emulandice2/glaciers/pipeline.yml @@ -63,7 +63,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_glaciers_project_data.tgz" + - "emulandice2_glaciers_project_2100_240223_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/shared/emulandice2 b/modules/emulandice2/shared/emulandice2 index fd319046..0a646f6c 160000 --- a/modules/emulandice2/shared/emulandice2 +++ b/modules/emulandice2/shared/emulandice2 @@ -1 +1 @@ -Subproject commit fd3190469636e14bf7e6c0d1ac18a14a052fd634 +Subproject commit 0a646f6c66a86d6ffbbea9ced676236ba3405501 diff --git a/modules/emulandice2/shared/grd_fingerprintmap.yml b/modules/emulandice2/shared/grd_fingerprintmap.yml index a539e8a0..0b6f47da 100644 --- a/modules/emulandice2/shared/grd_fingerprintmap.yml +++ b/modules/emulandice2/shared/grd_fingerprintmap.yml @@ -6,18 +6,41 @@ # Should go to modules-data when finished AIS: - AIS.ALL: - fingerprint: FPRINT/fprint_wais.nc -# AIS.EAIS: -# fingerprint: FPRINT/fprint_eais.nc -# AIS.WAIS: -# fingerprint: FPRINT/fprint_wais.nc -# AIS.AP: -# fingerprint: FPRINT/fprint_wais.nc + AIS.WAIS1: + fingerprint: FPRINT/fprint_ais11_1.nc + AIS.WAIS2: + fingerprint: FPRINT/fprint_ais11_2.nc + AIS.WAIS3: + fingerprint: FPRINT/fprint_ais11_3.nc + AIS.EAIS1: + fingerprint: FPRINT/fprint_ais11_5.nc + AIS.EAIS2: + fingerprint: FPRINT/fprint_ais11_6.nc + AIS.EAIS3: + fingerprint: FPRINT/fprint_ais11_7.nc + AIS.PEN: + fingerprint: FPRINT/fprint_ais11_4.nc + AIS.EAIS4: + fingerprint: FPRINT/fprint_ais11_8.nc + AIS.EAIS5: + fingerprint: FPRINT/fprint_ais11_9.nc + AIS.EAIS6: + fingerprint: FPRINT/fprint_ais11_10.nc + AIS.EAIS7: + fingerprint: FPRINT/fprint_ais11_11.nc GrIS: - GIS.ALL: - fingerprint: FPRINT/fprint_gis.nc - + GIS.NW: + fingerprint: FPRINT/fprint_gis01.nc + GIS.CW: + fingerprint: FPRINT/fprint_gis02.nc + GIS.SW: + fingerprint: FPRINT/fprint_gis03.nc + GIS.SE: + fingerprint: FPRINT/fprint_gis04.nc + GIS.NE: + fingerprint: FPRINT/fprint_gis05.nc + GIS.NO: + fingerprint: FPRINT/fprint_gis06.nc glaciers: GLA.RGI01: fingerprint: FPRINT/fprint_glac11.nc diff --git a/modules/oelsmann24/verticallandmotion/vlm_preprocess.py b/modules/oelsmann24/verticallandmotion/vlm_preprocess.py index b9598b4f..4014a76e 100644 --- a/modules/oelsmann24/verticallandmotion/vlm_preprocess.py +++ b/modules/oelsmann24/verticallandmotion/vlm_preprocess.py @@ -596,7 +596,7 @@ def interpolate_and_compute_weights(VLM_REC_MERGED,VLM_REC_MERGED_GRW,GIA_VLM,di lat_grid = np.arange((180*4))/4-90 distance_max = 150 - lon_ext = dist2coast['lon'].values + lon_ext = dist2coast['lon'].values.copy() lon_ext[lon_ext>180]=lon_ext[lon_ext>180]-360 dist2coast.assign_coords(lon=dist2coast.lon*0+lon_ext) @@ -631,7 +631,7 @@ def interpolate_and_compute_weights(VLM_REC_MERGED,VLM_REC_MERGED_GRW,GIA_VLM,di w_d = w_d.fillna(0) - lon_prt =GIA_VLM['lon'].values + lon_prt =GIA_VLM['lon'].values.copy() lon_prt[lon_prt>180] = lon_prt[lon_prt>180]-360 GIA_VLM['lon'] = lon_prt