save plotting data to netcdf files

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'

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`
@@ -226,11 +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())
 
-# 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)
@@ -245,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)
@@ -295,6 +334,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}

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,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_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)')
@@ -165,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)
@@ -210,6 +237,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}

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`
@@ -182,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)
@@ -233,6 +269,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}

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`
@@ -122,6 +130,10 @@ file_in = xr.open_dataset(path+filename)
 
 TMQ = file_in['TMQ'][0,:,:]
 
+# ----- 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))
 
@@ -168,6 +180,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}