diff --git a/PE_CBL.py b/PE_CBL.py
index 4acb54260673465a0f832e3b20958e3d2af65396..5db83c7d35027e694825952991c09acd6cf38f10 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 from matplotlib import pyplot as p
+from matplotlib import colors
 import pickle
 import os
 import json
@@ -25,85 +26,45 @@ def load_or_run(settings):
 
 if __name__ == '__main__':
 
+    workdir = "./runs/seed=181612_enssize=200_EAKF_6hrs/"
+
     ##################################################################################
     # Check if working directory exists, if not, create it.
-    workdir = "./runs/seed=181612_enssize=20_EAKF_6hrs/"
     if not os.path.exists(workdir):
         os.makedirs(workdir)
 
     ##################################################################################
-    # Load the experiments or run them
-
-    with open('./config_files/da_exp_A1.json', 'r') as fp:
-        da_settings_A1 = json.load(fp)
-    exp_A1 = load_or_run(da_settings_A1)
-
-    with open('./config_files/da_exp_A2.json', 'r') as fp:
-        da_settings_A2 = json.load(fp)
-    exp_A2 = load_or_run(da_settings_A2)
-
-    with open('./config_files/da_exp_B1.json', 'r') as fp:
-        da_settings_B1 = json.load(fp)
-    exp_B1 = load_or_run(da_settings_B1)
-
-    with open('./config_files/da_exp_B2.json', 'r') as fp:
-        da_settings_B2 = json.load(fp)
-    exp_B2 = load_or_run(da_settings_B2)
-
-    with open('./config_files/da_exp_B3.json', 'r') as fp:
-        da_settings_B3 = json.load(fp)
-    exp_B3 = load_or_run(da_settings_B3)
-
-    with open('./config_files/da_exp_B4.json', 'r') as fp:
-        da_settings_B4 = json.load(fp)
-    exp_B4 = load_or_run(da_settings_B4)
-
-    with open('./config_files/da_exp_C1.json', 'r') as fp:
-        da_settings_C1 = json.load(fp)
-    exp_C1 = load_or_run(da_settings_C1)
-
-    with open('./config_files/da_exp_C2.json', 'r') as fp:
-        da_settings_C2 = json.load(fp)
-    exp_C2 = load_or_run(da_settings_C2)
-
-    with open('./config_files/da_exp_D1.json', 'r') as fp:
-        da_settings_D1 = json.load(fp)
-    exp_D1 = load_or_run(da_settings_D1)
-
-    with open('./config_files/da_exp_D2.json', 'r') as fp:
-        da_settings_D2 = json.load(fp)
-    exp_D2 = load_or_run(da_settings_D2)
-
-    noPE_runs = False
-    if noPE_runs:
-        with open('./config_files/da_exp_B1_noPE.json', 'r') as fp:
-            da_settings_B1_noPE = json.load(fp)
-        exp_B1_noPE = load_or_run(da_settings_B1_noPE)
+    # Check if the necessary config files exist in the working directory.
+    # If not, copy them.
+    # for src in args[:-1]:
+    #   os.symlink(src, os.path.join(dst, os.path.dirname(src)))
 
     ########################################################################
     # Manuscript figures
     fig01 = False # CHECKED OK # Sensitivity to p
     fig02 = False # CHECKED OK # Successful PE experiment
-    fig03 = False # CHECKED OK # DA diagnostics
+    fig03 = False # CHECKED OK # Consistency checks
     fig04 = False # CHECKED OK # Experiments with varying weights
     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 # Comparison successful vs unsuccessful PE
-    fig09 = False # Comparison EAKF vs LETKF
+    fig09 = True  # CHECKED OK # Comparison EAKF vs LETKF
+    fig10 = True  # CHECKED OK # DA diagnostics (profiles) for perfect-obs experiments
     
     # Other figures
     opt01 = False # CHECKED OK # assimilation of a single observation
     opt02 = False # CHECKED OK # assimilation of profiles at two times
-    opt03 = False # CHECKED OK # NOT EXPECTED TO WORK # assimilation of real observations
-    opt04 = True # NOT EXPECTED TO WORK # assimilation of real observations, different spinup times
+    opt03 = False # CHECKED OK # assimilation of real observations
+    opt04 = False # CHECKED OK # assimilation of real observations, different spinup times
     opt05 = False # CHECKED OK # Check impact of localization cutoff
     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
-    opt11 = False # CHECKED OK # Consistency checks
+    opt11 = False # CHECKED OK # An experiment without any variance inflation
+    opt12 = False # CHECKED OK # Old version of Figure 08
 
     # A test to check if results are exactly reproducible
     check_reproducibility = False
@@ -112,7 +73,7 @@ if __name__ == '__main__':
 
         # Create a copy of the default settings
         # and update as necessary
-        with open('./config_files/default_cbl.json', 'r') as fp:
+        with open(workdir+'/default_cbl.json', 'r') as fp:
             cbl_settings = json.load(fp)
         maxtime = 21600
         cbl_settings['maxtime'] = maxtime
@@ -140,7 +101,7 @@ if __name__ == '__main__':
         # Panel c) spread induced by p
         # Do a free ensemble run (ensemble size set expliclity)
         # Set size to 200 members
-        with open('./config_files/default_cbl.json', 'r') as fp:
+        with open(workdir+'/default_cbl.json', 'r') as fp:
             cbl_settings_free = json.load(fp)
         cbl_settings_free['maxtime'] = maxtime
         cbl_settings_free['perturb_ensemble_state'] = False
@@ -191,9 +152,9 @@ if __name__ == '__main__':
         ax4.set_xlim([0,0.5])
         ax4.legend(loc=4,frameon=False)
 
-        ax5,c5 = plot_free_ensemble(cbl_free,show='cov',ax=ax5,cmap='RdBu_r')
+        ax5,c5 = plot_free_ensemble(cbl_free,show='corr',ax=ax5,cmap='RdBu_r')
         ax5.set_ylabel(r'Height (m)')
-        ax5.set_title(r'e) cov$(p,\overline{\theta})$ (K)')
+        ax5.set_title(r'e) $\rho(p,\theta\prime)$')
         ax5.set_xlabel('Time (h)')
         ax5.set_xticks(np.arange((maxtime+1)/3600))
         p.colorbar(c5,orientation='horizontal')
@@ -210,62 +171,101 @@ if __name__ == '__main__':
 
     if fig02:
 
-        exp_A1 = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A1 = load_or_run(da_settings_A1)
+
         plot_overview(exp_A1,"A$_1$",'fig02.png')
 
     if fig03:
 
-        exp1 = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
-        exp2 = pickle.load(open(workdir+"exp_A2.pickle", "rb"))
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A1 = load_or_run(da_settings_A1)
 
-        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
-        fig.set_size_inches(6,4)
-        plot_diagnostics(exp1,show='errors',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor='#1f77b4')
-        plot_diagnostics(exp2,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='#9EC9E9')
-        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)')
-        #ax1.set_xlim([0,0.2])
-        #ax2.set_xlim([0,0.2])
-        #ax3.set_xlim([0,0.02])
-        #ax4.set_xlim([0,10])
-        fig.savefig('fig03_with_errors.png',format='png',dpi=300)
-        p.close(fig)
+        with open(workdir+'/da_exp_A2.json', 'r') as fp:
+            da_settings_A2 = json.load(fp)
+        exp_A2 = load_or_run(da_settings_A2)
 
-        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
-        fig.set_size_inches(6,4)
-        plot_diagnostics(exp1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor=u'#1f77b4')
-        plot_diagnostics(exp2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor=u'#9EC9E9')
-        ax1.set_title('a) Analysis error')
-        ax1.set_xlabel(r'RMSD$^a_\theta$')
-        ax2.set_title('b) First-guess error')
-        ax2.set_xlabel(r'RMSD$^b_\theta$')
-        ax3.set_title('c) Error reduction')
-        ax3.set_xlabel(r'RMSD$^b_\theta-$RMSD$^a_\theta$')
-        ax4.set_title('d) Spread-error consistency')
-        ax4.set_xlabel(r'$\sigma^b_\theta$/RMSD$^b_\theta$')
-        ax1.set_ylabel('height (m)')
-        ax3.set_ylabel('height (m)')
-        #ax1.set_xlim([0,0.2])
-        #ax2.set_xlim([0,0.2])
-        #ax3.set_xlim([0,0.02])
-        #ax4.set_xlim([0,10])
-        fig.savefig('fig03_with_deviations.png',format='png',dpi=300)
+        with open(workdir+'/da_exp_B1.json', 'r') as fp:
+            da_settings_B1 = json.load(fp)
+        exp_B1 = load_or_run(da_settings_B1)
+
+        with open(workdir+'/da_exp_B2.json', 'r') as fp:
+            da_settings_B2 = json.load(fp)
+        exp_B2 = load_or_run(da_settings_B2)
+
+        with open(workdir+'/da_exp_B3.json', 'r') as fp:
+            da_settings_B3 = json.load(fp)
+        exp_B3 = load_or_run(da_settings_B3)
+
+        with open(workdir+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp_B4 = load_or_run(da_settings_B4)
+
+        with open(workdir+'/da_exp_B4_noPE.json', 'r') as fp:
+            da_settings_B4_noPE = json.load(fp)
+        exp_B4_noPE = load_or_run(da_settings_B4_noPE)
+
+        with open(workdir+'/da_exp_C1.json', 'r') as fp:
+            da_settings_C1 = json.load(fp)
+        exp_C1 = load_or_run(da_settings_C1)
+
+        with open(workdir+'/da_exp_C2.json', 'r') as fp:
+            da_settings_C2 = json.load(fp)
+        exp_C2 = load_or_run(da_settings_C2)
+
+        with open(workdir+'/da_exp_D1.json', 'r') as fp:
+            da_settings_D1 = json.load(fp)
+        exp_D1 = load_or_run(da_settings_D1)
+
+        with open(workdir+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D2 = load_or_run(da_settings_D2)
+
+        # Plots: Scatter plots of Desroziers statistics
+        #experiments1 = [exp_A1, exp_B1, exp_B2, exp_B3, exp_B4, exp_C1, exp_C2, exp_D1, exp_D2 ]
+        #labels1 = [r'A$_1$', r'B$_1$', r'B$_2$', r'B$_3$', r'B$_4$', r'C$_1$', r'C$_2$', r'D$_1$', r'D$_2$' ]
+        #experiments2 = [exp_B2, exp_B3, exp_D1, exp_D2 ]
+        #labels2      = [r'B$_2$', r'B$_3$', r'D$_1$', r'D$_2$' ]
+        #experiments1 = [exp_A1, exp_A2, exp_B1, exp_B4, exp_C1, exp_C2 ]
+        #labels1      = [r'A$_1$', r'A$_2$',r'B$_1$', r'B$_4$', r'C$_1$', r'C$_2$' ]
+        experiments1 = [exp_A1, exp_B1, exp_B4, exp_C1, exp_C2 ]
+        labels1      = [r'A$_1$', r'B$_1$', r'B$_4$', r'C$_1$', r'C$_2$' ]
+        experiments2 = [exp_B2, exp_B3, exp_D2 ]
+        labels2      = [r'B$_2$', r'B$_3$', r'D' ]
+
+        fig = p.figure(0)
+        fig.set_size_inches(6,3)
+        ax1 = fig.add_subplot(1,2,1)
+        ax2 = fig.add_subplot(1,2,2)
+
+        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_title(r"a) Exp. with $\sigma_{y_o}^\mathsf{assim}=0.1$ K")
+        ax2.set_title(r"b) Exp. with $\sigma_{y_o}^\mathsf{assim}=1$ K")
+
+        fig.tight_layout()
+        fig.savefig('fig03.png',format='png',dpi=150)
         p.close(fig)
 
     if fig04:
 
-        exp_B1 = pickle.load(open(workdir+"exp_B1.pickle", "rb"))
-        exp_B2 = pickle.load(open(workdir+"exp_B2.pickle", "rb"))
-        exp_B3 = pickle.load(open(workdir+"exp_B3.pickle", "rb"))
-        exp_B4 = pickle.load(open(workdir+"exp_B4.pickle", "rb"))
+        with open(workdir+'/da_exp_B1.json', 'r') as fp:
+            da_settings_B1 = json.load(fp)
+        exp_B1 = load_or_run(da_settings_B1)
+
+        with open(workdir+'/da_exp_B2.json', 'r') as fp:
+            da_settings_B2 = json.load(fp)
+        exp_B2 = load_or_run(da_settings_B2)
+
+        with open(workdir+'/da_exp_B3.json', 'r') as fp:
+            da_settings_B3 = json.load(fp)
+        exp_B3 = load_or_run(da_settings_B3)
+
+        with open(workdir+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp_B4 = load_or_run(da_settings_B4)
 
         fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
         fig.set_size_inches(6,4)
@@ -300,8 +300,13 @@ if __name__ == '__main__':
 
     if fig05:
 
-        exp_C1 = pickle.load(open(workdir+"exp_C1.pickle", "rb"))
-        exp_C2 = pickle.load(open(workdir+"exp_C2.pickle", "rb"))
+        with open(workdir+'/da_exp_C1.json', 'r') as fp:
+            da_settings_C1 = json.load(fp)
+        exp_C1 = load_or_run(da_settings_C1)
+
+        with open(workdir+'/da_exp_C2.json', 'r') as fp:
+            da_settings_C2 = json.load(fp)
+        exp_C2 = load_or_run(da_settings_C2)
 
         fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
         fig.set_size_inches(6,2)
@@ -324,8 +329,13 @@ if __name__ == '__main__':
 
     if fig06:
 
-        exp_C1 = pickle.load(open(workdir+"exp_C1.pickle", "rb"))
-        exp_C2 = pickle.load(open(workdir+"exp_C2.pickle", "rb"))
+        with open(workdir+'/da_exp_C1.json', 'r') as fp:
+            da_settings_C1 = json.load(fp)
+        exp_C1 = load_or_run(da_settings_C1)
+
+        with open(workdir+'/da_exp_C2.json', 'r') as fp:
+            da_settings_C2 = json.load(fp)
+        exp_C2 = load_or_run(da_settings_C2)
 
         npar2 = exp_C1.nr.parameter_number
         npar3 = exp_C2.nr.parameter_number
@@ -360,7 +370,9 @@ if __name__ == '__main__':
 
     if fig07:
 
-        exp_D = pickle.load(open(workdir+"exp_D2.pickle", "rb"))
+        with open(workdir+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D = load_or_run(da_settings_D2)
 
         fig, [ax0, ax1] = p.subplots(1,2,constrained_layout=True)
         fig.set_size_inches(6,3)
@@ -385,38 +397,25 @@ if __name__ == '__main__':
 
     if fig08:
 
-        exp_A = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
-        exp_D = pickle.load(open(workdir+"exp_D2.pickle", "rb"))
+        with open(workdir+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp_B4 = load_or_run(da_settings_B4)
 
-        fig, [[ax0, ax1],[ax2,ax3]] = p.subplots(2,2,constrained_layout=True)
-        fig.set_size_inches(6,6)
+        fig, [ax3,ax2] = p.subplots(1,2,constrained_layout=True)
+        fig.set_size_inches(6,3)
         #
-        [ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A1['obs_error_sdev_assimilate'],ax=[ax1,ax0])
-        ax0.set_title(r'a) Exp. A$_1$, $\delta\overline{y}$ (K)')
-        ax0.set_xlabel('Time (h)')
-        ax0.set_ylabel('Height (m)')
-        ax1.set_title(r'b) Exp. A$_1$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
-        ax1.set_xlabel('Time (h)')
-        ax1.set_ylabel('Height (m)')
-        maxtime = exp_A.trun
-        ax0.set_xlim([0,maxtime/3600])
-        ax1.set_xlim([0,maxtime/3600])
-        ax0.set_xticks(np.arange((maxtime+1)/3600))
-        ax1.set_xticks(np.arange((maxtime+1)/3600))
-        [ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D2['obs_error_sdev_assimilate'],ax=[ax3,ax2])
-        ax2.set_title(r'c) Exp. D$_2$, $\delta\overline{y}$ (K)')
+        [ax2,ax3],c2,c3 = plot_CBL_identifiability(exp_B4,da_settings_B4['obs_error_sdev_assimilate'],ax=[ax3,ax2])
+        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'd) Exp. D$_2$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
+        ax3.set_title(r'b) Exp. B$_4$, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
         ax3.set_xlabel('Time (h)')
         ax3.set_ylabel('Height (m)')
-        maxtime = exp_D.trun
+        maxtime = exp_B4.trun
         ax2.set_xlim([0,maxtime/3600])
         ax3.set_xlim([0,maxtime/3600])
         ax2.set_xticks(np.arange((maxtime+1)/3600))
         ax3.set_xticks(np.arange((maxtime+1)/3600))
-        p.colorbar(c0,orientation='horizontal')
-        p.colorbar(c1,orientation='horizontal')
         p.colorbar(c2,orientation='horizontal')
         p.colorbar(c3,orientation='horizontal')
         #
@@ -425,9 +424,74 @@ if __name__ == '__main__':
 
     if fig09:
 
-        path_EAKF = ""
-        path_LETKF = "./runs/seed=181612_enssize=20_LETKF_6hrs/"
+        workdir_eakf = "./runs/seed=181612_enssize=200_EAKF_6hrs/"
+        workdir_letkf = "./runs/seed=181612_enssize=200_LETKF_6hrs/"
+
+        with open(workdir_eakf+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A1a = load_or_run(da_settings_A1)
+
+        with open(workdir_eakf+'/da_exp_B1.json', 'r') as fp:
+            da_settings_B1 = json.load(fp)
+        exp_B1a = load_or_run(da_settings_B1)
+
+        with open(workdir_eakf+'/da_exp_B2.json', 'r') as fp:
+            da_settings_B2 = json.load(fp)
+        exp_B2a = load_or_run(da_settings_B2)
+
+        with open(workdir_eakf+'/da_exp_B3.json', 'r') as fp:
+            da_settings_B3 = json.load(fp)
+        exp_B3a = load_or_run(da_settings_B3)
 
+        with open(workdir_eakf+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp_B4a = load_or_run(da_settings_B4)
+
+        with open(workdir_eakf+'/da_exp_C1.json', 'r') as fp:
+            da_settings_C1 = json.load(fp)
+        exp_C1a = load_or_run(da_settings_C1)
+
+        with open(workdir_eakf+'/da_exp_C2.json', 'r') as fp:
+            da_settings_C2 = json.load(fp)
+        exp_C2a = load_or_run(da_settings_C2)
+
+        with open(workdir_eakf+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D2a = load_or_run(da_settings_D2)
+
+        with open(workdir_letkf+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A1b = load_or_run(da_settings_A1)
+
+        with open(workdir_letkf+'/da_exp_B1.json', 'r') as fp:
+            da_settings_B1 = json.load(fp)
+        exp_B1b = load_or_run(da_settings_B1)
+
+        with open(workdir_letkf+'/da_exp_B2.json', 'r') as fp:
+            da_settings_B2 = json.load(fp)
+        exp_B2b = load_or_run(da_settings_B2)
+
+        with open(workdir_letkf+'/da_exp_B3.json', 'r') as fp:
+            da_settings_B3 = json.load(fp)
+        exp_B3b = load_or_run(da_settings_B3)
+
+        with open(workdir_letkf+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp_B4b = load_or_run(da_settings_B4)
+
+        with open(workdir_letkf+'/da_exp_C1.json', 'r') as fp:
+            da_settings_C1 = json.load(fp)
+        exp_C1b = load_or_run(da_settings_C1)
+
+        with open(workdir_letkf+'/da_exp_C2.json', 'r') as fp:
+            da_settings_C2 = json.load(fp)
+        exp_C2b = load_or_run(da_settings_C2)
+
+        with open(workdir_letkf+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D2b = load_or_run(da_settings_D2)
+
+        """
         exp_Aa = pickle.load(open(path_EAKF+"exp_A1.pickle", "rb"))
         exp_B1a = pickle.load(open(path_EAKF+"exp_B1.pickle", "rb"))
         exp_B2a = pickle.load(open(path_EAKF+"exp_B2.pickle", "rb"))
@@ -435,7 +499,6 @@ if __name__ == '__main__':
         exp_B4a = pickle.load(open(path_EAKF+"exp_B4.pickle", "rb"))
         exp_C1a = pickle.load(open(path_EAKF+"exp_C1.pickle", "rb"))
         exp_C2a = pickle.load(open(path_EAKF+"exp_C2.pickle", "rb"))
-        exp_D1a = pickle.load(open(path_EAKF+"exp_D1.pickle", "rb"))
         exp_D2a = pickle.load(open(path_EAKF+"exp_D2.pickle", "rb"))
 
         exp_Ab = pickle.load(open(path_LETKF+"exp_A1.pickle", "rb"))
@@ -445,16 +508,16 @@ if __name__ == '__main__':
         exp_B4b = pickle.load(open(path_LETKF+"exp_B4.pickle", "rb"))
         exp_C1b = pickle.load(open(path_LETKF+"exp_C1.pickle", "rb"))
         exp_C2b = pickle.load(open(path_LETKF+"exp_C2.pickle", "rb"))
-        exp_D1b = pickle.load(open(path_LETKF+"exp_D1.pickle", "rb"))
         exp_D2b = pickle.load(open(path_LETKF+"exp_D2.pickle", "rb"))
+        """
 
         fig, [ax1, ax2] = p.subplots(1,2)
         fig.set_size_inches(6,3)
         #
         linecolors = ['darkviolet','chocolate','darkorange','orange','gold','darkgreen','seagreen','royalblue','black']
-        labels = [r'Exp. A',r'Exp. B$_1$',r'Exp. B$_2$',r'Exp. B$_3$',r'Exp. B$_4$',r'Exp. C$_1$',r'Exp. C$_2$',r'Exp. D$_1$',r'Exp. D$_2$']
-        experiments1 = [exp_Aa,exp_B1a,exp_B2a,exp_B3a,exp_B4a,exp_C1a,exp_C2a,exp_D1a,exp_D2a]
-        experiments2 = [exp_Ab,exp_B1b,exp_B2b,exp_B3b,exp_B4b,exp_C1b,exp_C2b,exp_D1b,exp_D2b]
+        labels = [r'Exp. A$_1$',r'Exp. B$_1$',r'Exp. B$_2$',r'Exp. B$_3$',r'Exp. B$_4$',r'Exp. C$_1$',r'Exp. C$_2$',r'Exp. D']
+        experiments1 = [exp_A1a,exp_B1a,exp_B2a,exp_B3a,exp_B4a,exp_C1a,exp_C2a,exp_D2a]
+        experiments2 = [exp_A1b,exp_B1b,exp_B2b,exp_B3b,exp_B4b,exp_C1b,exp_C2b,exp_D2b]
         #
         for exp,label,linecolor in zip(experiments1,labels,linecolors):
             # Experiment settings
@@ -512,6 +575,119 @@ if __name__ == '__main__':
         fig.savefig('fig09.png',format='png',dpi=300)
         p.close(fig)
 
+    if fig10:
+
+        # For the answer to reviewer A
+        """
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp1 = load_or_run(da_settings_A1)
+        label = "A$_1$"
+
+        with open(workdir+'/da_exp_A2.json', 'r') as fp:
+            da_settings_A2 = json.load(fp)
+        exp2 = load_or_run(da_settings_A2)
+        """
+
+        # For the revised manuscript
+        with open(workdir+'/da_exp_B4.json', 'r') as fp:
+            da_settings_B4 = json.load(fp)
+        exp1 = load_or_run(da_settings_B4)
+        label = "B$_4$"
+
+        with open(workdir+'/da_exp_B4_noPE.json', 'r') as fp:
+            da_settings_B4_noPE = json.load(fp)
+        exp2 = load_or_run(da_settings_B4_noPE)
+        
+        # Diagnostics wrt truths
+        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+        fig.set_size_inches(6,4)
+        plot_diagnostics(exp1,show='errors',axes=[ax1,ax2,ax3,ax4],label='B$_4$ with PE',linecolor='#1f77b4',zmax=1500)
+        plot_diagnostics(exp2,show='errors',axes=[ax1,ax2,ax3,ax4],label='B$_4$ without PE',linecolor='#9EC9E9',zmax=1500)
+        ax1.set_title('a) Analysis error')
+        ax1.set_xlabel(r'RMSE$_{y_a}$')
+        ax2.set_title('b) First-guess error')
+        ax2.set_xlabel(r'RMSE$_{y_b}$')
+        ax3.set_title('c) Error reduction')
+        ax3.set_xlabel(r'RMSE$_{y_b}$-RMSE$_{y_a}$')
+        ax4.set_title('d) Consistency of innovations')
+        ax4.set_xlabel(r'$(\sigma_{y_b}+\sigma_{y_o})$/RMSE$_{y_b}$')
+        ax1.set_ylabel('height (m)')
+        ax3.set_ylabel('height (m)')
+        ax1.legend(frameon=False,fontsize=8,loc=7)
+        #ax1.set_xlim([0,0.2])
+        #ax2.set_xlim([0,0.2])
+        #ax3.set_xlim([0,0.02])
+        #ax4.set_xlim([0,10])
+        fig.savefig('fig10.png',format='png',dpi=300)
+        p.close(fig)
+
+        # Diagnostics wrt truths
+        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+        fig.set_size_inches(6,4)
+        plot_diagnostics(exp1,show='errors',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor='#1f77b4')
+        plot_diagnostics(exp2,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='#9EC9E9')
+        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)')
+        #ax1.set_xlim([0,0.2])
+        #ax2.set_xlim([0,0.2])
+        #ax3.set_xlim([0,0.02])
+        #ax4.set_xlim([0,10])
+        fig.savefig('diagnostics_profiles_with_errors.png',format='png',dpi=300)
+        p.close(fig)
+
+        # Diagnostics wrt observations
+        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+        fig.set_size_inches(6,4)
+        plot_diagnostics(exp1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor=u'#1f77b4')
+        plot_diagnostics(exp2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor=u'#9EC9E9')
+        ax1.set_title('a) Analysis error')
+        ax1.set_xlabel(r'RMSD$^a_\theta$')
+        ax2.set_title('b) First-guess error')
+        ax2.set_xlabel(r'RMSD$^b_\theta$')
+        ax3.set_title('c) Error reduction')
+        ax3.set_xlabel(r'RMSD$^b_\theta-$RMSD$^a_\theta$')
+        ax4.set_title('d) Spread-error consistency')
+        ax4.set_xlabel(r'$\sigma^b_\theta$/RMSD$^b_\theta$')
+        ax1.set_ylabel('height (m)')
+        ax3.set_ylabel('height (m)')
+        #ax1.set_xlim([0,0.2])
+        #ax2.set_xlim([0,0.2])
+        #ax3.set_xlim([0,0.02])
+        #ax4.set_xlim([0,10])
+        fig.savefig('diagnostics_profiles_with_deviations.png',format='png',dpi=300)
+        p.close(fig)
+
+        # Diagnostics wrt truths and observations compared
+        if exp1.da.obs_coordinates.ndim == 1:
+            z = np.sort(exp1.da.obs_coordinates)
+        elif exp1.da.obs_coordinates.ndim == 2:
+            z = np.sort(exp1.da.obs_coordinates[0,:])
+
+        sigma_o = exp1.dg.varo**0.5
+        brmsd = exp1.dg.bRMSD_t
+        bsprd = exp1.dg.bSprd_t
+        brmse = exp1.dg.bRMSE_t
+
+        fig, ax1 = p.subplots(1,1,constrained_layout=True)
+        fig.set_size_inches(4,3)
+        ax1.plot(np.sqrt(bsprd**2+sigma_o**2)/brmse,z,color=u'#1f77b4',label=r'$\sigma^b_\theta$/RMSE$^b_\theta$; reference is truth')
+        ax1.plot(np.sqrt(bsprd**2+sigma_o**2)/brmsd,z,color=u'#9EC9E9',label=r'$\sigma^b_\theta$/RMSD$^b_\theta$; reference is observations')
+        ax1.set_title('Exp. %s'%label)
+        ax1.set_xlabel(r'spread/error ratio')
+        ax1.set_ylabel('height (m)')
+        ax1.legend(frameon=False,fontsize=8,loc=7)
+        fig.savefig('diagnostics_profiles_errors_vs_deviations.png',format='png',dpi=300)
+        p.close(fig)
+
     if opt01:
 
         with open('./default_da.json', 'r') as fp:
@@ -522,6 +698,10 @@ if __name__ == '__main__':
 
     if opt02:
 
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp1 = load_or_run(da_settings_A1)
+
         plot_profiles(exp_A1,da_settings_A1,'profiles.png')
 
     if opt03:
@@ -541,8 +721,8 @@ if __name__ == '__main__':
         for exp, label in zip([exp_RC, exp_R],['RC','R']):
             plot_overview(exp,label,'exp_%s.png'%label)
 
-        plot_profiles(exp_RC,da_settings_RC,'expRC_profiles.png')
-        plot_profiles(exp_R,da_settings_R,'expR_profiles.png')
+        plot_profiles(exp_RC,da_settings_RC,'exp_RC_profiles.png')
+        plot_profiles(exp_R,da_settings_R,'exp_R_profiles.png')
 
     if opt04:
 
@@ -570,20 +750,20 @@ if __name__ == '__main__':
 
     if opt05:
 
-        with open('./config_files/da_exp_L1.json', 'r') as fp:
+        with open(workdir+'/da_exp_L1.json', 'r') as fp:
             da_settings_localization_1 = json.load(fp)
         exp_localization_1 = load_or_run(da_settings_localization_1)
 
-        with open('./config_files/da_exp_A1.json', 'r') as fp:
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
             da_settings_localization_2 = json.load(fp)
         exp_localization_2 = load_or_run(da_settings_localization_2)
         exp_localization_2.label = 'L2'
 
-        with open('./config_files/da_exp_L3.json', 'r') as fp:
+        with open(workdir+'/da_exp_L3.json', 'r') as fp:
             da_settings_localization_3 = json.load(fp)
         exp_localization_3 = load_or_run(da_settings_localization_3)
 
-        with open('./config_files/da_exp_L4.json', 'r') as fp:
+        with open(workdir+'/da_exp_L4.json', 'r') as fp:
             da_settings_localization_4 = json.load(fp)
         exp_localization_4 = load_or_run(da_settings_localization_4)
 
@@ -596,20 +776,20 @@ if __name__ == '__main__':
 
     if opt06:
 
-        with open('./config_files/da_exp_I1.json', 'r') as fp:
+        with open(workdir+'/da_exp_I1.json', 'r') as fp:
             da_settings_inflation_1 = json.load(fp)
         exp_inflation_1 = load_or_run(da_settings_inflation_1)
 
-        with open('./config_files/da_exp_I2.json', 'r') as fp:
+        with open(workdir+'/da_exp_I2.json', 'r') as fp:
             da_settings_inflation_2 = json.load(fp)
         exp_inflation_2 = load_or_run(da_settings_inflation_2)
 
-        with open('./config_files/da_exp_A1.json', 'r') as fp:
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
             da_settings_inflation_3 = json.load(fp)
         exp_inflation_3 = load_or_run(da_settings_inflation_3)
         exp_inflation_3.label = 'I3'
 
-        with open('./config_files/da_exp_I4.json', 'r') as fp:
+        with open(workdir+'/da_exp_I4.json', 'r') as fp:
             da_settings_inflation_4 = json.load(fp)
         exp_inflation_4 = load_or_run(da_settings_inflation_4)
 
@@ -622,20 +802,20 @@ if __name__ == '__main__':
 
     if opt07:
 
-        with open('./config_files/da_exp_I5.json', 'r') as fp:
+        with open(workdir+'/da_exp_I5.json', 'r') as fp:
             da_settings_inflation_5 = json.load(fp)
         exp_inflation_5 = load_or_run(da_settings_inflation_5)
 
-        with open('./config_files/da_exp_A1.json', 'r') as fp:
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
             da_settings_inflation_6 = json.load(fp)
         exp_inflation_6 = load_or_run(da_settings_inflation_6)
         exp_inflation_6.label = 'I6'
 
-        with open('./config_files/da_exp_I7.json', 'r') as fp:
+        with open(workdir+'/da_exp_I7.json', 'r') as fp:
             da_settings_inflation_7 = json.load(fp)
         exp_inflation_7 = load_or_run(da_settings_inflation_7)
 
-        with open('./config_files/da_exp_I8.json', 'r') as fp:
+        with open(workdir+'/da_exp_I8.json', 'r') as fp:
             da_settings_inflation_8 = json.load(fp)
         exp_inflation_8 = load_or_run(da_settings_inflation_8)
 
@@ -648,12 +828,12 @@ if __name__ == '__main__':
 
     if opt08:
 
-        with open('./config_files/da_exp_A1.json', 'r') as fp:
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
             da_settings_O1 = json.load(fp)
         exp_O1 = load_or_run(da_settings_O1)
         exp_O1.label = 'O1'
 
-        with open('./config_files/da_exp_O2.json', 'r') as fp:
+        with open(workdir+'/da_exp_O2.json', 'r') as fp:
             da_settings_O2 = json.load(fp)
         exp_O2 = load_or_run(da_settings_O2)
 
@@ -661,7 +841,7 @@ if __name__ == '__main__':
         for exp,settings,label in zip([exp_O1,exp_O2],\
                                 [da_settings_O1,da_settings_O2],\
                                 [r"O$_1$",r"O$_2$"]):
-            print(exp.randomize_obs)
+            print('Random obs order? %s'%exp.randomize_obs)
 
             plot_overview(exp,label,'exp_%s.png'%exp.label)
 
@@ -697,7 +877,7 @@ if __name__ == '__main__':
 
         fig, ax = p.subplots(1,1,constrained_layout=True)
         fig.set_size_inches(4,3)
-        c0 = ax.pcolormesh(times/3600.,zobs1,delta_theta,cmap='RdBu_r')
+        c0 = ax.pcolormesh(times/3600.,zobs1,delta_theta,cmap='RdBu_r',norm=colors.CenteredNorm())
         ax.set_xlim([0,maxtime/3600])
         ax.set_xticks(np.arange((maxtime+1)/3600))
         ax.set_xlabel('Time (h)')
@@ -709,31 +889,25 @@ if __name__ == '__main__':
 
     if opt09:
 
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A1 = load_or_run(da_settings_A1)
+
         plot_B_matrix(exp_A1,timeid=71)
 
     if opt10:
 
-        # Plots: Scatter plots of Desroziers statistics
-        #experiments = [exp_A1, exp_B1, exp_B2, exp_B3, exp_B4, exp_C1, exp_C2, exp_D1, exp_D2 ]
-        #labels = [r'A$_1$', r'B$_1$', r'B$_2$', r'B$_3$', r'B$_4$', r'C$_1$', r'C$_2$', r'D$_1$', r'D$_2$' ]
-        experiments1 = [exp_A1, exp_B1, exp_B4, exp_C1, exp_C2 ]
-        labels1      = [r'A$_1$', r'B$_1$', r'B$_4$', r'C$_1$', r'C$_2$' ]
-        experiments2 = [exp_B2, exp_B3, exp_D1, exp_D2 ]
-        labels2      = [r'B$_2$', r'B$_3$', r'D$_1$', r'D$_2$' ]
+        with open(workdir+'/da_exp_B3.json', 'r') as fp:
+            da_settings_B3 = json.load(fp)
+        exp_B3 = load_or_run(da_settings_B3)
 
-        fig = p.figure(0)
-        fig.set_size_inches(6,3)
-        ax1 = fig.add_subplot(1,2,1)
-        ax2 = fig.add_subplot(1,2,2)
+        with open(workdir+'/da_exp_D1.json', 'r') as fp:
+            da_settings_D1 = json.load(fp)
+        exp_D1 = load_or_run(da_settings_D1)
 
-        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)
-
-        fig.tight_layout()
-        fig.savefig('Diag_desroziers.png',format='png',dpi=150)
-        p.close(fig)
-
-    if opt11:
+        with open(workdir+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D2 = load_or_run(da_settings_D2)
 
         # Plots
         plot_overview(exp_D1,r"D$_1$",'exp_D1.png')
@@ -742,9 +916,9 @@ if __name__ == '__main__':
 
         fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
         fig.set_size_inches(6,4)
-        plot_diagnostics(exp_D1,show='errors',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor='#1f77b4')
-        plot_diagnostics(exp_D2,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='#9EC9E9')
-        plot_diagnostics(exp_B3,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='orange')
+        plot_diagnostics(exp_D1,show='errors',axes=[ax1,ax2,ax3,ax4],label='D$_1$',linecolor='#1f77b4')
+        plot_diagnostics(exp_D2,show='errors',axes=[ax1,ax2,ax3,ax4],label='D$_2$',linecolor='#9EC9E9')
+        plot_diagnostics(exp_B3,show='errors',axes=[ax1,ax2,ax3,ax4],label='B$_3$',linecolor='orange')
         ax1.set_title('a) Analysis error')
         ax1.set_xlabel(r'RMSE$^a_\theta$')
         ax2.set_title('b) First-guess error')
@@ -755,6 +929,7 @@ if __name__ == '__main__':
         ax4.set_xlabel(r'$\sigma^b_\theta$/RMSE$^b_\theta$')
         ax1.set_ylabel('height (m)')
         ax3.set_ylabel('height (m)')
+        ax1.legend(frameon=False,loc=1)
         #ax1.set_xlim([0,0.2])
         #ax2.set_xlim([0,0.2])
         #ax3.set_xlim([0,0.02])
@@ -764,9 +939,9 @@ if __name__ == '__main__':
 
         fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
         fig.set_size_inches(6,4)
-        plot_diagnostics(exp_D1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor=u'#1f77b4')
-        plot_diagnostics(exp_D2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor=u'#9EC9E9')
-        plot_diagnostics(exp_B3,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='orange')
+        plot_diagnostics(exp_D1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='D$_1$',linecolor=u'#1f77b4')
+        plot_diagnostics(exp_D2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='D$_2$',linecolor=u'#9EC9E9')
+        plot_diagnostics(exp_B3,show='deviations',axes=[ax1,ax2,ax3,ax4],label='B$_3$',linecolor='orange')
         ax1.set_title('a) Analysis error')
         ax1.set_xlabel(r'RMSD$^a_\theta$')
         ax2.set_title('b) First-guess error')
@@ -777,6 +952,7 @@ if __name__ == '__main__':
         ax4.set_xlabel(r'$\sigma^b_\theta$/RMSD$^b_\theta$')
         ax1.set_ylabel('height (m)')
         ax3.set_ylabel('height (m)')
+        ax1.legend(frameon=False,loc=1)
         #ax1.set_xlim([0,0.2])
         #ax2.set_xlim([0,0.2])
         #ax3.set_xlim([0,0.02])
@@ -784,6 +960,60 @@ if __name__ == '__main__':
         fig.savefig('profiles_deviations.png',format='png',dpi=300)
         p.close(fig)
 
+    if opt11:
+
+        with open(workdir+'/da_exp_I9.json', 'r') as fp:
+            da_settings_I9 = json.load(fp)
+        exp_I9 = load_or_run(da_settings_I9)
+
+        # Plots
+        plot_overview(exp_I9,r"I$_9$",'exp_I9.png')
+
+    if opt12:
+
+        with open(workdir+'/da_exp_A1.json', 'r') as fp:
+            da_settings_A1 = json.load(fp)
+        exp_A = load_or_run(da_settings_A1)
+
+        with open(workdir+'/da_exp_D2.json', 'r') as fp:
+            da_settings_D2 = json.load(fp)
+        exp_D = load_or_run(da_settings_D2)
+
+        fig, [[ax0, ax1],[ax2,ax3]] = p.subplots(2,2,constrained_layout=True)
+        fig.set_size_inches(6,6)
+        #
+        [ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A1['obs_error_sdev_assimilate'],ax=[ax1,ax0])
+        ax0.set_title(r'a) Exp. A$_1$, $\delta\overline{y}$ (K)')
+        ax0.set_xlabel('Time (h)')
+        ax0.set_ylabel('Height (m)')
+        ax1.set_title(r'b) Exp. A$_1$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
+        ax1.set_xlabel('Time (h)')
+        ax1.set_ylabel('Height (m)')
+        maxtime = exp_A.trun
+        ax0.set_xlim([0,maxtime/3600])
+        ax1.set_xlim([0,maxtime/3600])
+        ax0.set_xticks(np.arange((maxtime+1)/3600))
+        ax1.set_xticks(np.arange((maxtime+1)/3600))
+        [ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D2['obs_error_sdev_assimilate'],ax=[ax3,ax2])
+        ax2.set_title(r'c) Exp. D$_2$, $\delta\overline{y}$ (K)')
+        ax2.set_xlabel('Time (h)')
+        ax2.set_ylabel('Height (m)')
+        ax3.set_title(r'd) Exp. D$_2$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
+        ax3.set_xlabel('Time (h)')
+        ax3.set_ylabel('Height (m)')
+        maxtime = exp_D.trun
+        ax2.set_xlim([0,maxtime/3600])
+        ax3.set_xlim([0,maxtime/3600])
+        ax2.set_xticks(np.arange((maxtime+1)/3600))
+        ax3.set_xticks(np.arange((maxtime+1)/3600))
+        p.colorbar(c0,orientation='horizontal')
+        p.colorbar(c1,orientation='horizontal')
+        p.colorbar(c2,orientation='horizontal')
+        p.colorbar(c3,orientation='horizontal')
+        #
+        fig.savefig('fig08_old.png',format='png',dpi=300)
+        p.close(fig)
+    
     if check_reproducibility:
 
         for i in ["exp_A1.pickle",