''' calculate global values from JSBACH gridded output ''' import pandas as pd import copy import numpy as np import matplotlib.pyplot as plt import netCDF4 as nc import sys ################# var_name_list = ['temp2_JJA','soil_temperature','sst','co2_flx_land','co2_flx_ocean','temp2','nep','cLand','cProduct','nbp','gpp','ra','npp','rh','fire','harvest','lcc', 'crop_harvest','cLitter'] iplotvar = int(sys.argv[1]) var_name = var_name_list[iplotvar-1] print(f'plotting variable: {var_name}') #exp = 'hist' #exp = 'OAE_STOP' #exp = 'OAE' #exp = 'AFF' #exp = 'SynTraAR' #exp = 'SynTraOAE' #exp = 'SynTraBOTH' #exp = 'SynTraHalfOAE' exp = 'SynTraHalfBOTH' #exp = 'SynTraHalfAR' #exp = 'VC003' #exp = 'constLu' #exp = 'CC021' #exp = 'ssp126' gridarea_T63 = nc.Dataset('/home/shkhwwey/OAE/gridarea_T63.nc').variables['cell_area'][:] endyear = 2100 #if exp=='OAE': # expname = ['CC104','CC105','CC106','HW002','HW005','HW006'] # startyear = 2015 #if exp=='AFF': # expname = ['CC104','CC105','CC106','HW003','HW007','HW008'] # startyear = 2015 #if exp=='OAE_STOP': # expname =['HW004','HW009','HW010'] # startyear = 2070 #if exp in['ssp126']: # expname =['NM004'] # startyear = 2015 #if exp=='CC021': # expname=['CC021'] # startyear = 1850 # endyear = 1925 #if exp=='VC003': # expname = [exp] # startyear = 2015 flux = False mass = False save_to_csv = True model = 'jsbach' dirlist = [] if exp =='OAE': 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/') if exp =='AFF': 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-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/') if exp =='OAE_STOP': 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 exp =='ssp126': dirlist.append ('/scratch/usr/shkhwwey/output_others/FOCI2.0-NM004_FM_ESM_SSP126os_2014CC105/nc/') if exp=='constLu': expname = ['HW013','HW012','HW011'] startyear = 2015 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 exp=='SynTraOAE': expname = ['HW028','HW032','HW033'] dirlist.append ('/scratch/usr/shkhwwey/models/foci2.0_maps_wh/experiments/FOCI20.0-HW028_esm_ssp370_OAE/nc/') 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/') if exp=='SynTraAR': expname = ['HW022','HW023','HW024'] startyear = 2015 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 exp=='SynTraBOTH': expname = ['HW029','HW030','HW031'] startyear = 2015 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-HW030_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 exp=='SynTraHalfOAE': expname = ['HW034','HW035','HW036'] startyear = 2015 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 exp=='SynTraHalfBOTH': #expname = ['HW037','HW038','HW039'] expname = ['HW037','HW038'] #expname = ['HW039'] startyear = 2015 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 exp=='SynTraHalfAR': expname = ['HW040','HW041','HW042'] startyear = 2015 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/') #### if exp=='CC021': #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-CC021_1pctCO2_spinup2329/nc/') if exp=='VC003': dirlist.append('/scratch/usr/shkhwwey/output_others/FOCI2.0-VS003_ssp534os_2099/nc/') if exp=='hist': expname = ['CC104_hist','CC105_hist','CC106_hist'] startyear,endyear = 1850,2014 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/') ######### print(exp) print(expname) print(dirlist) ####### # (nexp, nyear) Var_timeseries = [] Var_accu_timeseries = [] for i in range(len(dirlist)): Var_timeseries.append([]) Var_accu_timeseries.append([]) ######### molC_to_kgC = 0.01201 kgCO2_to_kgC = 0.272912 ### #cLand_2014_arr = nc.Dataset(dirlist[0]+'veg_2014.nc').variables['box_Cpools_total'][:] * molC_to_kgC #cLand_2014 = np.ma.sum(gridarea_T63 * cLand_2014_arr[-1]) /10**12 ## for year in range(startyear,endyear): monthly_by_exp = [] 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/': # # there is no 2067 for HW037 # monthly_by_exp.append(np.ones(12)*-999) # continue if var_name =='temp2_JJA': model = 'echam' temp2_JJA = nc.Dataset(dirname+'echam6_echam_' + str(year) + '_JJA.nc').variables['temp2'][:] -273.15 monthly_by_exp.append(temp2_JJA) slm = 1 if var_name =='temp2': model = 'echam' temp2 = nc.Dataset(dirname+'echam6_echam_' + str(year) + '.nc').variables['temp2'][:] -273.15 monthly_by_exp.append(temp2) slm = 1 if var_name =='sst': slm = nc.Dataset('/home/shkhwwey/OAE/slm.nc').variables['slm'][:] # sea=0; land=1 slm = slm -1 * -1 # now sea=1, land=0 model = 'echam' sst = nc.Dataset(dirname+'echam6_echam_temp_' + str(year) + '.nc').variables['tsurf'][:] -273.15 sst = sst * slm monthly_by_exp.append(sst) if var_name =='co2_flx_land': co2_flx_land = nc.Dataset(dirname + 'echam6_co2_' + str(year) + '.nc').variables['co2_flx_land'][:] * kgCO2_to_kgC * -1. monthly_by_exp.append(co2_flx_land) flux = True if var_name =='co2_flx_ocean': co2_flx_ocean = nc.Dataset(dirname + 'echam6_co2_' + str(year) + '.nc').variables['co2_flx_ocean'][:] * kgCO2_to_kgC * -1. monthly_by_exp.append(co2_flx_ocean) flux = True if var_name =='soil_temperature': soil_temperature = nc.Dataset(dirname + 'jsbach_land_' + str(year) + '.nc').variables['soil_temperature'][:,0] - 273.15 monthly_by_exp.append(soil_temperature) model = 'echam' #slm = nc.Dataset('/home/shkhwwey/OAE/slm.nc').variables['slm'][:] # sea=0; land=1 slm = nc.Dataset('/home/shkhwwey/OAE/slm_ice_free.nc').variables['slm'][:] # sea=0; land=1 if model=='jsbach' and var_name not in ['co2_flx_ocean','co2_flx_land']: # nbp = gpp - ra - rh - fFire - fHarvest - fLuc (fGrazing is included in rh) # rh = rs + herb - crop harvest (rs output includes crop harvest, but after cmor rs should not include crop harvest) gpp = np.sum (nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_GPP_yDayMean'][:],axis=1) * molC_to_kgC npp = np.sum (nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_NPP_act_yDayMean'][:],axis=1) * molC_to_kgC ra = gpp - npp rs = np.sum (nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_soil_respiration'][:],axis=1) * molC_to_kgC * -1. herb = np.sum (nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_Cflux_herbivory_2_atm'][:],axis=1) * molC_to_kgC * -1. crop_harvest = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_Cflx_crop_harvest_2_atm'][:] * molC_to_kgC fire = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_fire_CO2_flux_2_atmos'][:] * kgCO2_to_kgC harvest = nc.Dataset(dirname+'jsbach_' + str(year) + '.nc').variables['CO2_emission_harvest'][:] * molC_to_kgC lcc = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['LCC_flux_box_C2atmos'][:] * molC_to_kgC rh = rs + herb - crop_harvest nep = npp - rh - fire nbp = npp - rh - fire - harvest - lcc - crop_harvest cLand = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['box_Cpools_total'][:] * molC_to_kgC if var_name =='ra': monthly_by_exp.append(ra) flux = True if var_name =='gpp': monthly_by_exp.append(gpp) flux = True if var_name =='cProduct': var218 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_onSite_LCC'][:] var220 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_paper_LCC'][:] var221 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_construction_LCC'][:] var222 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_paper_harvest'][:] var223 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_construction_harvest'][:] var224 = nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_onSite_harvest'][:] var225 = np.sum(nc.Dataset(dirname+'jsbach_cProduct_'+str(year)+'.nc').variables['boxC_crop_harvest'][:],axis=1) cProduct = (var218+var220+var221+var222+var223+var224+var225) * molC_to_kgC monthly_by_exp.append(cProduct) mass = True if var_name =='cLitter': var31 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_acid_ag1'][:] var33 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_water_ag1'][:] var35 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_ethanol_ag1'][:] var37 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_nonsoluble_ag1'][:] var41 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_acid_ag2'][:] var43 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_water_ag2'][:] var45 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_ethanol_ag2'][:] var47 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_nonsoluble_ag2'][:] cLitter = (var31+var33+var35+var37+var41+var43+var45+var47) * molC_to_kgC monthly_by_exp.append(cLitter) mass = True if var_name =='cSoil': var32 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_acid_bg1'][:] var34 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_water_bg1'][:] var36 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_ethanol_bg1'][:] var38 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_nonsoluble_bg1'][:] var39 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_humus_1'][:] var42 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_acid_bg2'][:] var44 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_water_bg2'][:] var46 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_ethanol_bg2'][:] var48 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_nonsoluble_bg2'][:] var49 = nc.Dataset(dirname+'jsbach_yasso_'+str(year)+'.nc').variables['boxYC_humus_2'][:] cSoil = (var32+var34+var36+var38+var39+var42+var44+var46+var48+var49) * molC_to_kgC monthly_by_exp.append(cSoil) mass = True if var_name =='cVeg': cGreen = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['boxC_green'][:] * molC_to_kgC cWoods = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['boxC_woods'][:] * molC_to_kgC cReserve = nc.Dataset(dirname+'veg_' + str(year) + '.nc').variables['boxC_reserve'][:] * molC_to_kgC cVeg = cGreen + cWoods + cReserve monthly_by_exp.append(cVeg) mass = True if var_name =='fire': monthly_by_exp.append(fire) flux = True #if var_name =='rs': # monthly_by_exp.append(rs) # flux = True if var_name =='rh': monthly_by_exp.append(rh) flux = True if var_name =='npp': monthly_by_exp.append(npp) flux = True if var_name =='nep': monthly_by_exp.append(nep) flux = True if var_name =='lcc': monthly_by_exp.append(lcc) flux = True if var_name =='harvest': monthly_by_exp.append(harvest) flux = True if var_name =='nbp': monthly_by_exp.append(nbp) flux = True if var_name =='crop_harvest': monthly_by_exp.append(crop_harvest) flux = True if var_name =='cLand': monthly_by_exp.append(cLand) mass = True print(f'---\n{year}') ################# # calculate annual value for iexp in range(len(monthly_by_exp)): if model =='echam': print(monthly_by_exp[iexp].shape) annual = np.average(monthly_by_exp[iexp],axis=0) global_sum = np.ma.sum(gridarea_T63 * annual * slm)/np.ma.sum(gridarea_T63 * slm) #print(f'{year} {global_sum:.2f}') if flux: #global_sum = np.ma.sum(gridarea_T63 * np.ma.average(monthly_by_exp[i],axis=0)) #global_sum = global_sum / 10**12 * 86400 * 365.25 global_sum = np.ma.sum(gridarea_T63 * monthly_by_exp[iexp],axis=(1,2)) 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 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 if mass: global_sum = np.ma.sum(gridarea_T63 * monthly_by_exp[iexp][-1]) global_sum = global_sum / 10**12 ######## Var_timeseries[iexp].append(global_sum) Var_accu_timeseries[iexp].append( sum(Var_timeseries[iexp]) ) # if flux: # print('{} {:.3f} Pg/yr, {:.3f} Pg accumulated'.format(expname[iexp],global_sum,sum(Var_timeseries[iexp]))) # if mass: # if len(Var_timeseries[iexp])>=2: # print('{} {:.3f} Pg, {:.3f} Pg diff compared with start, {:.3f} Pg diff compared with last year'.format(expname[iexp],global_sum,global_sum-cLand_2014,global_sum-Var_timeseries[iexp][-2])) Var_timeseries = np.array(Var_timeseries) Var_accu_timeseries = np.array(Var_accu_timeseries) Var_accu_timeseries = Var_accu_timeseries ####################### # save to csv if save_to_csv: for iexp in range(len(Var_timeseries)): # combine variables #df = pd.read_csv('csv/' + expname[iexp] + '.csv',index_col='year') #df[var_name] = Var_timeseries[iexp] #df.to_csv(expname[iexp]+'.csv') #print(df) # separate file for the current variables df = pd.DataFrame(data=Var_timeseries[iexp], index=range(startyear,endyear), columns=[var_name]) df.index.name = 'year' df.to_csv('temp_csv/' + expname[iexp] + '_' + var_name + '.csv') ######### # parameters for figures fontsize = 18 collist = ['C0','C1','C2','C3','C4','C5'] collist = ['C0','C0','C0','C1','C1','C1'] ##### # plot timeseries of all experiments fig = plt.figure() ax = fig.add_subplot(111) for i in range(len(Var_timeseries)): print(f'{expname[i]} {np.average(Var_timeseries[i,:]):.2f}') #if i==1: #continue ax.plot(range(startyear,endyear),Var_timeseries[i,:],'.-',color=collist[i],label = expname[i]) ax.grid() ax.legend(loc = 'best',fontsize=fontsize) if flux: ax.set_ylabel('PgC/yr',fontsize=fontsize) if mass: ax.set_ylabel('PgC',fontsize=fontsize) if var_name=='temp2': ax.set_ylabel('degC',fontsize=fontsize) ax.tick_params(labelsize=fontsize) if 'long_name' in locals(): fig.suptitle(long_name,fontsize=18) else: fig.suptitle(var_name,fontsize=18) #plt.show() plt.savefig('fig/FOCI_' + var_name + '.png',bbox_inches='tight') plt.close() if model=='echam': exit() ########## # plot timeseries of all experiments (accu.) fig = plt.figure() ax = fig.add_subplot(111) for i in range(len(Var_accu_timeseries)): if i==1: continue ax.plot(range(startyear,endyear),Var_accu_timeseries[i,:],'.-',color=collist[i],label = expname[i]) ax.grid() ax.legend(loc = 'best', fontsize=fontsize) ax.set_ylabel('PgC',fontsize=fontsize) ax.tick_params(labelsize=fontsize) if 'long_name' in locals(): fig.suptitle(long_name,fontsize=18) else: fig.suptitle(var_name,fontsize=18) #plt.show() plt.savefig('fig/FOCI_' + var_name + '_accu.png',bbox_inches='tight') plt.close() print(var_name)