''' calculate values from cmorized CDRMIP runs ''' import copy import numpy as np import matplotlib.pyplot as plt import xarray as xr import glob import pandas as pd ################# flux, mass = False, False weighted_mean = False save_to_csv = True #for ireal in range(1,11): for ireal in [1,2,3]: real = 'r' + str(ireal) ############# #model_list = ['ACCESS-ESM1-5','CanESM5','MPI-ESM1-2-LR','NorESM2-LM','UKESM1-0-LL'] #model_list = ['ACCESS-ESM1-5','CanESM5','CESM2']#['MPI-ESM1-2-LR','UKESM1-0-LL'] #model_list = ['CanESM5'] #model_list = ['ACCESS-ESM1-5'] #model_list = ['UKESM1-0-LL'] #model_list = ['MIROC-ES2L'] #model_list = ['CESM2'] #model_list = ['MPI-ESM1-2-LR'] model_list = ['NorESM2-LM'] #exp_list = ['esm-ssp585'] #exp_list = ['esm-ssp585-ssp126Lu'] #exp_list = ['esm-ssp585','esm-ssp585-ssp126Lu'] #exp_list = ['esm-ssp585','esm-ssp585-ocn-alk'] exp_list = ['esm-ssp585-ocn-alk'] #var_list = ['co2_concentration'] var_list = ['cLand'] #var_list = ['fgco2'] #var_list = ['tree'] #var_list = ['tas'] #var_list = ['nbp'] ##### molC_to_kgC = 0.01201 kgCO2_to_kgC = 0.272912 ### dirname = '/scratch/usr/shkhwwey/CDRMIP/' ### for model_name in model_list: for exp_name in exp_list: if exp_name in ['esm-ssp585-ocn-alk']: if model_name in ['NorESM2-LM']: yrstart = 2015 else: yrstart = 2020 if exp_name in ['esm-ssp585','esm-ssp585-ssp126Lu']: yrstart = 2015 for var_name in var_list: if var_name =='co2_concentration': if not(model_name in ['NorESM2-LM','CESM2','MIROC-ES2L','ACCESS-ESM1-5']): var_name = 'co2mass' elif model_name=='ACCESS-ESM1-5': var_name = 'co23D' else: var_name = 'co2' if var_name =='tree': var_name = 'treeFrac' full_dir = dirname + model_name + '/' + exp_name + '/' + var_name + '/' if real != '': full_dir = full_dir + 'real/' file_list = glob.glob(full_dir + '*' + real + 'i1*.nc') else: file_list = glob.glob(full_dir + '*.nc') file_list.sort() for file_name in file_list: print(file_name.split('/')[-1]) ds = xr.open_mfdataset(file_list) print(model_name,exp_name,ds[var_name].shape) if var_name in ['cLand','treeFrac']: # unit: kg/m2; % mass,flux = True, False file_name = glob.glob(dirname + model_name + '/areacella*.nc')[0] cell_area = xr.open_mfdataset(file_name)['areacella'].to_masked_array() if var_name in ['nbp']: # unit: kg m-2 s-1 flux,mass = True, False file_name = glob.glob(dirname + model_name + '/areacella*.nc')[0] cell_area = xr.open_mfdataset(file_name)['areacella'].to_masked_array() if var_name in ['fgco2']: # unit: kg m-2 s-1 flux,mass = True, False file_name = glob.glob(dirname + model_name + '/areacello*.nc')[0] cell_area = xr.open_mfdataset(file_name)['areacello'].to_masked_array() if var_name in ['co2mass']: # unit: ppmV flux,mass = False, False if var_name in ['co2', 'tas','co23D']: # unit: ppmV, degC, kg/kg #flux,mass = False, False weighted_mean = True file_name = glob.glob(dirname + model_name + '/areacella*.nc')[0] cell_area = xr.open_mfdataset(file_name)['areacella'].to_masked_array() ####################### Var_annual = [] for year in range(yrstart,2100): # 2015 0:12 # 2016 12:24 istart = (year-yrstart) * 12 # if (var_name=='cLand' and model_name=='NorESM2'): # # NorESM cLand starts from 201502 # if year==yrstart: # continue # else: # istart = istart -1 var = ds[var_name][istart:istart+12].to_masked_array() ################# if var_name in ['tas']: var = var - 273.15 if var_name in ['co2mass']: Var_annual.append(var[-1]*0.272912/(2.124*10**12)) if flux: global_sum = np.ma.sum(cell_area * var,axis=(1,2)) # global sum if year%4 == 0: if year%100 == 0: if year%400 == 0: days_in_Feb = 29 else: days_in_Feb = 28 else: days_in_Feb = 29 else: days_in_Feb = 28 print(global_sum.shape,'001') global_sum = (global_sum[0]*31 + global_sum[1]*days_in_Feb + global_sum[2]*31+ global_sum[3]*30 + global_sum[4]*31 + global_sum[5]*30+ global_sum[6]*31 + global_sum[7]*31 + global_sum[8]*30+ global_sum[9]*31 + global_sum[10]*30 + global_sum[11]*31) global_sum = global_sum / 10**12 * 86400 # kg/s to Pg/yr Var_annual.append(global_sum) if mass: global_sum = np.ma.sum(cell_area * var[-1,:,:]) if var_name in ['cLand']: global_sum = global_sum / 10**12 # kg to Pg if var_name in ['treeFrac']: global_sum = global_sum / 100. / 10000. / 1000000. # percentage, and then m2 to ha, and then to Mha Var_annual.append(global_sum) if weighted_mean: if var_name in ['co2','co23D']: var = var[:,0,:,:] annual = np.average(var,axis=0) cell_area.mask = annual.mask global_sum = np.ma.sum(cell_area * annual )/np.ma.sum(cell_area) if var_name in ['co2']: global_sum = global_sum * 1000000 if var_name in ['co23D']: global_sum = global_sum * 658244.5 Var_annual.append(global_sum) #print(f'global_sum {global_sum:.3f}') ###### #weights = np.cos(np.deg2rad(ds['lat'])) #weights.name = "weights" #var1 = ds[var_name][istart:istart+12] - 273.15 #print(var1.dims) #var_weighted = var1.weighted(weights) #print(var_weighted) #weighted_mean = var_weighted.mean(("lon", "lat")) #print(weighted_mean.mean().values) ######## if flux: if year %10==0: print(f'{year} {model_name} {exp_name} {var_name} {global_sum:.3f} Pg/yr, {sum(Var_annual):.3f} Pg accumulated') if mass: if year %10==0 and var_name in ['cLand']: print(f'{year} {model_name} {exp_name} {var_name} {global_sum:.3f} Pg, {global_sum-Var_annual[0]:.3f} Pg diff compared with start') if year %10==0 and var_name in ['treeFrac']: print(f'{year} {model_name} {exp_name} {var_name} {global_sum:.3f} Mha, {global_sum-Var_annual[0]:.3f} Mha diff compared with start') Var_annual = np.array(Var_annual) ##################### if var_name in ['co2mass','co2','co23D']: var_name = 'co2_concentration' if var_name =='treeFrac': var_name = 'tree' ##################### # save to csv if save_to_csv: df = pd.DataFrame(data=Var_annual, index=range(yrstart,2100), columns=[var_name]) df.index.name = 'year' #df.to_csv('LUMIP_csv/' + model_name + '_' + exp_name + '_' + var_name + '.csv') if real != '': df.to_csv('CDRMIP_csv/tt' + model_name + '_' + exp_name + '_' + var_name + '_' + real + '.csv') else: # if real == 'r1' df.to_csv('CDRMIP_csv/tt' + model_name + '_' + exp_name + '_' + var_name + '.csv') #df = pd.read_csv(expname + '.csv',index_col='year') #df[var_name] = Var_timeseries #df.to_csv(expname+'.csv') print(df)