import my_functions import copy import numpy as np import matplotlib.pyplot as plt import netCDF4 as nc import pandas as pd ####### # molar Mass CO2: 44.0095 # molar Mass C: 12.01 # => factor: 0.272912 ################# save_to_csv = True gridarea_T63 = nc.Dataset('/home/shkhwwey/OAE/gridarea_T63.nc').variables['cell_area'][:] dirlist = [] ########## #expset = 'SynTraAE' #expset = 'SynTraOAE' #expset = 'SynTraBOTH' #expset = 'SynTraHalfOAE' #expset = 'SynTraHalfAR' expset = 'SynTraHalfBOTH' #expname = ['CC104','CC105','CC106','HW002','HW005','HW006'] #expname = ['CC104','CC105','CC106','HW003','HW007','HW008'] #expname = ['HW004','HW009','HW010'] #expname = ['CC104'] #dirlist.append ('/scratch/usr/shkhwwey/output_others/FOCI1.20.0-CC104_RCP_ESM_spinup2089/nc/') #dirlist.append ('/scratch/usr/shkhwwey/output_others/FOCI1.20.0-CC105_RCP_ESM_spinup2099/nc/') #dirlist.append ('/scratch/usr/shkhwwey/output_others/FOCI1.20.0-CC106_RCP_ESM_spinup2079/nc/') ##dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW002_esm_ssp585_ocn_alk/nc/') ##dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW005_esm_ssp585_ocn_alk/nc/') ##dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW006_esm_ssp585_ocn_alk/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW003_esm_ssp585_ssp126Lu/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW007_esm_ssp585_ssp126Lu/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW008_esm_ssp585_ssp126Lu/nc/') ##dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW004_esm_ssp585_ocn_alk_stop/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW004_esm_ssp585_ocn_alk_stop/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW009_esm_ssp585_ocn_alk_stop/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW010_esm_ssp585_ocn_alk_stop/nc/') ################### if expset == 'SynTraREF': expname = ['HW013','HW012','HW011'] yrstart, yrend = 2015, 2100 dirlist.append('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW013_esm_ssp370_constLu/nc/') dirlist.append('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW012_esm_ssp370_constLu/nc/') dirlist.append('/scratch/usr/shkhwwey/models/foci2.0/experiments/FOCI20.0-HW011_esm_ssp370_constLu/nc/') if expset == 'SynTraAR': expname = ['HW022','HW023','HW024'] yrstart, yrend = 2015, 2100 dirlist.append('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW022_esm_ssp370_AR/nc/') dirlist.append('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW023_esm_ssp370_AR/nc/') dirlist.append('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW024_esm_ssp370_AR/nc/') if expset=='SynTraOAE': expname = ['HW032','HW033'] dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW032_esm_ssp370_OAE/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW033_esm_ssp370_OAE/nc/') #expname = ['HW026','HW027','HW028'] yrstart, yrend = 2015, 2100 #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW026_esm_ssp370_OAE/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW027_esm_ssp370_OAE/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW028_esm_ssp370_OAE/nc/') if expset =='SynTraAR_OAE': expname = ['HW030'] yrstart, yrend = 2015, 2100 dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW030_esm_ssp370_AR_and_OAE/nc/') #expname = ['HW029','HW031'] #yrstart, yrend = 2015, 2100 #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW029_esm_ssp370_AR_and_OAE/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW031_esm_ssp370_AR_and_OAE/nc/') if expset =='SynTraAE': expname = ['HW044'] yrstart, yrend = 2015, 2081 dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW044_esm_ssp370_OAE_avoided_emission/nc/') if expset=='SynTraHalfOAE': expname = ['HW034','HW035','HW036'] yrstart, yrend = 2015, 2100 dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW034_esm_ssp370_half_OAE/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW035_esm_ssp370_half_OAE/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW036_esm_ssp370_half_OAE/nc/') if expset=='SynTraHalfBOTH': expname = ['HW037','HW038'] #expname = ['HW039'] yrstart, yrend = 2015, 2100 dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW037_esm_ssp370_half_AR_and_OAE/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW038_esm_ssp370_half_AR_and_OAE/nc/') #dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW039_esm_ssp370_half_AR_and_OAE/nc/') if expset=='SynTraHalfAR': expname = ['HW040','HW041','HW042'] yrstart, yrend = 2015, 2100 dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW040_esm_ssp370_half_AR/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW041_esm_ssp370_half_AR/nc/') dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW042_esm_ssp370_half_AR/nc/') ################## co2_timeseries = [] for i in range(len(dirlist)): co2_timeseries.append([]) #for year in range(2015,2100): #for year in range(2070,2100): for year in range(yrstart,yrend): co2list = [] for dirname in dirlist: #if year==2067 and dirname=='/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW037_esm_ssp370_half_AR_and_OAE/nc/': # co2list.append(np.ones((12,96,192))*-999) # continue co2list.append (nc.Dataset(dirname+'echam6_tracer_' + str(year) + '.nc').variables['CO2'][:,94]) for i in range(len(co2list)): #if year==2067 and expset=='SynTraHalfBOTH' and i==0: # co2_timeseries[i].append(-999) # continue global_avg = np.ma.sum(gridarea_T63 * np.ma.average(co2list[i],axis=0))/np.ma.sum(gridarea_T63) # unit: kg/kg # convert CO2 from [kg/kg] to [ppmV]; # factor: 28.970 / 44.011 * 1000000 = 658244.5 global_avg = global_avg * 658244.5 co2_timeseries[i].append(global_avg) #print('{} {:.3f} ppm'.format(expname[i],global_avg)) co2_timeseries = np.array(co2_timeseries) print('co2_timeseries',co2_timeseries.shape,len(co2_timeseries)) #exit() ######### # plot timeseries of all experiments #fig = plt.figure() #ax = fig.add_subplot(111) #for i in range(len(Var_timeseries)): # ax.plot(range(2015,2100),Var_timeseries[i,:],'.-',label = expname[i]) #ax.grid() #ax.legend(loc = 'best') #ax.set_ylabel('PgC/yr') # #if 'long_name' in locals(): # fig.suptitle(long_name,fontsize=18) #else: # fig.suptitle(var_name,fontsize=18) ##plt.show() #plt.savefig('fig/' + var_name + '.png',bbox_inches='tight') #plt.close() ########## # plot timeseries of all experiments (accu.) #fig = plt.figure() #ax = fig.add_subplot(111) #for i in range(len(Var_accu_timeseries)): # ax.plot(range(2015,2100),Var_accu_timeseries[i,:],'.-',label = expname[i]) #ax.grid() #ax.legend(loc = 'best') #ax.set_ylabel('PgC') # #if 'long_name' in locals(): # fig.suptitle(long_name,fontsize=18) #else: # fig.suptitle(var_name,fontsize=18) ##plt.show() #plt.savefig('fig/' + var_name + '_accu.png',bbox_inches='tight') #plt.close() ######### # plot CO2 concentration ################ # save co2 in ppm to csv if save_to_csv: for iexp in range(len(co2_timeseries)): df = pd.DataFrame(data=co2_timeseries[iexp], index=range(yrstart,yrend), columns=['co2_concentration']) df.index.name = 'year' #print(df) df.to_csv('temp_csv/' + expname[iexp] + '_co2_concentration.csv') exit() ############ fig = plt.figure() ax = fig.add_subplot(111) for i in range(len(co2_timeseries)): print('{} {:.3f} ppm (start)'.format(expname[i],co2_timeseries[i,0])) print('{} {:.3f} ppm'.format(expname[i],co2_timeseries[i,-1])) ax.plot(range(1850,1860),co2_timeseries[i,:],'.-',label = expname[i]) #for i in range(len(co2_from_scenario)): # ax.plot(range(2015,2100),co2_from_scenario[i,2015:2100],'.-',label = scenario_list[i]) #ax.plot(range(2015,2090),co2_from_mon[2015-1850:2090-1850],'.-',label = 'monitor') #for yr in range(2015,2090): # print('yr {}, esm: {}, from_monitor: {}, from_input: {}'.format(yr,co2_timeseries[0,yr-2015],co2_from_mon[yr-1850],co2_from_scenario[3,yr])) ax.grid() ax.legend(loc = 'best') ax.set_ylabel('ppmV') fig.suptitle('CO2 concentration',fontsize=18) #plt.show() plt.savefig('fig/FOCI_co2_compare.png',bbox_inches='tight') exit() ############## # MPI-ESM print('===') varname = 'nep' filename = '/scratch/usr/shkhwwey/MPI-ESM/' + varname + '_Emon_MPI-ESM1-2-LR_esm-hist_r1i1p1f1_gn_201001-201412.nc' Var=nc.Dataset(filename).variables[varname][:] global_sum = np.ma.sum(gridarea_T63 * np.ma.average(Var[-12:],axis=0)) global_sum = global_sum / 10**12 * 86400 * 365.25 #* 12.01 / 44.0095 print('{} {} {:.3f} Pg/yr'.format('MPI-ESM',varname,global_sum))