diff --git a/PE_CBL.py b/PE_CBL.py
index 77978b20d1dc077c77754d38ebeacf92b41fca6f..3f18509833be39158857403ebe424f93caf7ba3d 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -56,7 +56,7 @@ if __name__ == '__main__':
fig05 = False # CHECKED OK # Experiments with systematic model errors
fig06 = False # CHECKED OK # Analysis increments in model error experiments
fig07 = False # CHECKED OK # Experiment with error growth
- fig08 = False # CHECKED OK # Perfect-obs experiments
+ fig08 = True # CHECKED OK # Perfect-obs experiments
fig09 = False # CHECKED OK # DA diagnostics (profiles) for perfect-obs experiments
fig10 = False # CHECKED OK # Comparison EAKF vs LETKF
@@ -69,11 +69,13 @@ if __name__ == '__main__':
opt06 = False # CHECKED OK # Check impact of RTPS variance inflation factor for parameters
opt07 = False # CHECKED OK # Check impact of RTPS variance inflation factor for state vector
opt08 = False # CHECKED OK # Check sensitivity to observation sorting
- opt09 = False # CHECKED OK # Plot covariance matrix
- opt10 = False # CHECKED OK # Consistency checks
+ opt09 = False # CHECKED OK # Background error covariance matrix
+ opt10 = False # CHECKED OK # Profiles of consistency diagnostics
opt11 = False # CHECKED OK # An experiment without any variance inflation
opt12 = False # CHECKED OK # Old version of Figure 08
opt13 = False # CHECKED OK # More profiles of consistency diagnostics
+ opt14 = False # CHECKED OK # More plots of consistency checks
+ opt15 = False # CHECKED OK # O-B for experiments C1 and C2
# A test to check if results are exactly reproducible
check_reproducibility = False
@@ -361,7 +363,7 @@ if __name__ == '__main__':
ax2.set_title(r'a) Exp. B$_4$, $\delta\overline{y}$ (K)')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
- ax3.set_title(r'b) Exp. B$_4$, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
+ ax3.set_title(r'b) Exp. B$_4$, $\delta\overline{p}\prime\prime$')
ax3.set_xlabel('Time (h)')
ax3.set_ylabel('Height (m)')
maxtime = exp_B4.trun
@@ -776,6 +778,123 @@ if __name__ == '__main__':
fig.savefig('diagnostics_profiles_errors_vs_deviations.png',format='png',dpi=300)
p.close(fig)
+ if opt14:
+
+ workdir = "./runs/seed=181612_enssize=20_EAKF_6hrs/"
+
+ [exp_I1,exp_I2,exp_I3,exp_I4,
+ exp_I5,exp_I6,exp_I7,exp_I8,] = map(get_experiment,
+ [workdir+'/da_exp_I1.json',
+ workdir+'/da_exp_I2.json',
+ workdir+'/da_exp_A1.json',
+ workdir+'/da_exp_I4.json',
+ workdir+'/da_exp_I5.json',
+ workdir+'/da_exp_A1.json',
+ workdir+'/da_exp_I7.json',
+ workdir+'/da_exp_I8.json'])
+ exp_I3.label = 'I3'
+ exp_I6.label = 'I6'
+
+ # Plots: Scatter plots of Desroziers statistics
+ experiments1 = [exp_I1, exp_I2, exp_I3, exp_I4]
+ experiments2 = [exp_I5, exp_I6, exp_I7, exp_I8]
+ labels1 = [exp.label for exp in experiments1]
+ labels2 = [exp.label for exp in experiments2]
+
+ fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
+ fig.set_size_inches(6,3)
+
+ ax1 = plot_desroziers([i.dg for i in experiments1],labels=labels1,ax=ax1)
+ ax2 = plot_desroziers([i.dg for i in experiments2],labels=labels2,ax=ax2)
+ ax1.set_xlim([0,0.018])
+ ax1.set_ylim([0,0.018])
+ ax2.set_xlim([0,0.018])
+ ax2.set_ylim([0,0.018])
+ ax1.set_title(r"a) Variable $\alpha_{RTPS}$ for parameter")
+ ax2.set_title(r"b) Variable $\alpha_{RTPS}$ for state")
+
+ fig.savefig('Diag_desroziers_1.png',format='png',dpi=150)
+ p.close(fig)
+
+ [exp_L1,exp_L2,exp_L3,exp_L4] = map(get_experiment,
+ [workdir+'/da_exp_L1.json',
+ workdir+'/da_exp_A1.json',
+ workdir+'/da_exp_L3.json',
+ workdir+'/da_exp_L4.json'])
+ exp_L2.label = 'L2'
+ experiments3 = [exp_L1, exp_L2, exp_L3, exp_L4]
+ labels3 = [exp.label for exp in experiments3]
+
+ fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
+ fig.set_size_inches(6,3)
+
+ ax1 = plot_desroziers([i.dg for i in experiments3],labels=labels3,ax=ax1)
+ ax2 = plot_desroziers([i.dg for i in experiments2],labels=labels2,ax=ax2)
+ ax1.set_xlim([0,0.018])
+ ax1.set_ylim([0,0.018])
+ ax2.set_xlim([0,0.018])
+ ax2.set_ylim([0,0.018])
+ ax1.set_title(r"a) Variable localicazion scale")
+ ax2.set_title(r"b) Variable $\alpha_{RTPS}$ for state")
+
+ fig.savefig('Diag_desroziers_2.png',format='png',dpi=150)
+ p.close(fig)
+
+ if opt15:
+
+ [exp_C1,exp_C2] = map(get_experiment, [workdir+'/da_exp_C1.json',
+ workdir+'/da_exp_C2.json'])
+
+ ntimes,nobs = exp_C1.da.observations.shape
+
+ # Expand times array to 2D (matching the number of observations)
+ times = exp_C1.da.time/3600
+ times = times[:,None]+np.zeros((ntimes,nobs))
+
+ # Sort coordinate and innovation arrays in regular order
+ ominb_C1 = np.zeros((ntimes,nobs))+np.nan
+ ominb_C2 = np.zeros((ntimes,nobs))+np.nan
+ z_C1 = np.zeros((ntimes,nobs))+np.nan
+ z_C2 = np.zeros((ntimes,nobs))+np.nan
+ for [z_sorted,innov_sorted,z,innov] in zip([z_C1,z_C2],[ominb_C1,ominb_C2],[exp_C1.da.obs_coordinates,exp_C2.da.obs_coordinates],[exp_C1.da.innov,exp_C2.da.innov]):
+ for k in range(ntimes):
+ if is_sorted(z[k,:]):
+ indices = np.arange(nobs)
+ else:
+ indices = np.argsort(z[k,:])
+ innov_sorted[k,:] = innov[k,indices]
+ z_sorted[k,:] = z[k,indices]
+
+ fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
+ fig.set_size_inches(6,3)
+ c1 = ax1.pcolormesh(times,z_C1,ominb_C1,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
+ c2 = ax2.pcolormesh(times,z_C2,ominb_C2,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
+ ax1.text(0.1,0.2,r'$\mu_d$ = %05.4f'%ominb_C1.mean(), transform = ax1.transAxes)
+ ax1.text(0.1,0.1,r'$\sigma_d$ = %05.4f'%ominb_C1.std(),transform = ax1.transAxes)
+ ax2.text(0.1,0.2,r'$\mu_d$ = %05.4f'%ominb_C2.mean(), transform = ax2.transAxes)
+ ax2.text(0.1,0.1,r'$\sigma_d$ = %05.4f'%ominb_C2.std(),transform = ax2.transAxes)
+ #
+ ax1.set_ylim([0,2000])
+ ax2.set_ylim([0,2000])
+ ax1.set_title(r'a) Exp. C$_1$, $\overline{\theta}^o-\overline{\theta}^b$')
+ ax2.set_title(r'b) Exp. C$_2$, $\overline{\theta}^o-\overline{\theta}^b$')
+ ax1.set_xlabel('Time (h)')
+ ax2.set_xlabel('Time (h)')
+ ax1.set_ylabel('Height (m)')
+ ax2.sharey(ax1)
+ maxtime = exp_C1.trun
+ ax1.set_xlim([0,maxtime/3600])
+ ax2.set_xlim([0,maxtime/3600])
+ ax1.set_xticks(np.arange((maxtime+1)/3600))
+ ax2.set_xticks(np.arange((maxtime+1)/3600))
+ p.setp(ax2.get_yticklabels(), visible=False)
+ #
+ p.colorbar(c1,orientation='horizontal')
+ p.colorbar(c2,orientation='horizontal')
+ #
+ fig.savefig('ominb.png',format='png',dpi=300)
+ p.close(fig)
+
if check_reproducibility:
def get_obs(f):