#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Wed Jul 5 12:03:56 EDT 2023 # run this script after finishing FigS01.py and FigS02.py if __name__ == '__main__': import sys,os from misc.timer import Timer tt = Timer(f'[{os.getcwd()}] start ' + ' '.join(sys.argv)) import sys, os.path, os, glob, datetime import xarray as xr, numpy as np, pandas as pd, matplotlib.pyplot as plt #more imports wython = '/tigress/wenchang/wython' if wython not in sys.path: sys.path.append(wython); print('added to python path:', wython) #from misc import get_kws_from_argv # if __name__ == '__main__': tt.check('end import') # #start from here forcings = ['co2x2', 'co2x4', 'co2xp5', 'co2xp25', 'p2p0solar', 'p6p0solar', 'm2p0solar', 'm6p0solar']#[:4] labels = ['2xCO2', '4xCO2', '0.5xCO2', '0.25xCO2', '$+2$% solar', '$+6$% solar', '$-2$% solar', '$-6$% solar']#[:4] forcingScales = [1, 2, -1, -2, 1, 3, -1, -3]#[:4] yearspan = slice(231,250) ofile = __file__.replace('.py', f'__{yearspan.start-100}-{yearspan.stop-100}.nc') # simulations start from year 101. make them start from year 1 if not os.path.exists(ofile) or 'od' in sys.argv: das_eta_cmip = [] #eta from coupled model experiments das_eta_approx = [] #eta from coupled model experiments estimated from AGCM experiments das_eta_cmip_tmsst = [] #eta defined as sensitivity of global mean precip to tropical mean sst for forcing,label,scale in zip(forcings,labels,forcingScales): #forcing, label, scale = forcings[0], labels[0], forcingScales[0] #prediction by gmst ifile = f'fig_lines_precip__FLOR_{forcing}_glbmeanTs.nc' print(ifile) ds = xr.open_dataset(ifile) #adjustment term if 'co2' in forcing: A = ds['A_co2'] else: A = ds['A_solar'] A = A*scale da_actual = ds['pr'].groupby('time.year').mean('time')/ds['pr_ctl_clim'].mean()*100 -100 #% of precip change dTs = ds['Ts'].groupby('time.year').mean('time') - ds['Ts_ctl_clim'].mean() dTs_gmst = dTs #da_predict = dTs * ds['eta'] + A eta_cmip = (da_actual.sel(year=yearspan).mean('year') - A)/dTs.sel(year=yearspan).mean('year') print(f'eta_cmip = {eta_cmip.item()}') #eta_cmip = (da_actual - A)/dTs das_eta_cmip.append(eta_cmip) #prediction by tropical mean SST ifile = f'fig_lines_precip__FLOR_{forcing}_tropoceanmeanTs.nc' ds = xr.open_dataset(ifile) da_actual = ds['pr'].groupby('time.year').mean('time')/ds['pr_ctl_clim'].mean()*100 -100 #% of precip change dTs = ds['Ts'].groupby('time.year').mean('time') - ds['Ts_ctl_clim'].mean() #da_predict_to = dTs * ds['eta'] + A #eta_cmip_tmsst = (da_actual - A)/dTs eta_cmip_tmsst = (da_actual.sel(year=yearspan).mean('year') - A)/dTs.sel(year=yearspan).mean('year') print(f'eta_cmip_tmsst = {eta_cmip_tmsst.item()}') das_eta_cmip_tmsst.append(eta_cmip_tmsst) # dTs_tmsst = dTs eta_approx = ds['eta'] * dTs_tmsst.sel(year=yearspan).mean('year')/dTs_gmst.sel(year=yearspan).mean('year') print(f'eta_approx = {eta_approx.item()}') #eta_approx = ds['eta'] * dTs_tmsst/dTs_gmst das_eta_approx.append(eta_approx) da_eta_cmip = xr.concat(das_eta_cmip, dim=pd.Index(labels, name='label')) da_eta_cmip_tmsst = xr.concat(das_eta_cmip_tmsst, dim=pd.Index(labels, name='label')) da_eta_approx = xr.concat(das_eta_approx, dim=pd.Index(labels, name='label')) ds = xr.Dataset(dict(eta_cmip=da_eta_cmip,eta_cmip_tmsst=da_eta_cmip_tmsst, eta_approx=da_eta_approx)) #save ds.to_netcdf(ofile) print('[saved]:', ofile, 'see', __file__) else: ds = xr.open_dataset(ofile) print('[loaded]:', ofile, 'see', __file__) def wyplot(ax=None, eta_type='gmsst'): if ax is None: fig,ax = plt.subplots() if eta_type == 'gmsst': #df = ds[['eta_cmip', 'eta_approx']].rename(eta_cmip='$\eta$ in FLOR experiments', eta_approx='equivalent $\eta$').to_pandas() df = ds[['eta_cmip']].rename(eta_cmip='FLOR experiments').to_pandas() else: #eta_tpye==tmsst df = ds[['eta_cmip_tmsst']].rename(eta_cmip_tmsst='FLOR experiments').to_pandas() #print(df) df.plot.bar(rot=15, ax=ax, color='C0') if eta_type == 'gmsst': ax.axhline(3.28, color='C1', ls='--', label='AGCM uniform warming ($\eta_{+2K}$)') else: #eta_type == 'tmsst' ax.axhline(3.52, color='C1', ls='--', label='AGCM uniform warming ($\eta^{TMSST}_{+2K}$)') ax.legend(loc='upper right', ncol=2, borderpad=0) #ax.legend(loc='upper right', ncol=1, borderpad=0) ax.set_ylim(None, 4) ax.set_xlabel('FLOR experiments') ax.set_ylabel('%/K') if eta_type=='gmsst': ax.set_title(f'$\eta$') else: #eta_type=='tmsst' ax.set_title(f'$\eta^{{TMSST}}$') def wyplot_gap_explained(ax=None): if ax is None: fig,ax = plt.subplots() eta_u = 3.28 #eta in uniform warming da = (eta_u - ds['eta_approx']) /(eta_u - ds['eta_cmip']) #* 100 df = da.to_pandas() df.plot.bar(rot=15, ax=ax, color='C2') ax.set_xlabel('FLOR experiments') #ax.set_ylabel('($\eta_{+2K} - \eta_{e}$)/($\eta_{+2K} - \eta_{FLOR} \\times 100\%$') ax.set_ylabel('($\eta_{+2K} - \eta_{e}$)/($\eta_{+2K} - \eta_{FLOR}$)') if __name__ == '__main__': from wyconfig import * #my plot settings wyplot() #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'.png') if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) wyplot(eta_type='tmsst') #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'__tmsst.png') if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) #wyplot_gap_explained() tt.check(f'**Done**') print() if 'notshowfig' in sys.argv or 'n' in sys.argv: pass else: if 'plt' in globals(): plt.show()