From b026ee7d970e91e9710889e33ab9e0065ad056ec Mon Sep 17 00:00:00 2001 From: Hsiang-He Lee Date: Fri, 22 Sep 2023 15:46:14 -0500 Subject: [PATCH 1/2] save plot data to netcdf files --- ChemDyg_example_script.cfg | 12 +------- chemdyg/templates/chemdyg_STE_flux_diags.bash | 22 +++++++++++++- chemdyg/templates/chemdyg_TOZ_eq_plot.bash | 18 +++++++++-- .../templates/chemdyg_cmip_comparison.bash | 15 ++++++++-- chemdyg/templates/chemdyg_lat_lon_plots.bash | 14 ++++++++- .../templates/chemdyg_noaa_co_comparison.bash | 16 ++++++++-- chemdyg/templates/chemdyg_nox_emis_plots.bash | 20 ++++++++++++- chemdyg/templates/chemdyg_o3_hole_diags.bash | 25 +++++++++++----- chemdyg/templates/chemdyg_pres_lat_plots.bash | 30 +++++++++++++++++-- chemdyg/templates/chemdyg_surf_o3_diags.bash | 18 +++++++++-- .../chemdyg_temperature_eq_plot.bash | 23 ++++++++++++-- 11 files changed, 180 insertions(+), 33 deletions(-) diff --git a/ChemDyg_example_script.cfg b/ChemDyg_example_script.cfg index c55a1ea..0dc04fe 100644 --- a/ChemDyg_example_script.cfg +++ b/ChemDyg_example_script.cfg @@ -59,6 +59,7 @@ active = False active = True short_name = '20230425.amip.e99b02.FourthSmoke_chem_corrected.chrysalis' years = "2013:2013:1", +ncfile_save = 'true' [[index]] grid = 'native' @@ -72,59 +73,48 @@ years = "2013:2013:1", [[climo_diags]] #need preprocessing data from "native_aave" grid = 'native' - years = "2013:2013:1", #should be the same as the preprocessing data [[cmip_comparison]] #need preprocessing data from "atm_monthly_1.0x1.25_aave" grid = '180x360_aave' nodes = 1 - years = "2013:2013:1", #should be the same as the preprocessing data [[noaa_co_comparison]] #need preprocessing data from "atm_monthly_180x360_aave" grid = '180x360_aave' nodes = 1 - years = "2013:2013:1", #should be the same as the preprocessing data [[o3_hole_diags]] grid = 'native' input_files = "eam.h3" frequency = "daily" - years = "2013:2013:1", [[TOZ_eq_native]] grid = 'native' input_files = "eam.h3" frequency = "daily" - years = "2013:2013:1", [[surf_o3_diags]] #need preprocessing data from "atm_hourly_US1.0x1.0_nco" and "atm_hourly_EU1.0x1.0_nco" grid1 = 'MDA8EU1.0x1.0_nco' grid2 = 'MDA8US1.0x1.0_nco' - years = "2013:2013:1", #should be the same as the preprocessing data [[STE_flux_native]] grid = 'native' input_files = "eam.h1" - years = "2013:2013:1", [[summary_table_native]] grid = 'native' input_files = "eam.h1" - years = "2013:2013:1", [[pres_lat_plots]] #need preprocessing data from "180x360_aave" grid = '180x360_aave' nodes = 1 - years = "2013:2013:1", #should be the same as the preprocessing data [[lat_lon_plots]] #need preprocessing data from "180x360_aave" grid = '180x360_aave' nodes = 1 - years = "2013:2013:1", #should be the same as the preprocessing data [[nox_emis_plots]] #need preprocessing data from "180x360_aave" grid = '180x360_aave' nodes = 1 - years = "2013:2013:1", #should be the same as the preprocessing data [[temperature_eq_native]] grid = 'native' diff --git a/chemdyg/templates/chemdyg_STE_flux_diags.bash b/chemdyg/templates/chemdyg_STE_flux_diags.bash index 036374d..20dcc46 100644 --- a/chemdyg/templates/chemdyg_STE_flux_diags.bash +++ b/chemdyg/templates/chemdyg_STE_flux_diags.bash @@ -35,7 +35,15 @@ www="{{ www }}" y1={{ year1 }} y2={{ year2 }} run_type="{{ run_type }}" -tag="{{ tag }}" +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -230,6 +238,15 @@ y_mean = 1.E-9*12*np.array(y.mean()) yN_mean = 1.E-9*12*np.array(yN.mean()) yS_mean = 1.E-9*12*np.array(yS.mean()) +y_xr = y.assign_attrs(units="kg/month", description="O3 flux") +yN_xr = yN.assign_attrs(units="kg/month", description="NH O3 flux") +yS_xr = yS.assign_attrs(units="kg/month", description="SH O3 flux") +ds1 = y_xr.to_dataset(name='O3_STE') +ds2 = yN_xr.to_dataset(name='O3_STE_NH') +ds3 = yS_xr.to_dataset(name='O3_STE_SH') +ds = xr.merge([ds1, ds2, ds3]) +ds.to_netcdf(pathout+'E3SM_O3_STE_${y1}-${y2}.nc') + # time series plot fig = plt.figure(figsize=(10,5)) plt.plot(time_range_month[1::],y*1.E-9*12) @@ -295,6 +312,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_TOZ_eq_plot.bash b/chemdyg/templates/chemdyg_TOZ_eq_plot.bash index cff8f6c..4ebcdf0 100644 --- a/chemdyg/templates/chemdyg_TOZ_eq_plot.bash +++ b/chemdyg/templates/chemdyg_TOZ_eq_plot.bash @@ -37,8 +37,15 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -157,6 +164,10 @@ for i in range(startindex,endindex): TOZ_index = result[0].min() O3_64S[ii] = (TOZ_sort[0:TOZ_index-1]*AREA_sort[0:TOZ_index-1]).mean()/(AREA_sort[0:TOZ_index-1].mean()) +O3_64S_xr = xr.DataArray(O3_64S, name = 'TOZ_64S',coords=[time_range[startindex:endindex]], dims=["time"]) +O3_64S_xr = O3_64S_xr.assign_attrs(units="DU", description='Total column ozone conc. with equivalent latitude (64S)') +O3_64S_xr.to_netcdf(pathout+'E3SM_TOZ_64S_${y1}-${y2}.nc') + fig = plt.figure(figsize=(10,5)) plt.plot(time_range[startindex:endindex],O3_64S) plt.title('Total column ozone conc. with equivalent latitude (64S)') @@ -210,6 +221,9 @@ mkdir -p ${f} if [ -d "${f}" ]; then mv ./*.png ${f} fi +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_cmip_comparison.bash b/chemdyg/templates/chemdyg_cmip_comparison.bash index 5faa104..01be92d 100644 --- a/chemdyg/templates/chemdyg_cmip_comparison.bash +++ b/chemdyg/templates/chemdyg_cmip_comparison.bash @@ -37,10 +37,17 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" +ncfile_save="{{ ncfile_save }}" # diagnostics_base_path is set by zppy using the mache package cmipDir="{{ diagnostics_base_path }}/observations/Atm/ChemDyg_inputs" -results_dir=${tag}_${Y1}-${Y2} +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -167,6 +174,7 @@ E3SM_long = np.zeros(t_end) E3SM_long[0:year_start] = 'NAN' E3SM_long[year_start:t_end] = E3SM_TCO[:].values E3SM_xr = xr.DataArray(E3SM_long, coords=[CESM_in['time'][0:t_end]], dims=["time"], name='TCO') +E3SM_xr.to_netcdf(pathout+'E3SM_cmip_'+startyear+'-'+endyear+'.nc') time_range_year = pd.date_range('1850-01-01', periods=timeperiod_year, freq='Y') CESM_ANN = CESM_TCO.groupby('time.year').mean('time') @@ -233,6 +241,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_lat_lon_plots.bash b/chemdyg/templates/chemdyg_lat_lon_plots.bash index b6f6cbb..fa7cc88 100644 --- a/chemdyg/templates/chemdyg_lat_lon_plots.bash +++ b/chemdyg/templates/chemdyg_lat_lon_plots.bash @@ -37,7 +37,15 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -121,6 +129,7 @@ filename = short_name+'_ANN_'+startyear+'01_'+endyear+'12_climo.nc' file_in = xr.open_dataset(path+filename) TMQ = file_in['TMQ'][0,:,:] +TMQ.to_netcdf(pathout+'TMQ_'+startyear+'-'+endyear+'.nc') fig = plt.figure(figsize=(20,10)) ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) @@ -168,6 +177,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_noaa_co_comparison.bash b/chemdyg/templates/chemdyg_noaa_co_comparison.bash index c654fd2..f6fc3ae 100644 --- a/chemdyg/templates/chemdyg_noaa_co_comparison.bash +++ b/chemdyg/templates/chemdyg_noaa_co_comparison.bash @@ -37,10 +37,17 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" # diagnostics_base_path is set by zppy using the mache package noaaDir="{{ diagnostics_base_path }}/observations/Atm/ChemDyg_inputs" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -189,6 +196,8 @@ for n in range(len(Sta)): lin_noaa = nmonth_noaa*slope_noaa+intercept lin_noaa_xa = xr.DataArray(lin_noaa, coords=[time_range_noaa], dims=["time"]) diff_noaa = CO_noaa - lin_noaa_xa + CO_sel_1D_xr = xr.DataArray(CO_sel_1D, name='CO', coords=[time_range_noaa], dims=["time"]) + CO_sel_1D_xr.to_netcdf(pathout+'E3SM_CO_'+startyear+'-'+endyear+'.nc') nmonth = np.arange(0,len(CO_sel_1D),1) mask = ~np.isnan(CO_sel_1D) @@ -261,6 +270,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_nox_emis_plots.bash b/chemdyg/templates/chemdyg_nox_emis_plots.bash index a7f816e..9d31e16 100644 --- a/chemdyg/templates/chemdyg_nox_emis_plots.bash +++ b/chemdyg/templates/chemdyg_nox_emis_plots.bash @@ -37,7 +37,15 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -136,6 +144,13 @@ for i in range(len(lev)): NOx_acf_3d[i,:,:] = NOx_acf[i,:,:] * AREA # Tg N/m2/year ->Tg N/year NOx_lgt_3d[i,:,:] = NOx_lgt[i,:,:] * AREA +NOx_acf_xr = NOx_acf_3d.assign_attrs(units="Tg N/year", description="NOx aircraft emission") +NOx_lgt_xr = NOx_lgt_3d.assign_attrs(units="Tg N/year", description="NOx lightning emission") +ds1 = NOx_acf_xr.to_dataset() +ds2 = NOx_lgt_xr.to_dataset() +ds = xr.merge([ds1, ds2]) +ds.to_netcdf(pathout+'E3SM_NOx_emission_'+startyear+'-'+endyear+'.nc') + NOx_acf_1d = NOx_acf_3d.sum(axis=1).sum(axis=1) NOx_acf_2d = NOx_acf_3d.sum(axis=0) NOx_lgt_1d = NOx_lgt_3d.sum(axis=1).sum(axis=1) @@ -216,6 +231,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_o3_hole_diags.bash b/chemdyg/templates/chemdyg_o3_hole_diags.bash index b75994e..91d480d 100644 --- a/chemdyg/templates/chemdyg_o3_hole_diags.bash +++ b/chemdyg/templates/chemdyg_o3_hole_diags.bash @@ -38,10 +38,17 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" # diagnostics_base_path is set by zppy using the mache package obsDir="{{ diagnostics_base_path }}/observations/Atm/ChemDyg_inputs" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -49,15 +56,11 @@ cd ${workdir} # Create local links to input climo files tsDir={{ output }}/post/atm/{{ grid }}/ts/monthly/{{ '%dyr' % (ypf) }} -#tsDir={{ output }}/post/atm/{{ grid }} mkdir -p ts mkdir -p figs mkdir -p data ln -s ${obsDir}/O3_hole/*_obs.nc ./ts -#cd ts -#ln -s ${tsDir}/TCO*.nc ./ts -#ln -s ${cmipDir}/*TCO.nc ./ts -#cd .. + # Create symbolic links to input files input={{ input }}/{{ input_subdir }} eamfile={{ input_files }} @@ -164,6 +167,11 @@ for i in range(startindex,endindex): else: O3_area[i] = 0. +ds1 = O3_area.to_dataset(name='O3_area') +ds2 = TOZ_min.to_dataset(name='TOZ_min') +ds1.to_netcdf(pathout+'E3SM_O3_area_${y1}-${y2}.nc') +ds2.to_netcdf(pathout+'E3SM_TOZ_min_${y1}-${y2}.nc') + O3_area_time = O3_area[startindex:endindex].sel(time=O3_area[startindex:endindex].time.dt.month.isin([7, 8, 9, 10, 11, 12])) TOZ_min_time = TOZ_min[startindex:endindex].sel(time=TOZ_min[startindex:endindex].time.dt.month.isin([7, 8, 9, 10, 11, 12])) @@ -225,6 +233,9 @@ mkdir -p ${f} if [ -d "${f}" ]; then mv ./*.png ${f} fi +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_pres_lat_plots.bash b/chemdyg/templates/chemdyg_pres_lat_plots.bash index 7bdd81a..a62598d 100644 --- a/chemdyg/templates/chemdyg_pres_lat_plots.bash +++ b/chemdyg/templates/chemdyg_pres_lat_plots.bash @@ -37,10 +37,17 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" # diagnostics_base_path is set by zppy using the mache package referDir="{{ diagnostics_base_path }}/observations/Atm/ChemDyg_inputs" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -141,6 +148,10 @@ o3_2d = o3.mean(axis=2) o3_refer_2d = o3_refer.mean(axis=2) o3_new = o3_2d.copy() o3_refer_new = o3_refer_2d.copy() +ds1 = o3_new.to_dataset() +ds2 = o3_refer_new.to_dataset(name="O3_v2") +ds = xr.merge([ds1, ds2],compat='override') +ds.to_netcdf(pathout+'E3SM_O3_PLatcomparison_'+startyear+'-'+endyear+'.nc') Q = file_in['Q'][0,:,:,:]*28.96/18 Q_refer = refer_in['Q'][0,:,:,:]*28.96/18 @@ -148,6 +159,10 @@ Q_2d = Q.mean(axis=2) Q_refer_2d = Q_refer.mean(axis=2) Q_new = Q_2d.copy() Q_refer_new = Q_refer_2d.copy() +ds1 = Q_new.to_dataset() +ds2 = Q_refer_new.to_dataset(name="Q_v2") +ds = xr.merge([ds1, ds2],compat='override') +ds.to_netcdf(pathout+'E3SM_Q_PLatcomparison_'+startyear+'-'+endyear+'.nc') T = file_in['T'][0,:,:,:] T_refer = refer_in['T'][0,:,:,:] @@ -155,6 +170,10 @@ T_2d = T.mean(axis=2) T_refer_2d = T_refer.mean(axis=2) T_new = T_2d.copy() T_refer_new = T_refer_2d.copy() +ds1 = T_new.to_dataset() +ds2 = T_refer_new.to_dataset(name="T_v2") +ds = xr.merge([ds1, ds2],compat='override') +ds.to_netcdf(pathout+'E3SM_Temp_PLatcomparison_'+startyear+'-'+endyear+'.nc') theda = T*(100000/lev)**0.286 theda_refer = T_refer*(100000/lev)**0.286 @@ -162,6 +181,10 @@ theda_2d = theda.mean(axis=2) theda_refer_2d = theda_refer.mean(axis=2) theda_new = theda_2d.copy() theda_refer_new = theda_refer_2d.copy() +ds1 = theda_new.to_dataset(name="theda") +ds2 = theda_refer_new.to_dataset(name="theda_v2") +ds = xr.merge([ds1, ds2],compat='override') +ds.to_netcdf(pathout+'E3SM_theda_PLatcomparison_'+startyear+'-'+endyear+'.nc') for i in range(180): f = interpolate.interp1d(lev_2d[:,i], o3_2d[:,i]) @@ -470,6 +493,9 @@ fi # Copy files mv *.png ${f} +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_surf_o3_diags.bash b/chemdyg/templates/chemdyg_surf_o3_diags.bash index 6aa2ad7..2adb8f2 100644 --- a/chemdyg/templates/chemdyg_surf_o3_diags.bash +++ b/chemdyg/templates/chemdyg_surf_o3_diags.bash @@ -37,10 +37,17 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" # diagnostics_base_path is set by zppy using the mache package obsDir="{{ diagnostics_base_path }}/observations/Atm/ChemDyg_inputs" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -193,6 +200,10 @@ O3local_US[O3local_US == 0.] = 'nan' # ## calculate diurnal cycle O3EU_xr = xr.DataArray(O3local_EU, coords=[time_EU,lat_EU,lon_EU], dims=["time","lat","lon"]) O3US_xr = xr.DataArray(O3local_US, coords=[time_US,lat_US,lon_US], dims=["time","lat","lon"]) +ds1 = O3EU_xr.to_dataset(name='sfcO3_EU') +ds2 = O3US_xr.to_dataset(name='sfcO3_US') +ds = xr.merge([ds1, ds2]) +ds.to_netcdf(pathout+'E3SM_sfcO3_${y1}-${y2}.nc') O3EU_JJA = O3EU_xr.sel(time=O3EU_xr.time.dt.month.isin([6, 7, 8])) O3EU_DJF = O3EU_xr.sel(time=O3EU_xr.time.dt.month.isin([12, 1, 2])) @@ -778,6 +789,9 @@ mkdir -p ${f} if [ -d "${f}" ]; then mv ./*.png ${f} fi +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} diff --git a/chemdyg/templates/chemdyg_temperature_eq_plot.bash b/chemdyg/templates/chemdyg_temperature_eq_plot.bash index c49fa94..f7305b8 100644 --- a/chemdyg/templates/chemdyg_temperature_eq_plot.bash +++ b/chemdyg/templates/chemdyg_temperature_eq_plot.bash @@ -37,8 +37,15 @@ y2={{ year2 }} Y1="{{ '%04d' % (year1) }}" Y2="{{ '%04d' % (year2) }}" run_type="{{ run_type }}" -tag="{{ tag }}" -results_dir=${tag}_${Y1}-${Y2} +ncfile_save="{{ ncfile_save }}" +if [[ "${ncfile_save}" == "true" ]]; then + results_dir={{ output }}/post/atm/ncfiles + if [[ -d ${results_dir} ]]; then + echo "directory exists." + else + mkdir -p ${results_dir} + fi +fi # Create temporary workdir workdir=`mktemp -d tmp.${id}.XXXX` @@ -198,6 +205,15 @@ for i in range(0,20,2): temp_lat_l[ii] = temp_l.mean() latlist = np.arange(60,80,2) +temp_lat_h_xr = xr.DataArray(temp_lat_h, name = 'temp_lat_h',coords=[latlist], dims=["lat"]) +temp_lat_m_xr = xr.DataArray(temp_lat_m, name = 'temp_lat_m',coords=[latlist], dims=["lat"]) +temp_lat_l_xr = xr.DataArray(temp_lat_l, name = 'temp_lat_l',coords=[latlist], dims=["lat"]) +ds1 = temp_lat_h_xr.to_dataset() +ds2 = temp_lat_m_xr.to_dataset() +ds3 = temp_lat_l_xr.to_dataset() +ds = xr.merge([ds1, ds2, ds3]) +ds.to_netcdf(pathout+'E3SM_O3_equi_lat_${y1}-${y2}.nc') + fig = plt.figure(figsize=(10,5)) plt.plot(latlist,temp_lat_h) plt.plot(latlist,temp_lat_m) @@ -230,6 +246,9 @@ mkdir -p ${f} if [ -d "${f}" ]; then mv ./*.png ${f} fi +if [[ "${ncfile_save}" == "true" ]]; then + mv *.nc ${results_dir} +fi # Change file permissions chmod -R go+rX,go-w ${f} From 22f5e811e3ea278bbcbcc07b998a5ba281d12c84 Mon Sep 17 00:00:00 2001 From: Hsiang-He Lee Date: Sat, 23 Sep 2023 18:37:06 -0500 Subject: [PATCH 2/2] saving plotting data (1D or 2D array) in the netcdf format --- chemdyg/templates/chemdyg_STE_flux_diags.bash | 44 ++++++--- chemdyg/templates/chemdyg_TOZ_eq_plot.bash | 38 +++++--- .../templates/chemdyg_cmip_comparison.bash | 30 +++++- chemdyg/templates/chemdyg_lat_lon_plots.bash | 5 +- .../templates/chemdyg_noaa_co_comparison.bash | 48 +++++++--- chemdyg/templates/chemdyg_nox_emis_plots.bash | 20 ++-- chemdyg/templates/chemdyg_o3_hole_diags.bash | 31 ++++-- chemdyg/templates/chemdyg_pres_lat_plots.bash | 65 +++++++++---- chemdyg/templates/chemdyg_surf_o3_diags.bash | 96 +++++++++++++------ .../chemdyg_temperature_eq_plot.bash | 13 ++- 10 files changed, 286 insertions(+), 104 deletions(-) diff --git a/chemdyg/templates/chemdyg_STE_flux_diags.bash b/chemdyg/templates/chemdyg_STE_flux_diags.bash index 20dcc46..7ccb6c2 100644 --- a/chemdyg/templates/chemdyg_STE_flux_diags.bash +++ b/chemdyg/templates/chemdyg_STE_flux_diags.bash @@ -234,20 +234,42 @@ else: yS_ann[0] = 'nan' yS_ann[1:12] = y[0:11]*1.E-9*12 +# ----- writing ncfile ----- +y_xr = xr.DataArray(y, coords=[time_range_month[1::]], dims=["time"]) +yN_xr = xr.DataArray(yN, coords=[time_range_month[1::]], dims=["time"]) +yS_xr = xr.DataArray(yS, coords=[time_range_month[1::]], dims=["time"]) +y_xr = y_xr.assign_attrs(units="kg/month", description="O3 flux") +yN_xr = yN_xr.assign_attrs(units="kg/month", description="NH O3 flux") +yS_xr = yS_xr.assign_attrs(units="kg/month", description="SH O3 flux") +ds1 = y_xr.to_dataset(name='y') +ds2 = yN_xr.to_dataset(name='yN') +ds3 = yS_xr.to_dataset(name='yS') +y_ann_xr = xr.DataArray(y_ann, coords=[np.arange(1,13)], dims=["month"]) +y_std_xr = xr.DataArray(y_std, coords=[np.arange(1,13)], dims=["month"]) +y_ann_xr = y_ann_xr.assign_attrs(units="Tg/year", description='O3 flux mean') +y_std_xr = y_std_xr.assign_attrs(units="Tg/year", description='O3 flux standard deviation') +yN_ann_xr = xr.DataArray(yN_ann, coords=[np.arange(1,13)], dims=["month"]) +yN_std_xr = xr.DataArray(yN_std, coords=[np.arange(1,13)], dims=["month"]) +yN_ann_xr = yN_ann_xr.assign_attrs(units="Tg/year", description='NH O3 flux mean') +yN_std_xr = yN_std_xr.assign_attrs(units="Tg/year", description='NH O3 flux standard deviation') +yS_ann_xr = xr.DataArray(yS_ann, coords=[np.arange(1,13)], dims=["month"]) +yS_std_xr = xr.DataArray(yS_std, coords=[np.arange(1,13)], dims=["month"]) +yS_ann_xr = yS_ann_xr.assign_attrs(units="Tg/year", description='SH O3 flux mean') +yS_std_xr = yS_std_xr.assign_attrs(units="Tg/year", description='SH O3 flux standard deviation') +ds4 = y_ann_xr.to_dataset(name='y_ann') +ds5 = yN_ann_xr.to_dataset(name='yN_ann') +ds6 = yS_ann_xr.to_dataset(name='yS_ann') +ds7 = y_ann_xr.to_dataset(name='y_std') +ds8 = yN_ann_xr.to_dataset(name='yN_std') +ds9 = yS_ann_xr.to_dataset(name='yS_std') +ds = xr.merge([ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8, ds9]) +ds.to_netcdf(pathout+'E3SM_O3_STE_${y1}-${y2}.nc') + +# ----- time series plot ----- y_mean = 1.E-9*12*np.array(y.mean()) yN_mean = 1.E-9*12*np.array(yN.mean()) yS_mean = 1.E-9*12*np.array(yS.mean()) -y_xr = y.assign_attrs(units="kg/month", description="O3 flux") -yN_xr = yN.assign_attrs(units="kg/month", description="NH O3 flux") -yS_xr = yS.assign_attrs(units="kg/month", description="SH O3 flux") -ds1 = y_xr.to_dataset(name='O3_STE') -ds2 = yN_xr.to_dataset(name='O3_STE_NH') -ds3 = yS_xr.to_dataset(name='O3_STE_SH') -ds = xr.merge([ds1, ds2, ds3]) -ds.to_netcdf(pathout+'E3SM_O3_STE_${y1}-${y2}.nc') - -# time series plot fig = plt.figure(figsize=(10,5)) plt.plot(time_range_month[1::],y*1.E-9*12) plt.plot(time_range_month[1::],yN*1.E-9*12) @@ -262,7 +284,7 @@ line3 = 'SH mean:'+str(np.round(yS_mean,2)) plt.legend( [line1,line2,line3]) pylab.savefig(pathout+'O3_STE_flux.png', dpi=300) -# annual plot +# ----- annual plot ----- fig = plt.figure(figsize=(10,5)) month_txt = np.arange(1,13) plt.plot(month_txt,y_ann, linewidth = 1) diff --git a/chemdyg/templates/chemdyg_TOZ_eq_plot.bash b/chemdyg/templates/chemdyg_TOZ_eq_plot.bash index 4ebcdf0..dbc5099 100644 --- a/chemdyg/templates/chemdyg_TOZ_eq_plot.bash +++ b/chemdyg/templates/chemdyg_TOZ_eq_plot.bash @@ -164,10 +164,33 @@ for i in range(startindex,endindex): TOZ_index = result[0].min() O3_64S[ii] = (TOZ_sort[0:TOZ_index-1]*AREA_sort[0:TOZ_index-1]).mean()/(AREA_sort[0:TOZ_index-1].mean()) -O3_64S_xr = xr.DataArray(O3_64S, name = 'TOZ_64S',coords=[time_range[startindex:endindex]], dims=["time"]) -O3_64S_xr = O3_64S_xr.assign_attrs(units="DU", description='Total column ozone conc. with equivalent latitude (64S)') -O3_64S_xr.to_netcdf(pathout+'E3SM_TOZ_64S_${y1}-${y2}.nc') +O3_array = O3_64S.reshape((years,365)) +O3_mean = O3_array.mean(axis=0) +O3_std = O3_array.std(axis=0) +TOZ_array = TOZ_min[startindex:endindex].values.reshape((years,365)) +TOZ_mean = TOZ_array.mean(axis=0) +TOZ_std = TOZ_array.std(axis=0) +# ----- writing ncfile ----- +O3_64S_xr = xr.DataArray(O3_64S, name = 'TOZ_64S',coords=[time_range[startindex:endindex]], dims=["time"], + attrs=dict(units="DU", description='Total column ozone conc. with equivalent latitude (64S)')) +TOZ_mean_xr = xr.DataArray(TOZ_mean, name = 'TOZ_mean', coords=[time_range[0:365]], dims=["day"], + attrs=dict(units="DU", description='Climatological mean of daily SH minimum E3SM total column ozone')) +TOZ_std_xr = xr.DataArray(TOZ_std, name = 'TOZ_std', coords=[time_range[0:365]], dims=["day"], + attrs=dict(units="DU", description='Standard deviation of daily SH minimum E3SM total column ozone')) +O3_mean_xr = xr.DataArray(O3_mean, name = 'O3_mean', coords=[time_range[0:365]], dims=["day"], + attrs=dict(units="DU", description='Mean of total column ozone conc. with equivalent latitude (64S)')) +O3_std_xr = xr.DataArray(O3_std, name = 'O3_std', coords=[time_range[0:365]], dims=["day"], + attrs=dict(units="DU", description='Standard deviation of total column ozone conc. with equivalent latitude (64S)')) +ds1 = TOZ_mean_xr.to_dataset() +ds2 = TOZ_std_xr.to_dataset() +ds3 = O3_mean_xr.to_dataset() +ds4 = O3_std_xr.to_dataset() +ds5 = O3_64S_xr.to_dataset() +ds = xr.merge([ds1, ds2, ds3, ds4,ds5]) +ds.to_netcdf(pathout+'E3SM_TOZ_64S_${y1}-${y2}.nc') + +# ----- plotting ----- fig = plt.figure(figsize=(10,5)) plt.plot(time_range[startindex:endindex],O3_64S) plt.title('Total column ozone conc. with equivalent latitude (64S)') @@ -176,16 +199,9 @@ plt.ylabel('O3 conc. (DU)',fontsize='large') pylab.savefig(pathout+'TOZ_PDF_timeseries.png', dpi=300) #-------- climo plot --------------- -O3_array = O3_64S.reshape((years,365)) -O3_mean = O3_array.mean(axis=0) -O3_std = O3_array.std(axis=0) -TOZ_array = TOZ_min[startindex:endindex].values.reshape((years,365)) -TOZ_mean = TOZ_array.mean(axis=0) -TOZ_std = TOZ_array.std(axis=0) - npdate = np.array(time_range[0:365]) fig, ax = plt.subplots(figsize=(10, 5)) -plt.plot(npdate,TOZ_mean, label ='E3SMv2') +plt.plot(npdate,TOZ_mean, label ='E3SM') plt.fill_between(npdate,TOZ_mean+TOZ_std,TOZ_mean-TOZ_std, alpha=.5, linewidth=0) plt.plot(npdate,O3_mean, label ='TOZ(64S)') plt.fill_between(npdate,O3_mean+O3_std,O3_mean-O3_std, alpha=.5, linewidth=0) diff --git a/chemdyg/templates/chemdyg_cmip_comparison.bash b/chemdyg/templates/chemdyg_cmip_comparison.bash index 01be92d..3fa7c9d 100644 --- a/chemdyg/templates/chemdyg_cmip_comparison.bash +++ b/chemdyg/templates/chemdyg_cmip_comparison.bash @@ -174,7 +174,6 @@ E3SM_long = np.zeros(t_end) E3SM_long[0:year_start] = 'NAN' E3SM_long[year_start:t_end] = E3SM_TCO[:].values E3SM_xr = xr.DataArray(E3SM_long, coords=[CESM_in['time'][0:t_end]], dims=["time"], name='TCO') -E3SM_xr.to_netcdf(pathout+'E3SM_cmip_'+startyear+'-'+endyear+'.nc') time_range_year = pd.date_range('1850-01-01', periods=timeperiod_year, freq='Y') CESM_ANN = CESM_TCO.groupby('time.year').mean('time') @@ -190,6 +189,35 @@ UKESM_std = UKESM_TCO.groupby('time.year').std('time') E3SM_ANN = E3SM_xr.groupby('time.year').mean('time') E3SM_std = E3SM_xr.groupby('time.year').std('time') +# ----- writing ncfile ----- +CESM_ANN = CESM_ANN.assign_attrs(units="Tg", description='CESM Tropospheric-ozone burden (TCO) mean') +CESM_std = CESM_std.assign_attrs(units="Tg", description='CESM Tropospheric-ozone burden (TCO) std') +CESM_ANN_xr = CESM_ANN.to_dataset(name='CESM_ANN') +CESM_std_xr = CESM_std.to_dataset(name='CESM_std') +GFDL_ANN = GFDL_ANN.assign_attrs(units="Tg", description='GFDL Tropospheric-ozone burden (TCO) mean') +GFDL_std = GFDL_std.assign_attrs(units="Tg", description='GFDL Tropospheric-ozone burden (TCO) std') +GFDL_ANN_xr = GFDL_ANN.to_dataset(name='GFDL_ANN') +GFDL_std_xr = GFDL_std.to_dataset(name='GFDL_std') +GISS_ANN = GISS_ANN.assign_attrs(units="Tg", description='GISS Tropospheric-ozone burden (TCO) mean') +GISS_std = GISS_std.assign_attrs(units="Tg", description='GISS Tropospheric-ozone burden (TCO) std') +GISS_ANN_xr = GISS_ANN.to_dataset(name='GISS_ANN') +GISS_std_xr = GISS_std.to_dataset(name='GISS_std') +MRI_ANN = MRI_ANN.assign_attrs(units="Tg", description='MRI Tropospheric-ozone burden (TCO) mean') +MRI_std = MRI_std.assign_attrs(units="Tg", description='MRI Tropospheric-ozone burden (TCO) std') +MRI_ANN_xr = MRI_ANN.to_dataset(name='MRI_ANN') +MRI_std_xr = MRI_std.to_dataset(name='MRI_std') +UKESM_ANN = UKESM_ANN.assign_attrs(units="Tg", description='UKESM Tropospheric-ozone burden (TCO) mean') +UKESM_std = UKESM_std.assign_attrs(units="Tg", description='UKESM Tropospheric-ozone burden (TCO) std') +UKESM_ANN_xr = UKESM_ANN.to_dataset(name='UKESM_ANN') +UKESM_std_xr = UKESM_std.to_dataset(name='UKESM_std') +E3SM_ANN = E3SM_ANN.assign_attrs(units="Tg", description='E3SM Tropospheric-ozone burden (TCO) mean') +E3SM_std = E3SM_std.assign_attrs(units="Tg", description='E3SM Tropospheric-ozone burden (TCO) std') +E3SM_ANN_xr = E3SM_ANN.to_dataset(name='E3SM_ANN') +E3SM_std_xr = E3SM_std.to_dataset(name='E3SM_std') +ds = xr.merge([CESM_ANN_xr, CESM_std_xr,GFDL_ANN_xr, GFDL_std_xr,GISS_ANN_xr, GISS_std_xr, + MRI_ANN_xr, MRI_std_xr, UKESM_ANN_xr, UKESM_std_xr, E3SM_ANN_xr, E3SM_std_xr]) +ds.to_netcdf(pathout+'E3SM_cmip_'+startyear+'-'+endyear+'.nc') + #----- plotting ----- fig = plt.figure(figsize=(10,5)) plt.plot(time_range_year,CESM_ANN, label='CESM', linewidth = 1) diff --git a/chemdyg/templates/chemdyg_lat_lon_plots.bash b/chemdyg/templates/chemdyg_lat_lon_plots.bash index fa7cc88..47c9f03 100644 --- a/chemdyg/templates/chemdyg_lat_lon_plots.bash +++ b/chemdyg/templates/chemdyg_lat_lon_plots.bash @@ -129,8 +129,11 @@ filename = short_name+'_ANN_'+startyear+'01_'+endyear+'12_climo.nc' file_in = xr.open_dataset(path+filename) TMQ = file_in['TMQ'][0,:,:] -TMQ.to_netcdf(pathout+'TMQ_'+startyear+'-'+endyear+'.nc') +# ----- writing ncfile ----- +TMQ.to_netcdf(pathout+'E3SM_TMQ_'+startyear+'-'+endyear+'.nc') + +# ----- plotting ----- fig = plt.figure(figsize=(20,10)) ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) diff --git a/chemdyg/templates/chemdyg_noaa_co_comparison.bash b/chemdyg/templates/chemdyg_noaa_co_comparison.bash index f6fc3ae..40c3f63 100644 --- a/chemdyg/templates/chemdyg_noaa_co_comparison.bash +++ b/chemdyg/templates/chemdyg_noaa_co_comparison.bash @@ -194,10 +194,8 @@ for n in range(len(Sta)): mask_noaa = ~np.isnan(CO_noaa) slope_noaa, intercept, r_value, p_value, std_err = linregress(nmonth_noaa[mask_noaa], CO_noaa[mask_noaa]) lin_noaa = nmonth_noaa*slope_noaa+intercept - lin_noaa_xa = xr.DataArray(lin_noaa, coords=[time_range_noaa], dims=["time"]) - diff_noaa = CO_noaa - lin_noaa_xa - CO_sel_1D_xr = xr.DataArray(CO_sel_1D, name='CO', coords=[time_range_noaa], dims=["time"]) - CO_sel_1D_xr.to_netcdf(pathout+'E3SM_CO_'+startyear+'-'+endyear+'.nc') + lin_noaa_xr = xr.DataArray(lin_noaa, coords=[time_range_noaa], dims=["time"]) + diff_noaa = CO_noaa - lin_noaa_xr nmonth = np.arange(0,len(CO_sel_1D),1) mask = ~np.isnan(CO_sel_1D) @@ -205,15 +203,32 @@ for n in range(len(Sta)): if (len(nmonth[mask]) != 0): slope_e3sm, intercept, r_value, p_value, std_err = linregress(nmonth[mask], CO_sel_1D[mask]) lin_e3sm = nmonth*slope_e3sm+intercept - lin_e3sm_xa = xr.DataArray(lin_e3sm, coords=[time_range_noaa], dims=["time"]) - diff = CO_sel_1D - lin_e3sm_xa - - # plotting + lin_e3sm_xr = xr.DataArray(lin_e3sm, coords=[time_range_noaa], dims=["time"]) + diff = CO_sel_1D - lin_e3sm_xr + + # ----- writing ncfile ----- + CO_sel_1D_xr = xr.DataArray(CO_sel_1D, name='CO_sel_1D', coords=[CO_noaa['time']], dims=["time"]) + CO_sel_1D_xr = CO_sel_1D_xr.assign_attrs(units="ppb", description='E3SM surface CO') + diff_xr = xr.DataArray(diff, name='diff', coords=[CO_noaa['time']], dims=["time"]) + diff_xr = diff_xr.assign_attrs(units="ppb", description='E3SM surface CO anomaly') + lin_e3sm_xr = lin_e3sm_xr.assign_attrs(units="ppb", description='Trend E3SM surface CO') + lin_noaa_xr = lin_noaa_xr.assign_attrs(units="ppb", description='Trend NOAA surface CO') + diff_noaa = diff_noaa.assign_attrs(units="ppb", description='NOAA surface CO anomaly') + ds1 = CO_sel_1D_xr.to_dataset() + ds2 = CO_noaa.to_dataset(name='CO_noaa') + ds3 = lin_e3sm_xr.to_dataset(name='lin_e3sm_xr') + ds4 = lin_noaa_xr.to_dataset(name='lin_noaa_xr') + ds5 = diff_xr.to_dataset() + ds6 = diff_noaa.to_dataset(name='diff_noaa') + ds = xr.merge([ds1,ds2,ds3,ds4,ds5,ds6]) + ds.to_netcdf(pathout+'E3SM_CO_'+Sta[n]+'_'+startyear+'-'+endyear+'.nc') + + # ----- plotting ----- fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5)) ax1.plot(CO_noaa['time'],CO_sel_1D,'k') - ax1.plot(CO_noaa['time'][mask],lin_e3sm_xa[mask],'k--') + ax1.plot(CO_noaa['time'][mask],lin_e3sm_xr[mask],'k--') ax1.plot(CO_noaa['time'],CO_noaa,'r') - ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xa[mask_noaa],'r--') + ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xr[mask_noaa],'r--') ax1.set_title('Surface CO at '+Sta[n]+' (Lat '+str(lat_noaa[0].values)+', Lon '+str(lon_noaa[0].values)+')') line1 = 'E3SM mean:'+str(np.round(CO_sel_1D[mask].mean(),2)) line2 = 'E3SM trend:'+str(np.round(slope_e3sm*12,2))+' ppb/yr' @@ -228,10 +243,19 @@ for n in range(len(Sta)): ax2.set_ylabel('Anomalies') pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600) else: - # plotting + # ----- writing ncfile ----- + lin_noaa_xr = lin_noaa_xr.assign_attrs(units="ppb", description='Trend NOAA surface CO') + diff_noaa = diff_noaa.assign_attrs(units="ppb", description='NOAA surface CO anomaly') + ds1 = CO_noaa.to_dataset(name='CO_noaa') + ds2 = lin_noaa_xr.to_dataset(name='lin_noaa_xr') + ds3 = diff_noaa.to_dataset(name='diff_noaa') + ds = xr.merge([ds1,ds2,ds3]) + ds.to_netcdf(pathout+'E3SM_CO_'+Sta[n]+'_'+startyear+'-'+endyear+'.nc') + + # ----- plotting ------ fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5)) ax1.plot(CO_noaa['time'],CO_noaa,'r') - ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xa[mask_noaa],'r--') + ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xr[mask_noaa],'r--') ax1.set_title('Surface CO at '+Sta[n]+' (Lat '+str(lat_noaa[0].values)+', Lon '+str(lon_noaa[0].values)+')') line3 = 'NOAA mean:'+str(np.round(CO_noaa[mask_noaa].mean().values,2)) line4 = 'NOAA trend:'+str(np.round(slope_noaa*12,2))+' ppb/yr' diff --git a/chemdyg/templates/chemdyg_nox_emis_plots.bash b/chemdyg/templates/chemdyg_nox_emis_plots.bash index 9d31e16..421a670 100644 --- a/chemdyg/templates/chemdyg_nox_emis_plots.bash +++ b/chemdyg/templates/chemdyg_nox_emis_plots.bash @@ -144,20 +144,24 @@ for i in range(len(lev)): NOx_acf_3d[i,:,:] = NOx_acf[i,:,:] * AREA # Tg N/m2/year ->Tg N/year NOx_lgt_3d[i,:,:] = NOx_lgt[i,:,:] * AREA -NOx_acf_xr = NOx_acf_3d.assign_attrs(units="Tg N/year", description="NOx aircraft emission") -NOx_lgt_xr = NOx_lgt_3d.assign_attrs(units="Tg N/year", description="NOx lightning emission") -ds1 = NOx_acf_xr.to_dataset() -ds2 = NOx_lgt_xr.to_dataset() -ds = xr.merge([ds1, ds2]) -ds.to_netcdf(pathout+'E3SM_NOx_emission_'+startyear+'-'+endyear+'.nc') - NOx_acf_1d = NOx_acf_3d.sum(axis=1).sum(axis=1) NOx_acf_2d = NOx_acf_3d.sum(axis=0) NOx_lgt_1d = NOx_lgt_3d.sum(axis=1).sum(axis=1) NOx_lgt_2d = NOx_lgt_3d.sum(axis=0) -# plotting +# ----- writing ncfile ----- +NOx_acf_1d = NOx_acf_1d.assign_attrs(units="Tg N/year", description="NOx aircraft emission") +NOx_acf_2d = NOx_acf_2d.assign_attrs(units="Tg N/year", description="NOx aircraft emission") +NOx_lgt_1d = NOx_lgt_1d.assign_attrs(units="Tg N/year", description="NOx lightning emission") +NOx_lgt_2d = NOx_lgt_2d.assign_attrs(units="Tg N/year", description="NOx lightning emission") +ds1 = NOx_acf_1d.to_dataset(name='NOx_acf_1d') +ds2 = NOx_acf_2d.to_dataset(name='NOx_acf_2d') +ds3 = NOx_lgt_1d.to_dataset(name='NOx_lgt_1d') +ds4 = NOx_lgt_2d.to_dataset(name='NOx_lgt_2d') +ds = xr.merge([ds1, ds2, ds3, ds4]) +ds.to_netcdf(pathout+'E3SM_NOx_emission_'+startyear+'-'+endyear+'.nc') +# ----- plotting ----- fig = plt.figure(figsize=(18,12)) ax1 = fig.add_subplot(221) plt.plot(NOx_lgt_1d,lev) #levels=np.arange(0, 1000, 20), diff --git a/chemdyg/templates/chemdyg_o3_hole_diags.bash b/chemdyg/templates/chemdyg_o3_hole_diags.bash index 91d480d..1d27be5 100644 --- a/chemdyg/templates/chemdyg_o3_hole_diags.bash +++ b/chemdyg/templates/chemdyg_o3_hole_diags.bash @@ -167,24 +167,39 @@ for i in range(startindex,endindex): else: O3_area[i] = 0. -ds1 = O3_area.to_dataset(name='O3_area') -ds2 = TOZ_min.to_dataset(name='TOZ_min') -ds1.to_netcdf(pathout+'E3SM_O3_area_${y1}-${y2}.nc') -ds2.to_netcdf(pathout+'E3SM_TOZ_min_${y1}-${y2}.nc') - O3_area_time = O3_area[startindex:endindex].sel(time=O3_area[startindex:endindex].time.dt.month.isin([7, 8, 9, 10, 11, 12])) TOZ_min_time = TOZ_min[startindex:endindex].sel(time=TOZ_min[startindex:endindex].time.dt.month.isin([7, 8, 9, 10, 11, 12])) -#-------- climo plot --------------- O3_array = O3_area_time.values.reshape((years,184)) -#O3_array[O3_array==0] = 'nan' O3_mean = O3_array.mean(axis=0) *1.e-12 O3_std = O3_array.std(axis=0) *1.e-12 TOZ_array = TOZ_min_time.values.reshape((years,184)) TOZ_mean = TOZ_array.mean(axis=0) TOZ_std = TOZ_array.std(axis=0) -npdate = np.array(time_range[181:365]) +# ----- writiing ncfile ----- +oz_avg_xr = oz_avg.to_dataset(name='oz_avg') +oz_std_xr = oz_std.to_dataset(name='oz_std') +area_avg_xr = area_avg.to_dataset(name='area_avg') +area_std_xr = area_std.to_dataset(name='area_std') +TOZ_mean_xr = xr.DataArray(TOZ_mean, coords=[oz_avg['time']], dims=["time"]) +TOZ_mean_xr = TOZ_mean_xr.assign_attrs(units="DU", description='Climatological mean of daily SH minimum E3SM total column ozone') +TOZ_std_xr = xr.DataArray(TOZ_std, coords=[oz_avg['time']], dims=["time"]) +TOZ_std_xr = TOZ_std_xr.assign_attrs(units="DU", description='Standard deviation of daily SH minimum E3SM total column ozone') +TOZ_mean_xr = TOZ_mean_xr.to_dataset(name='TOZ_mean') +TOZ_std_xr = TOZ_std_xr.to_dataset(name='TOZ_std') +O3_mean_xr = xr.DataArray(O3_mean, coords=[oz_avg['time']], dims=["time"]) +O3_mean_xr = O3_mean_xr.assign_attrs(units="million of km2", description='Climatological mean of daily mean E3SM O3 hole area') +O3_std_xr = xr.DataArray(O3_std, coords=[oz_avg['time']], dims=["time"]) +O3_std_xr = O3_std_xr.assign_attrs(units="million of km2", description='Standard deviation of daily mean O3 hole area') +O3_mean_xr = O3_mean_xr.to_dataset(name='O3_mean') +O3_std_xr = O3_std_xr.to_dataset(name='O3_std') +ds = xr.merge([oz_avg_xr, oz_std_xr, TOZ_mean_xr, TOZ_std_xr, + area_avg_xr, area_std_xr, O3_mean_xr, O3_std_xr]) +ds.to_netcdf(pathout+'E3SM_O3_hole_${y1}-${y2}.nc') + +#-------- climo plot --------------- +npdate = np.array(time_range[182:366]) fig, ax = plt.subplots(figsize=(10, 5)) plt.plot(npdate,oz_avg, label ='Obs.') plt.fill_between(npdate,oz_avg+oz_std,oz_avg-oz_std, alpha=.5, linewidth=0) diff --git a/chemdyg/templates/chemdyg_pres_lat_plots.bash b/chemdyg/templates/chemdyg_pres_lat_plots.bash index a62598d..0778335 100644 --- a/chemdyg/templates/chemdyg_pres_lat_plots.bash +++ b/chemdyg/templates/chemdyg_pres_lat_plots.bash @@ -148,10 +148,6 @@ o3_2d = o3.mean(axis=2) o3_refer_2d = o3_refer.mean(axis=2) o3_new = o3_2d.copy() o3_refer_new = o3_refer_2d.copy() -ds1 = o3_new.to_dataset() -ds2 = o3_refer_new.to_dataset(name="O3_v2") -ds = xr.merge([ds1, ds2],compat='override') -ds.to_netcdf(pathout+'E3SM_O3_PLatcomparison_'+startyear+'-'+endyear+'.nc') Q = file_in['Q'][0,:,:,:]*28.96/18 Q_refer = refer_in['Q'][0,:,:,:]*28.96/18 @@ -159,10 +155,6 @@ Q_2d = Q.mean(axis=2) Q_refer_2d = Q_refer.mean(axis=2) Q_new = Q_2d.copy() Q_refer_new = Q_refer_2d.copy() -ds1 = Q_new.to_dataset() -ds2 = Q_refer_new.to_dataset(name="Q_v2") -ds = xr.merge([ds1, ds2],compat='override') -ds.to_netcdf(pathout+'E3SM_Q_PLatcomparison_'+startyear+'-'+endyear+'.nc') T = file_in['T'][0,:,:,:] T_refer = refer_in['T'][0,:,:,:] @@ -170,10 +162,6 @@ T_2d = T.mean(axis=2) T_refer_2d = T_refer.mean(axis=2) T_new = T_2d.copy() T_refer_new = T_refer_2d.copy() -ds1 = T_new.to_dataset() -ds2 = T_refer_new.to_dataset(name="T_v2") -ds = xr.merge([ds1, ds2],compat='override') -ds.to_netcdf(pathout+'E3SM_Temp_PLatcomparison_'+startyear+'-'+endyear+'.nc') theda = T*(100000/lev)**0.286 theda_refer = T_refer*(100000/lev)**0.286 @@ -181,10 +169,6 @@ theda_2d = theda.mean(axis=2) theda_refer_2d = theda_refer.mean(axis=2) theda_new = theda_2d.copy() theda_refer_new = theda_refer_2d.copy() -ds1 = theda_new.to_dataset(name="theda") -ds2 = theda_refer_new.to_dataset(name="theda_v2") -ds = xr.merge([ds1, ds2],compat='override') -ds.to_netcdf(pathout+'E3SM_theda_PLatcomparison_'+startyear+'-'+endyear+'.nc') for i in range(180): f = interpolate.interp1d(lev_2d[:,i], o3_2d[:,i]) @@ -227,7 +211,54 @@ Tdiff_relate = Tdiff/T_refer_new thedadiff = theda_new - theda_refer_new thedadiff_relate = thedadiff/theda_refer_new -# plotting +# ----- writing ncffiles ----- +o3_new = o3_new.assign_attrs(units="ppt", description='E3SM O3 concentration') +o3_new_xr = o3_new.to_dataset(name = 'o3_new') +o3_refer_new = o3_refer_new.assign_attrs(units="ppt", description='E3SMv2 O3 concentration') +o3_refer_new_xr = o3_refer_new.to_dataset(name="o3_refer_new") +diff = diff.assign_attrs(units="ppt", description='O3 concentration differences between E3SM and E3SMv2') +diff_xr = diff.to_dataset(name = 'diff') +diff_relate = diff_relate.assign_attrs(description='O3 concentration relative differences between E3SM and E3SMv2') +diff_relate_xr = diff.to_dataset(name = 'diff_relate') +theda_new = theda_new.assign_attrs(units="K", description='E3SM potential temperature') +theda_new_xr = theda_new.to_dataset(name="theda_new") +theda_refer_new = theda_refer_new.assign_attrs(units="K", description='E3SMv2 potential temperature') +theda_refer_new_xr = theda_refer_new.to_dataset(name="theda_refer_new") +thedadiff = thedadiff.assign_attrs(units="K", description='Potential temperature differences between E3SM and E3SMv2') +thedadiff_xr = thedadiff.to_dataset(name = 'thedadiff') +thedadiff_relate = thedadiff_relate.assign_attrs(description='Potential temperature relative differences between E3SM and E3SMv2') +thedadiff_relate_xr = thedadiff_relate.to_dataset(name = 'thedadiff_relate') +T_new = T_new.assign_attrs(units="K", description='E3SM temperature') +T_new_xr = T_new.to_dataset(name="T_new") +T_refer_new = T_refer_new.assign_attrs(units="K", description='E3SMv2 temperature') +T_refer_new_xr = T_refer_new.to_dataset(name="T_refer_new") +Tdiff = Tdiff.assign_attrs(units="K", description='Temperature differences between E3SM and E3SMv2') +Tdiff_xr = Tdiff.to_dataset(name = 'Tdiff') +Tdiff_relate = Tdiff_relate.assign_attrs(description='Temperature relative differences between E3SM and E3SMv2') +Tdiff_relate_xr = Tdiff_relate.to_dataset(name = 'Tdiff_relate') +Q_new = Q_new.assign_attrs(units="ppt", description='E3SM specific humidity') +Q_new_xr = Q_new.to_dataset(name="Q_new") +Q_refer_new = Q_refer_new.assign_attrs(units="ppt", description='E3SMv2 specific humidity') +Q_refer_new_xr = Q_refer_new.to_dataset(name="Q_refer_new") +Qdiff = Qdiff.assign_attrs(units="ppt", description='Specific humidity differences between E3SM and E3SMv2') +Qdiff_xr = Qdiff.to_dataset(name = 'Qdiff') +Qdiff_relate = Qdiff_relate.assign_attrs(description='Specific humidity relative differences between E3SM and E3SMv2') +Qdiff_relate_xr = Qdiff_relate.to_dataset(name = 'Qdiff_relate') +tpp = tpp.assign_attrs(units="Pa", description='Tropopause Pressure') +tpp_xr = tpp.to_dataset(name = 'tpp') +tpp_3d = tpp_3d.assign_attrs(units="Pa", description='Tropopause Pressure (E90 3D)') +tpp_3d_xr = tpp_3d.to_dataset(name = 'tpp_3d') +tpp_refer = tpp_refer.assign_attrs(units="Pa", description='E3SMv2 Tropopause Pressure') +tpp_refer_xr = tpp_refer.to_dataset(name = 'tpp_refer') + +ds = xr.merge([o3_new_xr, o3_refer_new_xr, diff_xr, diff_relate_xr, + theda_new_xr, theda_refer_new_xr, thedadiff_xr, thedadiff_relate_xr, + T_new_xr, T_refer_new_xr, Tdiff_xr, Tdiff_relate_xr, + Q_new_xr, Q_refer_new_xr, Qdiff_xr, Qdiff_relate_xr, + tpp_xr, tpp_3d_xr, tpp_refer_xr],compat='override') +ds.to_netcdf(pathout+'E3SM_Pres_Lat_comparison_'+startyear+'-'+endyear+'.nc') + +# ----- plotting ----- fig = plt.figure(figsize=(18,12)) plt.subplot(2, 2, 1) diff --git a/chemdyg/templates/chemdyg_surf_o3_diags.bash b/chemdyg/templates/chemdyg_surf_o3_diags.bash index 2adb8f2..c629e41 100644 --- a/chemdyg/templates/chemdyg_surf_o3_diags.bash +++ b/chemdyg/templates/chemdyg_surf_o3_diags.bash @@ -200,11 +200,6 @@ O3local_US[O3local_US == 0.] = 'nan' # ## calculate diurnal cycle O3EU_xr = xr.DataArray(O3local_EU, coords=[time_EU,lat_EU,lon_EU], dims=["time","lat","lon"]) O3US_xr = xr.DataArray(O3local_US, coords=[time_US,lat_US,lon_US], dims=["time","lat","lon"]) -ds1 = O3EU_xr.to_dataset(name='sfcO3_EU') -ds2 = O3US_xr.to_dataset(name='sfcO3_US') -ds = xr.merge([ds1, ds2]) -ds.to_netcdf(pathout+'E3SM_sfcO3_${y1}-${y2}.nc') - O3EU_JJA = O3EU_xr.sel(time=O3EU_xr.time.dt.month.isin([6, 7, 8])) O3EU_DJF = O3EU_xr.sel(time=O3EU_xr.time.dt.month.isin([12, 1, 2])) O3US_JJA = O3US_xr.sel(time=O3US_xr.time.dt.month.isin([6, 7, 8])) @@ -287,18 +282,43 @@ for i in range(len(O3EU_DJF)): ndays = int(len(O3EU_JJA)/24) n_d = int(ndays*24) -O3_ENA_JJA_24h = O3_ENA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_WNA_JJA_24h = O3_WNA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_NEU_JJA_24h = O3_NEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_SEU_JJA_24h = O3_SEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_ENA_JJA_24h = 1.e9*O3_ENA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_WNA_JJA_24h = 1.e9*O3_WNA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_NEU_JJA_24h = 1.e9*O3_NEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_SEU_JJA_24h = 1.e9*O3_SEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0) ndays = int(len(O3EU_DJF)/24) n_d = int(ndays*24) -O3_ENA_DJF_24h = O3_ENA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_WNA_DJF_24h = O3_WNA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_NEU_DJF_24h = O3_NEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) -O3_SEU_DJF_24h = O3_SEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) - +O3_ENA_DJF_24h = 1.e9*O3_ENA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_WNA_DJF_24h = 1.e9*O3_WNA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_NEU_DJF_24h = 1.e9*O3_NEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) +O3_SEU_DJF_24h = 1.e9*O3_SEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0) + +# ----- writing ncfile ----- +O3_ENA_JJA_24h_xr = xr.DataArray(O3_ENA_JJA_24h, name= 'O3_ENA_JJA_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="JJA O3 abundance in ENA") ) +O3_WNA_JJA_24h_xr = xr.DataArray(O3_WNA_JJA_24h, name= 'O3_WNA_JJA_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="JJA O3 abundance in WNA") ) +O3_NEU_JJA_24h_xr = xr.DataArray(O3_NEU_JJA_24h, name= 'O3_NEU_JJA_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="JJA O3 abundance in NEU") ) +O3_SEU_JJA_24h_xr = xr.DataArray(O3_SEU_JJA_24h, name= 'O3_SEU_JJA_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="JJA O3 abundance in SEU") ) +O3_ENA_DJF_24h_xr = xr.DataArray(O3_ENA_DJF_24h, name= 'O3_ENA_DJF_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="DJF O3 abundance in ENA") ) +O3_WNA_DJF_24h_xr = xr.DataArray(O3_WNA_DJF_24h, name= 'O3_WNA_DJF_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="DJF O3 abundance in WNA") ) +O3_NEU_DJF_24h_xr = xr.DataArray(O3_NEU_DJF_24h, name= 'O3_NEU_DJF_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="DJF O3 abundance in NEU") ) +O3_SEU_DJF_24h_xr = xr.DataArray(O3_SEU_DJF_24h, name= 'O3_SEU_DJF_24h', coords=[np.arange(0,24)], dims=["hour"], + attrs=dict(units="ppb", description="DJF O3 abundance in SEU") ) +ds1 = O3_ENA_JJA_24h_xr.to_dataset() +ds2 = O3_WNA_JJA_24h_xr.to_dataset() +ds3 = O3_NEU_JJA_24h_xr.to_dataset() +ds4 = O3_SEU_JJA_24h_xr.to_dataset() +ds5 = O3_ENA_DJF_24h_xr.to_dataset() +ds6 = O3_WNA_DJF_24h_xr.to_dataset() +ds7 = O3_NEU_DJF_24h_xr.to_dataset() +ds8 = O3_SEU_DJF_24h_xr.to_dataset() # ## read observations obs_WNA_JJA = [30.5009, 29.4010, 28.0499, 26.9680, 25.3009, 23.5642, 24.1381, 28.3088, 33.9598, 39.3651, @@ -422,28 +442,28 @@ fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=4, sharex=True, figsize=(24, 5)) ax0.plot(obs_WNA_JJA,'ko-') ax0.plot(model_WNA_JJA) -ax0.plot(O3_WNA_JJA_24h*1.e9) +ax0.plot(O3_WNA_JJA_24h) ax0.set_ylim(0,100) ax0.set_ylabel('JJA O3 abundance (ppb)',fontsize='x-large') ax0.set_xlabel('WNA local time (hr)',fontsize='x-large') ax1.plot(obs_ENA_JJA,'ko-') ax1.plot(model_ENA_JJA) -ax1.plot(O3_ENA_JJA_24h*1.e9) +ax1.plot(O3_ENA_JJA_24h) ax1.set_ylim(0,100) #ax1.set_ylabel('JJA O3 abundance (ppb)',fontsize='x-large') ax1.set_xlabel('ENA local time (hr)',fontsize='x-large') ax2.plot(obs_SEU_JJA,'ko-') ax2.plot(model_SEU_JJA) -ax2.plot(O3_SEU_JJA_24h*1.e9) +ax2.plot(O3_SEU_JJA_24h) ax2.set_ylim(0,100) #ax2.set_ylabel('JJA O3 abundance (ppb)',fontsize='x-large') ax2.set_xlabel('SEU local time (hr)',fontsize='x-large') ax3.plot(obs_NEU_JJA,'ko-') ax3.plot(model_NEU_JJA) -ax3.plot(O3_NEU_JJA_24h*1.e9) +ax3.plot(O3_NEU_JJA_24h) ax3.set_ylim(0,100) #ax3.set_ylabel('JJA O3 abundance (ppb)',fontsize='x-large') ax3.set_xlabel('NEU local time (hr)',fontsize='x-large') @@ -573,28 +593,28 @@ fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=4, sharex=True, figsize=(24, 5)) ax0.plot(obs_WNA_DJF,'ko-') ax0.plot(model_WNA_DJF) -ax0.plot(O3_WNA_DJF_24h*1.e9) +ax0.plot(O3_WNA_DJF_24h) ax0.set_ylim(0,100) ax0.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax0.set_xlabel('WNA local time (hr)',fontsize='x-large') ax1.plot(obs_ENA_DJF,'ko-') ax1.plot(model_ENA_DJF) -ax1.plot(O3_ENA_DJF_24h*1.e9) +ax1.plot(O3_ENA_DJF_24h) ax1.set_ylim(0,100) #ax1.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax1.set_xlabel('ENA local time (hr)',fontsize='x-large') ax2.plot(obs_SEU_DJF,'ko-') ax2.plot(model_SEU_DJF) -ax2.plot(O3_SEU_DJF_24h*1.e9) +ax2.plot(O3_SEU_DJF_24h) ax2.set_ylim(0,100) #ax2.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax2.set_xlabel('SEU local time (hr)',fontsize='x-large') ax3.plot(obs_NEU_DJF,'ko-') ax3.plot(model_NEU_DJF) -ax3.plot(O3_NEU_DJF_24h*1.e9) +ax3.plot(O3_NEU_DJF_24h) ax3.set_ylim(0,100) #ax3.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax3.set_xlabel('NEU local time (hr)',fontsize='x-large') @@ -666,10 +686,26 @@ O3_MDA8_WNA_month = np.zeros(12) O3_MDA8_SEU_month = np.zeros(12) O3_MDA8_NEU_month = np.zeros(12) for n in range(12): - O3_MDA8_ENA_month[n] = O3_MDA8_ENA_1D.sel(time=O3_MDA8_ENA_1D.time.dt.month.isin([n+1])).mean() - O3_MDA8_WNA_month[n] = O3_MDA8_WNA_1D.sel(time=O3_MDA8_WNA_1D.time.dt.month.isin([n+1])).mean() - O3_MDA8_NEU_month[n] = O3_MDA8_NEU_1D.sel(time=O3_MDA8_NEU_1D.time.dt.month.isin([n+1])).mean() - O3_MDA8_SEU_month[n] = O3_MDA8_SEU_1D.sel(time=O3_MDA8_SEU_1D.time.dt.month.isin([n+1])).mean() + O3_MDA8_ENA_month[n] = 1.e9*O3_MDA8_ENA_1D.sel(time=O3_MDA8_ENA_1D.time.dt.month.isin([n+1])).mean() + O3_MDA8_WNA_month[n] = 1.e9*O3_MDA8_WNA_1D.sel(time=O3_MDA8_WNA_1D.time.dt.month.isin([n+1])).mean() + O3_MDA8_NEU_month[n] = 1.e9*O3_MDA8_NEU_1D.sel(time=O3_MDA8_NEU_1D.time.dt.month.isin([n+1])).mean() + O3_MDA8_SEU_month[n] = 1.e9*O3_MDA8_SEU_1D.sel(time=O3_MDA8_SEU_1D.time.dt.month.isin([n+1])).mean() + +# ----- writing ncfile ----- +O3_MDA8_ENA_month_xr = xr.DataArray(O3_MDA8_ENA_month, name= 'O3_MDA8_ENA_month',coords=[np.arange(1,13)], dims=["month"], + attrs=dict(units="ppb", description="Mean MDA8 O3 in ENA") ) +O3_MDA8_WNA_month_xr = xr.DataArray(O3_MDA8_WNA_month, name= 'O3_MDA8_WNA_month',coords=[np.arange(1,13)], dims=["month"], + attrs=dict(units="ppb", description="Mean MDA8 O3 in WNA") ) +O3_MDA8_NEU_month_xr = xr.DataArray(O3_MDA8_NEU_month, name= 'O3_MDA8_NEU_month',coords=[np.arange(1,13)], dims=["month"], + attrs=dict(units="ppb", description="Mean MDA8 O3 in NEU") ) +O3_MDA8_SEU_month_xr = xr.DataArray(O3_MDA8_SEU_month, name= 'O3_MDA8_SEU_month',coords=[np.arange(1,13)], dims=["month"], + attrs=dict(units="ppb", description="Mean MDA8 O3 in SEU") ) +M_ds1 = O3_MDA8_ENA_month_xr.to_dataset() +M_ds2 = O3_MDA8_WNA_month_xr.to_dataset() +M_ds3 = O3_MDA8_NEU_month_xr.to_dataset() +M_ds4 = O3_MDA8_SEU_month_xr.to_dataset() +ds = xr.merge([ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8, M_ds1, M_ds2, M_ds3, M_ds4]) +ds.to_netcdf(pathout+'E3SM_surf_O3_${y1}-${y2}.nc') obs_MDA8_WNA = [31.0984, 36.6661, 43.0965, 48.3994, 49.6384, 49.6071, 50.3713, 49.5312, 44.6972, 36.8627, 31.6061, 29.5710] model_MDA8_WNA = [ @@ -737,28 +773,28 @@ fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=4, sharex=True, monthlist = np.arange(1,13) ax0.plot(monthlist,obs_MDA8_WNA,'ko-') ax0.plot(monthlist,model_MDA8_WNA) -ax0.plot(monthlist,O3_MDA8_WNA_month*1.e9) +ax0.plot(monthlist,O3_MDA8_WNA_month) ax0.set_ylim(0,100) ax0.set_ylabel('Mean MDA8 O3 (ppb)',fontsize='x-large') ax0.set_xlabel('WNA (month)',fontsize='x-large') ax1.plot(monthlist,obs_MDA8_ENA,'ko-') ax1.plot(monthlist,model_MDA8_ENA) -ax1.plot(monthlist,O3_MDA8_ENA_month*1.e9) +ax1.plot(monthlist,O3_MDA8_ENA_month) ax1.set_ylim(0,100) #ax1.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax1.set_xlabel('ENA (month)',fontsize='x-large') ax2.plot(monthlist,obs_MDA8_SEU,'ko-') ax2.plot(monthlist,model_MDA8_SEU) -ax2.plot(monthlist,O3_MDA8_SEU_month*1.e9) +ax2.plot(monthlist,O3_MDA8_SEU_month) ax2.set_ylim(0,100) #ax2.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax2.set_xlabel('SEU (month))',fontsize='x-large') ax3.plot(monthlist,obs_MDA8_NEU,'ko-') ax3.plot(monthlist,model_MDA8_NEU) -ax3.plot(monthlist,O3_MDA8_NEU_month*1.e9) +ax3.plot(monthlist,O3_MDA8_NEU_month) ax3.set_ylim(0,100) #ax3.set_ylabel('DJF O3 abundance (ppb)',fontsize='x-large') ax3.set_xlabel('NEU (month)',fontsize='x-large') diff --git a/chemdyg/templates/chemdyg_temperature_eq_plot.bash b/chemdyg/templates/chemdyg_temperature_eq_plot.bash index f7305b8..d7b4b3d 100644 --- a/chemdyg/templates/chemdyg_temperature_eq_plot.bash +++ b/chemdyg/templates/chemdyg_temperature_eq_plot.bash @@ -205,14 +205,17 @@ for i in range(0,20,2): temp_lat_l[ii] = temp_l.mean() latlist = np.arange(60,80,2) -temp_lat_h_xr = xr.DataArray(temp_lat_h, name = 'temp_lat_h',coords=[latlist], dims=["lat"]) -temp_lat_m_xr = xr.DataArray(temp_lat_m, name = 'temp_lat_m',coords=[latlist], dims=["lat"]) -temp_lat_l_xr = xr.DataArray(temp_lat_l, name = 'temp_lat_l',coords=[latlist], dims=["lat"]) +temp_lat_h_xr = xr.DataArray(temp_lat_h, name = 'temp_lat_h',coords=[latlist], dims=["lat"], + attrs=dict(units='K', description='Mean temp.(Jul. to Dec.) with equivalent latitude at 25km height')) +temp_lat_m_xr = xr.DataArray(temp_lat_m, name = 'temp_lat_m',coords=[latlist], dims=["lat"], + attrs=dict(units='K', description='Mean temp.(Jul. to Dec.) with equivalent latitude at 20km height')) +temp_lat_l_xr = xr.DataArray(temp_lat_l, name = 'temp_lat_l',coords=[latlist], dims=["lat"], + attrs=dict(units='K', description='Mean temp.(Jul. to Dec.) with equivalent latitude at 14km height')) ds1 = temp_lat_h_xr.to_dataset() ds2 = temp_lat_m_xr.to_dataset() ds3 = temp_lat_l_xr.to_dataset() ds = xr.merge([ds1, ds2, ds3]) -ds.to_netcdf(pathout+'E3SM_O3_equi_lat_${y1}-${y2}.nc') +ds.to_netcdf(pathout+'E3SM_temp_equi_lat_${y1}-${y2}.nc') fig = plt.figure(figsize=(10,5)) plt.plot(latlist,temp_lat_h) @@ -220,7 +223,7 @@ plt.plot(latlist,temp_lat_m) plt.plot(latlist,temp_lat_l) plt.title('Mean temp.(Jul. to Dec.) with equivalent latitude ${Y1}~${Y2}') plt.xlabel('Lat') -plt.ylabel('Temperature (k)',fontsize='large') +plt.ylabel('Temperature (K)',fontsize='large') plt.legend( ['14 km alt.','20 km alt.','25 km alt.']) pylab.savefig(pathout+'temp_PDF_climo.png', dpi=300)