diff --git a/PE_CBL.py b/PE_CBL.py
index faff346e82210537aa34a0b1377ddef9a663a0cb..e37393f049d3d61ef810a91503c29ae4183e3794 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -71,6 +71,7 @@ if __name__ == '__main__':
# In case of parameter estimation
'do_parameter_estimation' : True,
'parameter_inflation_rtps_alpha' : np.array([0.8]),
+ 'return_covariances_increments_and_innovations' : True
}
integration_dt = 0.25*default_cbl_settings["dz"]**2/default_cbl_settings["Kmax"]
@@ -95,12 +96,12 @@ if __name__ == '__main__':
"Assimilation interval must be an integer multiplier of model dt"
# Decide what figures to plot
- fig01 = True
+ fig01 = False
fig02 = True
- fig03 = True
- fig04 = True
- fig05 = True
- fig06 = True
+ fig03 = False
+ fig04 = False
+ fig05 = False
+ fig06 = False
fig07 = True
fig08 = True
@@ -140,9 +141,12 @@ if __name__ == '__main__':
cbl_settings_A_noPE['do_parameter_estimation'] = False
da_settings_A_noPE['cbl_settings'] = cbl_settings_A_noPE
- # Run and save to disk
- exp_A_noPE = experiment(da_settings_A_noPE)
- pickle.dump(exp_A_noPE, open('exp_A_noPE.pickle', 'wb'))
+ # Run it
+ try:
+ exp_A_noPE = pickle.load(open("exp_A_noPE.pickle", "rb"))
+ except:
+ exp_A_noPE = experiment(da_settings_A_noPE)
+ pickle.dump(exp_A_noPE, open('exp_A_noPE.pickle', 'wb'))
if fig01:
@@ -182,15 +186,15 @@ if __name__ == '__main__':
# Make plots
ncont = 13
- fig, [ax1, ax2, ax3] = p.subplots(1,3,constrained_layout=True)
- fig.set_size_inches(6,3)
+ fig, [[ax4, ax2],[ax1, ax3]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,6)
c1 = ax1.pcolormesh(cbl_det.history['time']/3600,cbl_det.zt,cbl_det.history['theta'],vmin=290,vmax=296)
ax1.set_ylim([0,zmax])
ax1.set_ylabel(r'Height (m)')
ax1.set_xlabel(r'Time (h)')
ax1.set_xticks(np.arange(4))
- ax1.set_title(r'a) $\theta$ (K)')
+ ax1.set_title(r'c) $\overline{\theta}$ (K)')
p.colorbar(c1,orientation='horizontal')
ax1.contour(cbl_det.history['time']/3600,cbl_det.zt,cbl_det.history['theta'],
np.linspace(cbl_det.theta_0,cbl_det.theta_0+cbl_det.gamma*zmax,ncont),
@@ -199,20 +203,32 @@ if __name__ == '__main__':
linewidths=0.75)
ax2 = plot_p(p_factors,theta_profiles,cbl_pf.zt,None,ax=ax2)
- ax2.set_xlabel(r'$\theta$ (K)')
+ ax2.set_ylabel(r'Height (m)')
+ ax2.set_xlabel(r'$\overline{\theta}$ (K)')
ax2.set_xlim([291,297])
ax2.set_ylim([0,zmax])
ax2.legend(loc=4,frameon=False)
- ax2.set_title(r'b) Sensitivity to $p$')
+ ax2.set_title(r'b) $\overline{\theta}$ sensitivity to $p$')
ax3,c3 = plot_spread(cbl_free,ax=ax3)
- ax3.set_title(r'c) $\sigma_\theta$ (K)')
+ ax3.set_ylabel(r'Height (m)')
+ ax3.set_title(r'd) $\sigma_\theta$ (K)')
ax3.set_xlabel('Time (h)')
ax3.set_xticks(np.arange(4))
p.colorbar(c3,orientation='horizontal')
- p.setp(ax2.get_yticklabels(), visible=False)
- p.setp(ax3.get_yticklabels(), visible=False)
+ zoverh = np.linspace(0,1,101)
+ for pfac in p_factors:
+ Koverkws = zoverh*(1-zoverh)**pfac
+ ax4.plot(Koverkws,zoverh,label='$p=%4.1f$'%pfac)
+ ax4.set_title(r'a) $K_h$ sensitivity to $p$')
+ ax4.set_xlabel('$K_h/(\kappa w_s h)$')
+ ax4.set_ylabel('$z/h$')
+ ax4.set_xlim([0,0.5])
+ ax4.legend(loc=4,frameon=False)
+
+ #p.setp(ax2.get_yticklabels(), visible=False)
+ #p.setp(ax3.get_yticklabels(), visible=False)
fig.savefig('fig01.png',format='png',dpi=300)
p.close(fig)
@@ -223,13 +239,13 @@ if __name__ == '__main__':
fig.set_size_inches(6,6)
#
[ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],None,ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. A, cov($p,y_b}$) (K)')
+ ax0.set_title(r'a) Exp. A, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. A, $\sigma^2_{y^b}}$ (K)')
+ ax1.set_title(r'b) Exp. A, $\delta y\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. A, $K_{p,y_b}$ (K$^{-1}$)')
+ ax2.set_title(r'c) Exp. A, $\delta p\prime\prime$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3 = plot_CBL_PE(exp_A,None,ax=ax3)
@@ -245,56 +261,14 @@ if __name__ == '__main__':
if fig03:
- def plotfig(exprange,filename):
- fig, [[ax1, ax2],[ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,4)
- z = exp_A.obs_coordinates
- z_pbl = z*1.
- z_pbl[z>1000] = np.nan
- for i in exprange:
- i1 = experiments_1[i].dg
- i2 = experiments_2[i].dg
- ax1.plot(i1.aRMSE_t,z,label=labels[i],color=linecolors[i])
- ax1.plot(i2.aRMSE_t,z,color=linecolors[i],dashes=[3,1],alpha=0.3)
- #
- ax2.plot(i1.bRMSE_t,z,label=labels[i],color=linecolors[i])
- ax2.plot(i2.bRMSE_t,z,color=linecolors[i],dashes=[3,1],alpha=0.3)
- #
- ax3.plot(i1.bRMSE_t-i1.aRMSE_t,z,label=labels[i],color=linecolors[i])
- ax3.plot(i2.bRMSE_t-i2.aRMSE_t,z,color=linecolors[i],dashes=[3,1],alpha=0.3)
- #
- ax4.plot(i1.bSprd_t/i1.bRMSE_t,z_pbl,label=labels[i],color=linecolors[i])
- ax4.plot(i2.bSprd_t/i2.bRMSE_t,z_pbl,color=linecolors[i],dashes=[3,1],alpha=0.3)
- ax1.set_title('a) Analysis error')
- ax1.set_xlabel(r'RMSE$^a_\theta$')
- ax2.set_title('b) First-guess error')
- ax2.set_xlabel(r'RMSE$^b_\theta$')
- ax3.set_title('c) Error reduction')
- ax3.set_xlabel(r'RMSE$^b_\theta-$RMSE$^a_\theta$')
- ax4.set_title('d) Spread-error consistency')
- ax4.set_xlabel(r'$\sigma^b_\theta$/RMSE$^b_\theta$')
- ax1.set_ylabel('height (m)')
- ax3.set_ylabel('height (m)')
- #
- #ax2.legend(frameon=False)
- ax4.axvline(x=1,color='k',linewidth=0.5,dashes=[3,1])
- ax2.sharey(ax1)
- ax4.sharey(ax3)
- p.setp(ax2.get_yticklabels(), visible=False)
- p.setp(ax4.get_yticklabels(), visible=False)
- #
- fig.savefig(filename,format='png',dpi=300)
- p.close(fig)
-
exp_A = pickle.load(open("exp_A.pickle", "rb"))
exp_A_noPE = pickle.load(open("exp_A_noPE.pickle", "rb"))
- experiments_1 = [exp_A]
- experiments_2 = [exp_A_noPE]
+ experiments_pe = [exp_A]
+ experiments_nope = [exp_A_noPE]
labels = ["Exp. A"]
- linecolors = p.rcParams['axes.prop_cycle'].by_key()['color']
- plotfig(range(1), 'fig03.png')
+ plot_diagnostics(experiments_pe,experiments_nope,labels,'fig03.png')
if fig04:
@@ -588,7 +562,7 @@ if __name__ == '__main__':
cbl_settings_D['Hmax'] = 0.15
cbl_settings_D['is_cgrad'] = False
cbl_settings_D['simulate_error_growth'] = True
- cbl_settings_D['error_growth_perturbations_amplitude'] = sigma_b_init*5
+ cbl_settings_D['error_growth_perturbations_amplitude'] = sigma_b_init*10
da_settings_D['cbl_settings'] = cbl_settings_D
da_settings_D['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*5
da_settings_D['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
@@ -598,29 +572,18 @@ if __name__ == '__main__':
setattr(exp_D,'label','D')
pickle.dump(exp_D, open('exp_D.pickle', 'wb'))
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_D_noPE = dict(cbl_settings_D)
- da_settings_D_noPE = dict(da_settings_D)
- cbl_settings_D_noPE['do_parameter_estimation'] = False
- da_settings_D_noPE['cbl_settings'] = cbl_settings_D_noPE
-
- # Run and save to disk
- exp_D_noPE = experiment(da_settings_D_noPE)
- pickle.dump(exp_D_noPE, open('exp_D_noPE.pickle', 'wb'))
-
# Make plots
fig, [[ax0, ax1],[ax2,ax3]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,6)
#
[ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],None,ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. D, cov($p,y_b}$) (K)')
+ ax0.set_title(r'a) Exp. D, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. D, $\sigma^2_{y^b}}$ (K)')
+ ax1.set_title(r'b) Exp. D, $\delta y\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. D, $K_{p,y_b}$ (K$^{-1}$)')
+ ax2.set_title(r'c) Exp. D, $\delta p\prime\prime$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3 = plot_CBL_PE(exp_D,None,ax=ax3)
@@ -635,6 +598,44 @@ if __name__ == '__main__':
fig.savefig('fig07.png',format='png',dpi=300)
p.close(fig)
+ if fig08:
+
+ # Create a copy of the default settings
+ cbl_settings_D = dict(default_cbl_settings)
+ da_settings_D = dict(default_da_settings)
+
+ # Change settings as necessary
+ # Changes include generation of observations, so the existing nature run
+ # can't be reused.
+ cbl_settings_D['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
+ cbl_settings_D['perturbations_theta_amplitude'] = sigma_b_init*10
+ cbl_settings_D['Hmax'] = 0.15
+ cbl_settings_D['is_cgrad'] = False
+ cbl_settings_D['simulate_error_growth'] = True
+ cbl_settings_D['error_growth_perturbations_amplitude'] = sigma_b_init*10
+ da_settings_D['cbl_settings'] = cbl_settings_D
+ da_settings_D['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*10
+ da_settings_D['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
+
+ # Experiment matching D, but without parameter estimation
+ cbl_settings_D_noPE = dict(cbl_settings_D)
+ da_settings_D_noPE = dict(da_settings_D)
+ cbl_settings_D_noPE['do_parameter_estimation'] = False
+ da_settings_D_noPE['cbl_settings'] = cbl_settings_D_noPE
+
+ exp_D = pickle.load(open("exp_D.pickle", "rb"))
+ try:
+ exp_D_noPE = pickle.load(open("exp_D_noPE.pickle", "rb"))
+ except:
+ exp_D_noPE = experiment(da_settings_D_noPE)
+ pickle.dump(exp_A_noPE, open('exp_D_noPE.pickle', 'wb'))
+
+ experiments_pe = [exp_D]
+ experiments_nope = [exp_D_noPE]
+ labels = ["Exp. D"]
+
+ plot_diagnostics(experiments_pe,experiments_nope,labels,'fig08.png')
+
if opt01:
da_settings = {'cbl_settings' : dict(default_cbl_settings),