#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Tue May 17 11:47:33 EDT 2022 if __name__ == '__main__': import sys from misc.timer import Timer tt = Timer('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 import xlinregress # if __name__ == '__main__': tt.check('end import') # #start from here ifile_T = sorted( glob.glob(f't_surf_FLOR_CTL1990_v201905_2xCO2_tigercpu_intelmpi_18_576PE_0101-????_glbmean.nc') )[-1] ifile_N = sorted( glob.glob(f'netrad_toa_FLOR_CTL1990_v201905_2xCO2_tigercpu_intelmpi_18_576PE_0101-????_glbmean.nc') )[-1] ifile_T_ctl = 't_surf_FLOR_CTL1990_v201905_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc' ifile_N_ctl = 'netrad_toa_FLOR_CTL1990_v201905_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc' n_iyears = 50# 15 n_years_fit = 200 da_T = xr.open_dataarray(ifile_T).groupby('time.year').mean('time') \ - xr.open_dataarray(ifile_T_ctl).sel(time=slice('0101', '0150')).mean('time') da_N = xr.open_dataarray(ifile_N).groupby('time.year').mean('time') \ - xr.open_dataarray(ifile_N_ctl).sel(time=slice('0101', '0150')).mean('time') ds = xr.Dataset(dict(TT=da_T.assign_attrs(units='K'), NN= da_N.assign_attrs(units='W/m^2'))) ds = ds.isel(year=slice(None, n_years_fit)).assign_coords(year=np.arange(1, n_years_fit+1)) #use 150 years ds0 = ds.isel(year=slice(None, n_iyears)) ds1 = ds.isel(year=slice(n_iyears, None)) #regression of slow process years rg = ds1.NN.linregress.on(ds1.TT, dim='year') N0 = rg.intercept.item() T0 = -rg.intercept.item()/rg.slope.item() #regression of initial years (fast process) rg0 = ds0.NN.linregress.on(ds0.TT, dim='year') N00 = rg0.intercept.item() T00 = -rg0.intercept.item()/rg0.slope.item() #regression of all years rg_all = ds.NN.linregress.on(ds.TT, dim='year') N0_all = rg_all.intercept.item() T0_all = -rg_all.intercept.item()/rg_all.slope.item() def cal_sensitivity(T, N, nyears=None, nyears_removed=None): """calculate climate sensitivity given T/N time series and # of years to be used.""" if nyears is not None: T = T.isel(year=slice(None, nyears)) N = N.isel(year=slice(None, nyears)) if nyears_removed is not None: T = T.isel(year=slice(nyears_removed, None)) N = N.isel(year=slice(nyears_removed, None)) rg = N.linregress.on(T, dim='year') F = rg.intercept.item() ECS = -rg.intercept.item()/rg.slope.item() return ECS, F #ECS/F as function of # of years used nyears_vec = np.arange(10,n_years_fit+1) TN0 = np.array([cal_sensitivity(ds.TT, ds.NN, n, None) for n in nyears_vec]) ECS = xr.DataArray(TN0[:,0], dims=['year'], coords=[nyears_vec,]).assign_attrs(units='K') F = xr.DataArray(TN0[:,1], dims=['year'], coords=[nyears_vec,]).assign_attrs(units='W/m^2') #as function of nyears_removed nyears_removed_vec = np.arange(0,150+1) TN0 = np.array([cal_sensitivity(ds.TT, ds.NN, n_years_fit, n) for n in nyears_removed_vec]) ECS_ = xr.DataArray(TN0[:,0], dims=['year'], coords=[nyears_removed_vec,]).assign_attrs(units='K') F_ = xr.DataArray(TN0[:,1], dims=['year'], coords=[nyears_removed_vec,]).assign_attrs(units='W/m^2') if __name__ == '__main__': from wyconfig import * #my plot settings #time series fig, axes = plt.subplots(2, 2, figsize=(10,6)) #fig, ax = plt.subplots() ax = axes[0,0] ds.NN.plot(ax=ax) ax.set_ylabel('N [W m$^{-2}$]', color='C0') ax.grid(False) ax1 = ax.twinx() ds.TT.plot(ax=ax1, color='C1') ax1.set_ylabel('T [K]', color='C1') ax1.grid(False) ax1.spines['right'].set_visible(True) ax1.axvline(n_iyears, color='gray', ls='--') #climate forcing and sensitivity #fig, ax = plt.subplots() ax = axes[0,1] #scatter plot alpha = 0.5 ds0.plot.scatter(ax=ax, x='TT', y='NN', color='none', edgecolor='C2', label=f'years 1-{n_iyears}', alpha=alpha) ds1.plot.scatter(ax=ax, x='TT', y='NN', color='none', edgecolor='C3', label=f'years {n_iyears+1}-{n_years_fit}', alpha=alpha) #regression of slow process years ax.axline((0,N0), (T0,0), color='C3', ls='--') ax.text(0, N0, f' F={N0:.2f}W m$^{{-2}}$', ha='left', va='top') ax.text(T0, 0, f' ECS={T0:.2f}K', ha='center', va='bottom') #regression of initial years ax.axline((0,N00), (T00,0), color='C2', ls='--') ax.text(0, N00, f' F={N00:.2f}W m$^{{-2}}$', ha='left') ax.text(T00, 0, f' ECS={T00:.2f}K', ha='right', va='bottom') #regression of all years ax.axline((0,N0_all), (T0_all,0), color='gray', ls='--') ax.text(0, N0_all, f' F={N0_all:.2f}W m$^{{-2}}$', ha='left') ax.text(T0_all, 0, f' ECS={T0_all:.2f}K', ha='right', va='top') ax.legend() ax.set_ylabel('N [W m$^{-2}$]') ax.set_xlabel('T [K]') #climate sensitivity/forcing as function of # of years used #fig, ax = plt.subplots() ax = axes[1,0] F.plot(ax=ax) ax.set_ylabel('F [W m$^{-2}$]', color='C0') ax.set_xlabel('# of years used in ECS calculation') ax.grid(False) ax.axhline(3.66, color='C0', ls='--', label='F from AGCM exp') ax.legend(loc='center left') ax1 = ax.twinx() ECS.plot(ax=ax1, color='C1') ax1.set_ylabel('ECS [K]', color='C1') ax1.grid(False) ax1.spines['right'].set_visible(True) ax1.axvline(n_iyears, color='gray', ls='--') #climate sensitivity/forcing as function of # of years removed #fig, ax = plt.subplots() ax = axes[1,1] F_.plot(ax=ax) ax.set_ylabel('F [W m$^{-2}$]', color='C0') ax.set_xlabel(f'# of years removed in ECS calculation ({n_years_fit} years in total)') ax.grid(False) ax.axhline(3.66, color='C0', ls='--', label='F from AGCM exp') ax.legend(loc='center left') ax1 = ax.twinx() ECS_.plot(ax=ax1, color='C1') ax1.set_ylabel('ECS [K]', color='C1') ax1.grid(False) ax1.spines['right'].set_visible(True) ax1.axvline(n_iyears, color='gray', ls='--') #savefig if len(sys.argv)>1 and 'savefig' in sys.argv[1:]: figname = __file__.replace('.py', f'.png') if 'overwritefig' in sys.argv[1:]: wysavefig(figname, overwritefig=True) else: wysavefig(figname) tt.check(f'**Done**') print() if 'notshowfig' in sys.argv: pass else: plt.show()