diff --git a/chemdyg/templates/chemdyg_index.bash b/chemdyg/templates/chemdyg_index.bash
index 76fcca7..92af929 100644
--- a/chemdyg/templates/chemdyg_index.bash
+++ b/chemdyg/templates/chemdyg_index.bash
@@ -111,7 +111,7 @@ import pandas as pd
pathout = './'
indexfile = open(pathout+'index.html',"w")
-index = ' E3SM chem diagnostics package v1.1 (September 2022) '
+index = ' E3SM chemistry diagnostics package (ChemDyg) '
index = index + ' Test: ${short} '
index = index + ' Reference: Observations and Reanalysis
'
index = index + '
'
diff --git a/chemdyg/templates/chemdyg_noaa_co_comparison.bash b/chemdyg/templates/chemdyg_noaa_co_comparison.bash
index 67e9378..529ad32 100644
--- a/chemdyg/templates/chemdyg_noaa_co_comparison.bash
+++ b/chemdyg/templates/chemdyg_noaa_co_comparison.bash
@@ -184,38 +184,51 @@ for n in range(len(Sta)):
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
nmonth = np.arange(0,len(CO_sel_1D),1)
mask = ~np.isnan(CO_sel_1D)
- if (len(nmonth[mask]) == 0):
- continue
- 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
- diff_noaa = CO_noaa - lin_noaa_xa
- # plotting
- fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5))
- ax1.plot(time_range_noaa,CO_sel_1D,'k')
- ax1.plot(time_range_noaa[mask],lin_e3sm_xa[mask],'k--')
- ax1.plot(time_range_noaa,CO_noaa,'r')
- ax1.plot(time_range_noaa[mask_noaa],lin_noaa_xa[mask_noaa],'r--')
- #ax1.set_ylim(0,200)
- 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'
- 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'
- ax1.legend([line1,line2, line3,line4])
- #ax1.grid(True)
-
- ax2.plot(time_range_noaa[mask],diff[mask],'k')
- ax2.plot(time_range_noaa[mask_noaa],diff_noaa[mask_noaa],'r')
- ax2.set_xlabel('time')
- #ax2.set_ylim(-60,60)
- ax2.set_ylabel('Anomalies')
- pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
+ 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
+ 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'],CO_noaa,'r')
+ ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xa[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'
+ 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'
+ ax1.legend([line1,line2, line3,line4])
+
+ ax2.plot(CO_noaa['time'][mask],diff[mask],'k')
+ ax2.plot(CO_noaa['time'][mask_noaa],diff_noaa[mask_noaa],'r')
+ ax2.set_xlabel('time')
+ #ax2.set_ylim(-60,60)
+ ax2.set_ylabel('Anomalies')
+ pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
+ else:
+ # 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.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'
+ ax1.legend([line3,line4])
+ print('E3SM data is not avaiable within this time period')
+ ax2.plot(CO_noaa['time'][mask_noaa],diff_noaa[mask_noaa],'r')
+ ax2.set_xlabel('time')
+ #ax2.set_ylim(-60,60)
+ ax2.set_ylabel('Anomalies')
+ pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
EOF
diff --git a/chemdyg/templates/chemdyg_o3_hole_diags.bash b/chemdyg/templates/chemdyg_o3_hole_diags.bash
index 265f075..256fe74 100644
--- a/chemdyg/templates/chemdyg_o3_hole_diags.bash
+++ b/chemdyg/templates/chemdyg_o3_hole_diags.bash
@@ -177,25 +177,30 @@ TOZ_array = TOZ_min_time.values.reshape((years,184))
TOZ_mean = TOZ_array.mean(axis=0)
TOZ_std = TOZ_array.std(axis=0)
-fig = plt.figure(figsize=(10,5))
-plt.plot(time_range[181:365],oz_avg, label ='Obs.')
-plt.fill_between(time_range[181:365],oz_avg+oz_std,oz_avg-oz_std, alpha=.5, linewidth=0)
-plt.plot(time_range[181:365],TOZ_mean, label ='E3SM')
-plt.fill_between(time_range[181:365],TOZ_mean+TOZ_std,TOZ_mean-TOZ_std, alpha=.5, linewidth=0)
+npdate = np.array(time_range[181:365])
+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)
+plt.plot(npdate,TOZ_mean, label ='E3SM')
+plt.fill_between(npdate,TOZ_mean+TOZ_std,TOZ_mean-TOZ_std, alpha=.5, linewidth=0)
+date_form = DateFormatter("%b-%d")
+ax.xaxis.set_major_formatter(date_form)
plt.legend(loc = 'upper left')
-plt.title('SH minimum total O3')
-plt.xlabel('Time')
+plt.title('SH minimum total O3 ('+str(startyear)+' - '+str(endyear)+')')
+plt.xlabel('Date')
plt.ylabel('O3 conc. (DU)',fontsize='large')
pylab.savefig(pathout+'to3mins_O3hole.png', dpi=600)
-fig = plt.figure(figsize=(10,5))
-plt.plot(time_range[181:365],area_avg, label ='Obs.')
-plt.fill_between(time_range[181:365],area_avg+area_std,area_avg-area_std, alpha=.5, linewidth=0)
-plt.plot(time_range[181:365],O3_mean, label ='E3SM')
-plt.fill_between(time_range[181:365],O3_mean+O3_std,O3_mean-O3_std, alpha=.5, linewidth=0)
+fig, ax = plt.subplots(figsize=(10, 5))
+plt.plot(npdate,area_avg, label ='Obs.')
+plt.fill_between(npdate,area_avg+area_std,area_avg-area_std, alpha=.5, linewidth=0)
+plt.plot(npdate,O3_mean, label ='E3SM')
+plt.fill_between(npdate,O3_mean+O3_std,O3_mean-O3_std, alpha=.5, linewidth=0)
+date_form = DateFormatter("%b-%d")
+ax.xaxis.set_major_formatter(date_form)
plt.legend(loc = 'upper left')
-plt.title('O3 hole area (million of km2)')
-plt.xlabel('Time')
+plt.title('O3 hole area (million of km2) ('+str(startyear)+' - '+str(endyear)+')')
+plt.xlabel('Date')
plt.ylabel('Area (million of km2)',fontsize='large')
pylab.savefig(pathout+'to3areas_O3hole.png', dpi=600)
diff --git a/chemdyg/templates/chemdyg_summary_table.bash b/chemdyg/templates/chemdyg_summary_table.bash
index 96152f7..2c60679 100644
--- a/chemdyg/templates/chemdyg_summary_table.bash
+++ b/chemdyg/templates/chemdyg_summary_table.bash
@@ -218,8 +218,6 @@ for var in range(len(varname)):
# write out annual chem tendency
line_ann = line_ann + ' '+ format('O3 burden (Tg)','28s')
line_ann = line_ann + ' '+"{0:+.3e}".format(np.array(TOZ_total)) +''
- line_ann = line_ann + ' '+ format('O3 emission (Tg/yr)','28s')
- line_ann = line_ann + ' '+"{0:+.3e}".format(np.array(TDS_total)*unit_covet) +''
line_ann = line_ann + ' '+ format('O3 deposition (Tg/yr)','28s')
line_ann = line_ann + ' '+"{0:+.3e}".format(np.array(TDD_total)*unit_covet) +''
line_ann = line_ann + ' '+ format('O3 production (Tg/yr)','28s')
diff --git a/chemdyg/templates/chemdyg_surf_o3_diags.bash b/chemdyg/templates/chemdyg_surf_o3_diags.bash
index 27c9797..d86869c 100644
--- a/chemdyg/templates/chemdyg_surf_o3_diags.bash
+++ b/chemdyg/templates/chemdyg_surf_o3_diags.bash
@@ -275,16 +275,18 @@ for i in range(len(O3EU_DJF)):
O3_SEU_DJF[i] = (O3EU_DJF_SEU[i,:,:]*area_SEU[:,:]).sum()/area_SEU[:,:].sum()
ndays = int(len(O3EU_JJA)/24)
-O3_ENA_JJA_24h = O3_ENA_JJA.reshape((ndays,24)).mean(axis=0)
-O3_WNA_JJA_24h = O3_WNA_JJA.reshape((ndays,24)).mean(axis=0)
-O3_NEU_JJA_24h = O3_NEU_JJA.reshape((ndays,24)).mean(axis=0)
-O3_SEU_JJA_24h = O3_SEU_JJA.reshape((ndays,24)).mean(axis=0)
+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)
ndays = int(len(O3EU_DJF)/24)
-O3_ENA_DJF_24h = O3_ENA_DJF.reshape((ndays,24)).mean(axis=0)
-O3_WNA_DJF_24h = O3_WNA_DJF.reshape((ndays,24)).mean(axis=0)
-O3_NEU_DJF_24h = O3_NEU_DJF.reshape((ndays,24)).mean(axis=0)
-O3_SEU_DJF_24h = O3_SEU_DJF.reshape((ndays,24)).mean(axis=0)
+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)
# ## read observations
@@ -614,8 +616,8 @@ for d in range(ndays):
O3US_MDA8 = O3US_navg.max(axis=1)
O3EU_MDA8 = O3EU_navg.max(axis=1)
-O3EU_MDA8_xr = xr.DataArray(O3EU_MDA8, coords=[day_timearray['time'],lat_EU,lon_EU], dims=["time","lat","lon"])
-O3US_MDA8_xr = xr.DataArray(O3US_MDA8, coords=[day_timearray['time'],lat_US,lon_US], dims=["time","lat","lon"])
+O3EU_MDA8_xr = xr.DataArray(O3EU_MDA8, coords=[day_timearray['time'][0:ndays],lat_EU,lon_EU], dims=["time","lat","lon"])
+O3US_MDA8_xr = xr.DataArray(O3US_MDA8, coords=[day_timearray['time'][0:ndays],lat_US,lon_US], dims=["time","lat","lon"])
# US subregion: Western (WNA) and eastern (ENA) North America is split by 264°E
lon1 = 234
@@ -636,16 +638,16 @@ O3_MDA8_ENA_1D = np.zeros(ndays)
O3_MDA8_WNA_1D = np.zeros(ndays)
O3_MDA8_NEU_1D = np.zeros(ndays)
O3_MDA8_SEU_1D = np.zeros(ndays)
-for i in range(len(day_timearray)):
+for i in range(ndays):
O3_MDA8_ENA_1D[i] = (O3_MDA8_ENA[i,:,:]*area_ENA[:,:]).sum()/area_ENA[:,:].sum()
O3_MDA8_WNA_1D[i] = (O3_MDA8_WNA[i,:,:]*area_WNA[:,:]).sum()/area_WNA[:,:].sum()
O3_MDA8_NEU_1D[i] = (O3_MDA8_NEU[i,:,:]*area_NEU[:,:]).sum()/area_NEU[:,:].sum()
O3_MDA8_SEU_1D[i] = (O3_MDA8_SEU[i,:,:]*area_SEU[:,:]).sum()/area_SEU[:,:].sum()
-O3_MDA8_ENA_1D = xr.DataArray(O3_MDA8_ENA_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_WNA_1D = xr.DataArray(O3_MDA8_WNA_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_NEU_1D = xr.DataArray(O3_MDA8_NEU_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_SEU_1D = xr.DataArray(O3_MDA8_SEU_1D, coords=[day_timearray['time']], dims=["time"])
+O3_MDA8_ENA_1D = xr.DataArray(O3_MDA8_ENA_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_WNA_1D = xr.DataArray(O3_MDA8_WNA_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_NEU_1D = xr.DataArray(O3_MDA8_NEU_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_SEU_1D = xr.DataArray(O3_MDA8_SEU_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
# calculate monthly mean
O3_MDA8_ENA_month = np.zeros(12)