#!/usr/bin/env python # Wenchang Yang (wenchang@princeton.edu) # Fri Mar 3 13:21:20 EST 2023 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 #from misc import get_kws_from_argv # if __name__ == '__main__': tt.check('end import') # #start from here model = 'FLOR' case = """ m1p7pctSolar 2xCO2_m1p7pctSolar 2xCO2 m2pctSolar 2xCO2_m2pctSolar """.split()[-1] ifiles = { '2xCO2': 'precip_FLOR_CTL1990_v201905_2xCO2_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc', 'm2pctSolar': 'precip_FLOR_CTL1990_v201905_m2pctSolar_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc', '2xCO2_m2pctSolar': 'precip_FLOR_CTL1990_v201905_2xCO2_m2pctSolar_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc', 'm1p7pctSolar': 'precip_FLOR_CTL1990_v201905_m1p7pctSolar_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc', '2xCO2_m1p7pctSolar': 'precip_FLOR_CTL1990_v201905_2xCO2_m1p7pctSolar_tigercpu_intelmpi_18_576PE_0101-0300_glbmean.nc', } #nonlinear_ratio = True b_gmst_am2p1 = 2.63 #2.628766826842829 # 2.6 # % per K b_fco2_am2p1 = -2.40 #-2.403719774284724 #-2.4 # % per 2xCO2 forcing #AM2.5 values: see /tigress/wenchang/analysis/misc/AM2.5/wyshow_dprecip.log b_gmst_am2p5 = 2.98 # % per K b_fco2_am2p5 = -2.27 # % per 2xCO2 forcing #HIRAM values: see /tigress/wenchang/analysis/misc/HIRAM/wyshow_dprecip.log b_gmst_hiram = 2.97 # % per K b_fco2_hiram = -2.32 # % per 2xCO2 forcing #AM2.5C360 values: see /tigress/wenchang/analysis/misc/AM2.5/wyshow_dprecip.log b_gmst_am2p5c360 = 3.29 # % per K b_fco2_am2p5c360 = -2.23 # % per 2xCO2 forcing #select the params b_gmst = b_gmst_am2p5 #b_gmst = 3.25 #%/K b_fco2 = b_fco2_am2p5 b_solar = -0.83 # % per 2% solar #precip ctl #ifile = 'precip_FLOR_CTL1860_newdiag_tigercpu_intelmpi_18_576PE_0001-2301_glbmean.nc' ifile = 'precip_FLOR_CTL1990_v201905_tigercpu_intelmpi_18_576PE_0001-1393_glbmean.nc' da = xr.open_dataarray(ifile) da_ctl = da da_ref = da.sel(time=slice('0101', '0200')).mean('time') #precip pert ifile = ifiles[case] print(ifile) da = xr.open_dataarray(ifile) da_precip = da #Ts ctl #ifile = 't_surf_FLOR_CTL1860_newdiag_tigercpu_intelmpi_18_576PE_0001-2301_glbmean.nc' ifile = 't_surf_FLOR_CTL1990_v201905_tigercpu_intelmpi_18_576PE_0001-1393_glbmean.nc' da = xr.open_dataarray(ifile) da_ctl_ts = da da_ref_ts = da.sel(time=slice('0101', '0200')).mean('time') #Ts pert ifile = ifiles[case].replace('precip', 't_surf') print(ifile) da_ts = xr.open_dataarray(ifile) #zz = np.linspace(0, 2, 140*12+1)[1:] # 140 years of monthly CO2 forcing, 1 at 2xCO2 and 2 at 4xCO2 #zz = np.hstack((zz, zz[-1::-1])) #da = xr.DataArray(zz, dims=['time',], coords=[da_ts.time,]).assign_attrs(units='2xCO2', long_name='CO2 forcing') if case == '2xCO2': da_fco2 = da_ts*0+1 da_solar = da_ts*0 elif case == 'm2pctSolar': da_fco2 = da_ts*0 da_solar = da_ts*0-1 elif case == '2xCO2_m2pctSolar': da_fco2 = da_ts*0+1 da_solar = da_ts*0-1 elif case == 'm1p7pctSolar': da_fco2 = da_ts*0 da_solar = da_ts*0-1.7/2 elif case == '2xCO2_m1p7pctSolar': da_fco2 = da_ts*0+1 da_solar = da_ts*0-1.7/2 #da_fco2 = 2 #4xCO2 forcing (2xCO2 forcing is the units) #da.plot(); plt.show(); sys.exit() if __name__ == '__main__': from wyconfig import * #my plot settings func_lp = lambda x: x.rolling(year=9, center=True, min_periods=1).mean('year') func_pct = lambda x: (x/da_ref*100-100).assign_attrs(units='%') fig, ax = plt.subplots() alpha = 0.2 #precip da_ = da_precip.groupby('time.year').mean('time') da_.pipe(func_pct).plot(color='k', alpha=alpha) da_.pipe(func_pct).pipe(func_lp).plot(label=f'precip ({model} {case})', color='k') ax.set_prop_cycle(None) y = da_.pipe(func_pct) #save data for later use #predicted ii = 0; cc = f'C{ii}' if case in ('m2pctSolar','m1p7pctSolar'): da_ = (da_ts.groupby('time.year').mean('time') - da_ref_ts)*b_gmst + da_solar.groupby('time.year').mean('time')*b_solar da_.plot(color=cc, alpha=alpha) da_.pipe(func_lp).plot(label=f'modeled by gmst({b_gmst:.1f}%/K) and solar({b_solar:.1f}%/2%solar)', color=cc) elif case in ('2xCO2_m2pctSolar', '2xCO2_m1p7pctSolar'): da_ = (da_ts.groupby('time.year').mean('time') - da_ref_ts)*b_gmst + da_solar.groupby('time.year').mean('time')*b_solar \ + da_fco2.groupby('time.year').mean('time')*b_fco2 da_.plot(color=cc, alpha=alpha) da_.pipe(func_lp).plot(label=f'modeled by gmst({b_gmst:.1f}%/K) and co2({b_fco2:.1f}%/2xCO2) and solar({b_solar:.1f}%/2%solar)', color=cc) else: da_ = (da_ts.groupby('time.year').mean('time') - da_ref_ts)*b_gmst + da_fco2.groupby('time.year').mean('time')*b_fco2 da_.plot(color=cc, alpha=alpha) da_.pipe(func_lp).plot(label=f'modeled by gmst({b_gmst:.1f}%/K) and co2({b_fco2:.1f}%/2xCO2)', color=cc) y_hat = da_ #save #from gmsta ii += 1; cc = f'C{ii}' da_ = (da_ts.groupby('time.year').mean('time') - da_ref_ts)*b_gmst da_.plot(color=cc, alpha=alpha) da_.pipe(func_lp).plot(label='gmst contribution', color=cc) #from direct rad forcing ii += 1; cc = f'C{ii}' if case in ('m2pctSolar', 'm1p7pctSolar'): da_ = -da_solar.groupby('time.year').mean('time')*b_solar #da_.plot(color=cc, alpha=alpha) da_.plot(label='solar contribution (mirrored)', ls='--', color=cc) elif case in ('2xCO2_m2pctSolar', '2xCO2_m1p7pctSolar'): da_ = da_solar.groupby('time.year').mean('time')*b_solar + da_fco2.groupby('time.year').mean('time')*b_fco2 #da_.plot(color=cc, alpha=alpha) da_.plot(label='both co2 and solar contribution', ls='--', color=cc) else: da_ = -da_fco2.groupby('time.year').mean('time')*b_fco2 #da_.plot(color=cc, alpha=alpha) da_.plot(label='co2 contribution (mirrored)', ls='--', color=cc) #da_ = (da_ts.groupby('time.year').mean('time') - da_ref_ts)*b_gmst_am2p5 + da_fco2.groupby('time.year').mean('time')*b_fco2_am2p5 #da_.pipe(func_lp).plot(label=f'modeled if AM2.5 coeff: gmst({b_gmst_am2p5}%/K) and co2({b_fco2_am2p5}%/2xCO2)') #residual ii += 1; cc = f'C{ii}' da_res = y - y_hat da_res.plot(color=cc, alpha=alpha) da_res.pipe(func_lp).plot(label='residual', color=cc) ax.legend(loc='upper left') ax.set_ylabel('global mean precip change [%]') #ax.set_xticks(range(1000, 1300, 40)) #ax.axvline(1140, color='gray', ls='--') ax.axhline(0, color='gray', ls='--') #if case in ('4xCO2',): # plt.xscale('log') #savefig if 'savefig' in sys.argv or 's' in sys.argv: figname = __file__.replace('.py', f'_{model}_{case}.png') if 'overwritefig' in sys.argv or 'o' in sys.argv: wysavefig(figname, overwritefig=True) else: wysavefig(figname) tt.check(f'**Done**') print() if 'notshowfig' in sys.argv or 'n' in sys.argv: pass else: if 'plt' in globals(): plt.show()