diff --git a/PE_CBL.py b/PE_CBL.py
index 20daccf1c0418fd455ca9295594d42de5b1f68cc..faff346e82210537aa34a0b1377ddef9a663a0cb 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -108,6 +108,9 @@ if __name__ == '__main__':
     opt01 = False # assimilation of a single observation
     opt02 = False # assimilation of profiles at two times
 
+    # Whether or not to run experiments without parameter estimation
+    noPE_runs = False
+
     # Default PE experiment
     # Create a copy of the default settings
     cbl_settings_A = dict(default_cbl_settings)
@@ -127,6 +130,20 @@ if __name__ == '__main__':
     truths = exp_A.truths
     observations = exp_A.observations
 
+    # Corresponding experiment without parameter estimation
+    cbl_settings_A_noPE = dict(cbl_settings_A)
+    da_settings_A_noPE = dict(da_settings_A)
+    da_settings_A_noPE['nr'] = deepcopy(nature_run)
+    da_settings_A_noPE['truths'] = truths
+    da_settings_A_noPE['observations'] = observations
+    da_settings_A_noPE['nr'].do_parameter_estimation = False
+    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')) 
+
     if fig01:
 
         # Create a copy of the default settings
@@ -201,20 +218,6 @@ if __name__ == '__main__':
 
     if fig02:
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_A_noPE = dict(cbl_settings_A)
-        da_settings_A_noPE = dict(da_settings_A)
-        da_settings_A_noPE['nr'] = deepcopy(nature_run)
-        da_settings_A_noPE['truths'] = truths
-        da_settings_A_noPE['observations'] = observations
-        da_settings_A_noPE['nr'].do_parameter_estimation = False
-        cbl_settings_A_noPE['do_parameter_estimation'] = False
-        da_settings_A_noPE['cbl_settings'] = cbl_settings_A_noPE
-        exp_A_noPE = experiment(da_settings_A_noPE)
-
-        # Save experiments to disk
-        pickle.dump(exp_A_noPE, open('exp_A_noPE.pickle', 'wb')) 
-
         # Make plots
         fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
         fig.set_size_inches(6,6)
@@ -242,6 +245,59 @@ 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]
+        labels = ["Exp. A"]
+        linecolors = p.rcParams['axes.prop_cycle'].by_key()['color']
+
+        plotfig(range(1), 'fig03.png')
+
+    if fig04:
+
         ########################################################################
         # Create a copy of the default settings
         # Then re-use the available nature run and derived information
@@ -262,17 +318,18 @@ if __name__ == '__main__':
         setattr(exp_B1,'label','B1')
         pickle.dump(exp_B1, open('exp_B1.pickle', 'wb')) 
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_B1_noPE = dict(cbl_settings_B1)
-        da_settings_B1_noPE = dict(da_settings_B1)
-        cbl_settings_B1_noPE['do_parameter_estimation'] = False
-        da_settings_B1_noPE['nr'] = deepcopy(nature_run)
-        da_settings_B1_noPE['nr'].do_parameter_estimation = False
-        da_settings_B1_noPE['cbl_settings'] = cbl_settings_B1_noPE
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_B1_noPE = dict(cbl_settings_B1)
+            da_settings_B1_noPE = dict(da_settings_B1)
+            cbl_settings_B1_noPE['do_parameter_estimation'] = False
+            da_settings_B1_noPE['nr'] = deepcopy(nature_run)
+            da_settings_B1_noPE['nr'].do_parameter_estimation = False
+            da_settings_B1_noPE['cbl_settings'] = cbl_settings_B1_noPE
 
-        # Run and save to disk
-        exp_B1_noPE = experiment(da_settings_B1_noPE)
-        pickle.dump(exp_B1_noPE, open('exp_B1_noPE.pickle', 'wb'))
+            # Run and save to disk
+            exp_B1_noPE = experiment(da_settings_B1_noPE)
+            pickle.dump(exp_B1_noPE, open('exp_B1_noPE.pickle', 'wb'))
 
         ########################################################################
         # Create a copy of the default settings
@@ -294,17 +351,18 @@ if __name__ == '__main__':
         setattr(exp_B2,'label','B2')
         pickle.dump(exp_B2, open('exp_B2.pickle', 'wb')) 
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_B2_noPE = dict(cbl_settings_B2)
-        da_settings_B2_noPE = dict(da_settings_B2)
-        cbl_settings_B2_noPE['do_parameter_estimation'] = False
-        da_settings_B2_noPE['nr'] = deepcopy(nature_run)
-        da_settings_B2_noPE['nr'].do_parameter_estimation = False
-        da_settings_B2_noPE['cbl_settings'] = cbl_settings_B2_noPE
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_B2_noPE = dict(cbl_settings_B2)
+            da_settings_B2_noPE = dict(da_settings_B2)
+            cbl_settings_B2_noPE['do_parameter_estimation'] = False
+            da_settings_B2_noPE['nr'] = deepcopy(nature_run)
+            da_settings_B2_noPE['nr'].do_parameter_estimation = False
+            da_settings_B2_noPE['cbl_settings'] = cbl_settings_B2_noPE
 
-        # Run and save to disk
-        exp_B2_noPE = experiment(da_settings_B2_noPE)
-        pickle.dump(exp_B2_noPE, open('exp_B2_noPE.pickle', 'wb'))
+            # Run and save to disk
+            exp_B2_noPE = experiment(da_settings_B2_noPE)
+            pickle.dump(exp_B2_noPE, open('exp_B2_noPE.pickle', 'wb'))
 
         ########################################################################
         # Create a copy of the default settings
@@ -326,79 +384,114 @@ if __name__ == '__main__':
         setattr(exp_B3,'label','B3')
         pickle.dump(exp_B3, open('exp_B3.pickle', 'wb')) 
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_B3_noPE = dict(cbl_settings_B3)
-        da_settings_B3_noPE = dict(da_settings_B3)
-        cbl_settings_B3_noPE['do_parameter_estimation'] = False
-        da_settings_B3_noPE['nr'] = deepcopy(nature_run)
-        da_settings_B3_noPE['nr'].do_parameter_estimation = False
-        da_settings_B3_noPE['cbl_settings'] = cbl_settings_B3_noPE
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_B3_noPE = dict(cbl_settings_B3)
+            da_settings_B3_noPE = dict(da_settings_B3)
+            cbl_settings_B3_noPE['do_parameter_estimation'] = False
+            da_settings_B3_noPE['nr'] = deepcopy(nature_run)
+            da_settings_B3_noPE['nr'].do_parameter_estimation = False
+            da_settings_B3_noPE['cbl_settings'] = cbl_settings_B3_noPE
 
-        # Run and save to disk
-        exp_B3_noPE = experiment(da_settings_B3_noPE)
-        pickle.dump(exp_B3_noPE, open('exp_B3_noPE.pickle', 'wb'))
+            # Run and save to disk
+            exp_B3_noPE = experiment(da_settings_B3_noPE)
+            pickle.dump(exp_B3_noPE, open('exp_B3_noPE.pickle', 'wb'))
 
         ########################################################################
         # Create a copy of the default settings
+        cbl_settings_B4 = dict(default_cbl_settings)
+        da_settings_B4 = dict(default_da_settings)
+
+        # Change settings as necessary
+        # Changes include generation of observations, so the existing nature run
+        # can't be reused.
+        cbl_settings_B4['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
+        da_settings_B4['cbl_settings'] = cbl_settings_B4
+        da_settings_B4['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*5
+        da_settings_B4['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
+
+        # Run and save to disk
+        exp_B4 = experiment(da_settings_B4)
+        setattr(exp_B4,'label','B4')
+        pickle.dump(exp_B4, open('exp_B4.pickle', 'wb')) 
+
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_B4_noPE = dict(cbl_settings_B4)
+            da_settings_B4_noPE = dict(da_settings_B4)
+            cbl_settings_B4_noPE['do_parameter_estimation'] = False
+            da_settings_B4_noPE['cbl_settings'] = cbl_settings_B4_noPE
+
+            # Run and save to disk
+            exp_B4_noPE = experiment(da_settings_B4_noPE)
+            pickle.dump(exp_B4_noPE, open('exp_B4_noPE.pickle', 'wb')) 
+
         ########################################################################
 
         # Make plots
-        fig, [ax1, ax2, ax3] = p.subplots(1,3,constrained_layout=True)
-        fig.set_size_inches(6,2)
+        fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+        fig.set_size_inches(6,4)
         #
-        #ax1,c1 = plot_CBL_identifiability(exp_B,da_settings_B['obs_error_sdev_assimilate'][0],None,ax=ax1)
-        #ax1.set_title(r'a) Exp. B, $K_{p,\theta}$ (K$^{-1}$)')
-        #ax1.set_xlabel('Time (h)')
-        #ax1.set_ylabel('Height (m)')
-        #p.colorbar(c1,orientation='horizontal')
         ax1 = plot_CBL_PE(exp_B1,None,ax=ax1)
         ax1.set_title(r'a) Exp. B$_1$')
         ax1.set_xlabel('Time (h)')
+        ax1.set_ylabel(r'$p$')
         #
         ax2 = plot_CBL_PE(exp_B2,None,ax=ax2)
         ax2.set_title(r'b) Exp. B$_2$')
         ax2.set_xlabel('Time (h)')
+        ax2.set_ylabel(r'$p$')
         #
         ax3 = plot_CBL_PE(exp_B3,None,ax=ax3)
         ax3.set_title(r'c) Exp. B$_3$')
         ax3.set_xlabel('Time (h)')
+        ax3.set_ylabel(r'$p$')
+        #
+        ax4 = plot_CBL_PE(exp_B4,None,ax=ax4)
+        ax4.set_title(r'd) Exp. B$_4$')
+        ax4.set_xlabel('Time (h)')
+        ax4.set_ylabel(r'$p$')
         #
         ax1.set_yticks([0,1,2,3,4,5])
         ax2.sharey(ax1)
-        ax3.sharey(ax1)
+        ax3.set_yticks([0,1,2,3,4,5])
+        ax4.sharey(ax3)
         #
-        fig.savefig('fig03.png',format='png',dpi=300)
+        fig.savefig('fig04.png',format='png',dpi=300)
         p.close(fig)
 
-    if fig04:
+    if fig05:
 
         ########################################################################
         # Create a copy of the default settings
         cbl_settings_C1 = dict(default_cbl_settings)
         da_settings_C1 = dict(default_da_settings)
 
+        # Re-use the available nature run and derived information
+        da_settings_C1['nr'] = deepcopy(nature_run)
+        da_settings_C1['truths'] = truths
+        da_settings_C1['observations'] = observations
+
         # Change settings as necessary
-        # Changes include generation of observations, so the existing nature run
-        # can't be reused.
         cbl_settings_C1['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
+        cbl_settings_C1['is_cgrad'] = False
         da_settings_C1['cbl_settings'] = cbl_settings_C1
-        da_settings_C1['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*5
-        da_settings_C1['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
 
         # Run and save to disk
         exp_C1 = experiment(da_settings_C1)
         setattr(exp_C1,'label','C1')
         pickle.dump(exp_C1, open('exp_C1.pickle', 'wb')) 
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_C1_noPE = dict(cbl_settings_C1)
-        da_settings_C1_noPE = dict(da_settings_C1)
-        cbl_settings_C1_noPE['do_parameter_estimation'] = False
-        da_settings_C1_noPE['cbl_settings'] = cbl_settings_C1_noPE
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_C1_noPE = dict(cbl_settings_C1)
+            da_settings_C1_noPE = dict(da_settings_C1)
+            cbl_settings_C1_noPE['do_parameter_estimation'] = False
+            da_settings_C1_noPE['cbl_settings'] = cbl_settings_C1_noPE
 
-        # Run and save to disk
-        exp_C1_noPE = experiment(da_settings_C1_noPE)
-        pickle.dump(exp_C1_noPE, open('exp_C1_noPE.pickle', 'wb')) 
+            # Run and save to disk
+            exp_C1_noPE = experiment(da_settings_C1_noPE)
+            pickle.dump(exp_C1_noPE, open('exp_C1_noPE.pickle', 'wb')) 
 
         ########################################################################
         # Create a copy of the default settings
@@ -412,7 +505,7 @@ if __name__ == '__main__':
 
         # Change settings as necessary
         cbl_settings_C2['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
-        cbl_settings_C2['is_cgrad'] = False
+        cbl_settings_C2['Hmax'] = 0.15
         da_settings_C2['cbl_settings'] = cbl_settings_C2
 
         # Run and save to disk
@@ -420,76 +513,68 @@ if __name__ == '__main__':
         setattr(exp_C2,'label','C2')
         pickle.dump(exp_C2, open('exp_C2.pickle', 'wb')) 
 
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_C2_noPE = dict(cbl_settings_C2)
-        da_settings_C2_noPE = dict(da_settings_C2)
-        cbl_settings_C2_noPE['do_parameter_estimation'] = False
-        da_settings_C2_noPE['cbl_settings'] = cbl_settings_C2_noPE
+        if noPE_runs:
+            # Corresponding experiment without parameter estimation
+            cbl_settings_C2_noPE = dict(cbl_settings_C2)
+            da_settings_C2_noPE = dict(da_settings_C2)
+            cbl_settings_C2_noPE['do_parameter_estimation'] = False
+            da_settings_C2_noPE['cbl_settings'] = cbl_settings_C2_noPE
 
-        # Run and save to disk
-        exp_C2_noPE = experiment(da_settings_C2_noPE)
-        pickle.dump(exp_C2_noPE, open('exp_C2_noPE.pickle', 'wb')) 
-
-        ########################################################################
-        # Create a copy of the default settings
-        cbl_settings_C3 = dict(default_cbl_settings)
-        da_settings_C3 = dict(default_da_settings)
-
-        # Re-use the available nature run and derived information
-        da_settings_C3['nr'] = deepcopy(nature_run)
-        da_settings_C3['truths'] = truths
-        da_settings_C3['observations'] = observations
-
-        # Change settings as necessary
-        cbl_settings_C3['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
-        cbl_settings_C3['Hmax'] = 0.15
-        da_settings_C3['cbl_settings'] = cbl_settings_C3
-
-        # Run and save to disk
-        exp_C3 = experiment(da_settings_C3)
-        setattr(exp_C3,'label','C3')
-        pickle.dump(exp_C3, open('exp_C3.pickle', 'wb')) 
-
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_C3_noPE = dict(cbl_settings_C3)
-        da_settings_C3_noPE = dict(da_settings_C3)
-        cbl_settings_C3_noPE['do_parameter_estimation'] = False
-        da_settings_C3_noPE['cbl_settings'] = cbl_settings_C3_noPE
-
-        # Run and save to disk
-        exp_C3_noPE = experiment(da_settings_C3_noPE)
-        pickle.dump(exp_C3_noPE, open('exp_C3_noPE.pickle', 'wb')) 
+            # Run and save to disk
+            exp_C2_noPE = experiment(da_settings_C2_noPE)
+            pickle.dump(exp_C2_noPE, open('exp_C2_noPE.pickle', 'wb')) 
 
         ########################################################################
         # Make plots
-        fig, [ax1, ax2, ax3] = p.subplots(1,3,constrained_layout=True)
+        fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
         fig.set_size_inches(6,2)
         #
-        #ax1,c1 = plot_CBL_identifiability(exp_C,da_settings_C['obs_error_sdev_assimilate'][0],None,ax=ax1)
-        #ax1.set_title(r'a) Exp. C, $K_{p,\theta}$ (K$^{-1}$)')
-        #ax1.set_xlabel('Time (h)')
-        #ax1.set_ylabel('Height (m)')
-        #p.colorbar(c1,orientation='horizontal')
         ax1 = plot_CBL_PE(exp_C1,None,ax=ax1)
         ax1.set_title(r'a) Exp. C$_1$')
         ax1.set_xlabel('Time (h)')
+        ax1.set_ylabel(r'$p$')
         #
         ax2 = plot_CBL_PE(exp_C2,None,ax=ax2)
         ax2.set_title(r'b) Exp. C$_2$')
         ax2.set_xlabel('Time (h)')
-        #
-        ax3 = plot_CBL_PE(exp_C3,None,ax=ax3)
-        ax3.set_title(r'c) Exp. C$_3$')
-        ax3.set_xlabel('Time (h)')
+        ax2.set_ylabel(r'$p$')
         #
         ax1.set_yticks([0,1,2,3,4,5])
         ax2.sharey(ax1)
-        ax3.sharey(ax1)
         #
-        fig.savefig('fig04.png',format='png',dpi=300)
+        fig.savefig('fig05.png',format='png',dpi=300)
         p.close(fig)
 
-    if fig05:
+    if fig06:
+
+        exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
+        exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
+
+        npar2 = exp_C1.nr.parameter_number
+        npar3 = exp_C2.nr.parameter_number
+        updates_C1 = exp_C1.da.analyses[:,:-npar2,:].mean(axis=2).T - exp_C1.da.backgrounds[:,:-npar2,:].mean(axis=2).T
+        updates_C2 = exp_C2.da.analyses[:,:-npar3,:].mean(axis=2).T - exp_C2.da.backgrounds[:,:-npar3,:].mean(axis=2).T
+
+        fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
+        fig.set_size_inches(6,3)
+        c1 = ax1.pcolormesh(exp_C1.da.time/3600,exp_C1.nr.zt,updates_C1,vmin=-0.03,vmax=0.03,cmap='RdBu_r')
+        c2 = ax2.pcolormesh(exp_C2.da.time/3600,exp_C2.nr.zt,updates_C2,vmin=-0.03,vmax=0.03,cmap='RdBu_r')
+        #
+        ax1.set_title(r'a) Exp. C$_1$, $\theta_a-\theta_b$')
+        ax2.set_title(r'b) Exp. C$_2$, $\theta_a-\theta_b$')
+        ax1.set_xlabel('Time (h)')
+        ax2.set_xlabel('Time (h)')
+        ax1.set_ylabel('Height (m)')
+        ax2.sharey(ax1)
+        p.setp(ax2.get_yticklabels(), visible=False)
+        #
+        p.colorbar(c1,orientation='horizontal')
+        p.colorbar(c2,orientation='horizontal')
+        #
+        fig.savefig('fig06.png',format='png',dpi=300)
+        p.close(fig)
+
+    if fig07:
 
         # Create a copy of the default settings
         cbl_settings_D = dict(default_cbl_settings)
@@ -508,20 +593,21 @@ if __name__ == '__main__':
         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
 
-        # Run
+        # Run and save to disk
         exp_D = experiment(da_settings_D)
         setattr(exp_D,'label','D')
+        pickle.dump(exp_D, open('exp_D.pickle', 'wb')) 
 
-        # Run the 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
-        exp_D_noPE = experiment(da_settings_D_noPE)
+        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
 
-        # Save experiments to disk
-        pickle.dump(exp_D, open('exp_D.pickle', 'wb')) 
-        pickle.dump(exp_D_noPE, open('exp_D_noPE.pickle', 'wb')) 
+            # 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)
@@ -546,94 +632,6 @@ if __name__ == '__main__':
         p.colorbar(c2,orientation='horizontal')
         #
 
-        fig.savefig('fig05.png',format='png',dpi=300)
-        p.close(fig)
-
-    if fig06:
-
-        def plotfig(exprange,filename):
-            fig, [ax1, ax2, ax3] = p.subplots(1,3,constrained_layout=True)
-            fig.set_size_inches(8,3)
-            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-i1.aRMSE_t,z,label=labels[i],color=linecolors[i])
-                ax2.plot(i2.bRMSE_t-i2.aRMSE_t,z,color=linecolors[i],dashes=[3,1],alpha=0.3)
-                #
-                ax3.plot(i1.bSprd_t/i1.bRMSE_t,z_pbl,label=labels[i],color=linecolors[i])
-                ax3.plot(i2.bSprd_t/i2.bRMSE_t,z_pbl,color=linecolors[i],dashes=[3,1],alpha=0.3)
-            ax1.set_xlabel(r'RMSE$^a_\theta$')
-            ax2.set_xlabel(r'RMSE$^b_\theta-$RMSE$^a_\theta$')
-            ax3.set_xlabel(r'$\sigma^b_\theta$/RMSE$^b_\theta$')
-            ax1.set_ylabel('height (m)')
-            #
-            ax2.legend(frameon=False)
-            ax3.axvline(x=1,color='k',linewidth=0.5,dashes=[3,1])
-            ax2.sharey(ax1)
-            ax3.sharey(ax1)
-            p.setp(ax2.get_yticklabels(), visible=False)
-            p.setp(ax3.get_yticklabels(), visible=False)
-            #
-            fig.savefig(filename,format='png',dpi=300)
-            p.close(fig)
-
-        exp_A = pickle.load(open("exp_A.pickle", "rb"))
-        exp_B1 = pickle.load(open("exp_B1.pickle", "rb"))
-        exp_B2 = pickle.load(open("exp_B2.pickle", "rb"))
-        exp_B3 = pickle.load(open("exp_B3.pickle", "rb"))
-        exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
-        exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
-        exp_C3 = pickle.load(open("exp_C3.pickle", "rb"))
-        exp_D = pickle.load(open("exp_D.pickle", "rb"))
-
-        exp_A_noPE = pickle.load(open("exp_A_noPE.pickle", "rb"))
-        exp_B1_noPE = pickle.load(open("exp_B1_noPE.pickle", "rb"))
-        exp_B2_noPE = pickle.load(open("exp_B2_noPE.pickle", "rb"))
-        exp_B3_noPE = pickle.load(open("exp_B3_noPE.pickle", "rb"))
-        exp_C1_noPE = pickle.load(open("exp_C1_noPE.pickle", "rb"))
-        exp_C2_noPE = pickle.load(open("exp_C2_noPE.pickle", "rb"))
-        exp_C3_noPE = pickle.load(open("exp_C3_noPE.pickle", "rb"))
-        exp_D_noPE = pickle.load(open("exp_D_noPE.pickle", "rb"))
-
-        experiments_1 = [exp_A,exp_B1,exp_B2,exp_B3,exp_C1,exp_C2,exp_C3,exp_D]
-        experiments_2 = [exp_A_noPE,exp_B1_noPE,exp_B2_noPE,exp_B3_noPE,exp_C1_noPE,exp_C2_noPE,exp_C3_noPE,exp_D_noPE]
-        labels = ["Exp. A","Exp. B1","Exp. B2","Exp. B3","Exp. C1","Exp. C2","Exp. C3", "Exp. D"]
-        linecolors = p.rcParams['axes.prop_cycle'].by_key()['color']
-
-        plotfig(range(4),  'fig06.png')
-
-    if fig07:
-
-        exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
-        exp_C3 = pickle.load(open("exp_C3.pickle", "rb"))
-
-        npar2 = exp_C2.nr.parameter_number
-        npar3 = exp_C3.nr.parameter_number
-        updates_C2 = exp_C2.da.analyses[:,:-npar2,:].mean(axis=2).T - exp_C2.da.backgrounds[:,:-npar2,:].mean(axis=2).T
-        updates_C3 = exp_C3.da.analyses[:,:-npar3,:].mean(axis=2).T - exp_C3.da.backgrounds[:,:-npar3,:].mean(axis=2).T
-
-        fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
-        fig.set_size_inches(6,3)
-        c1 = ax1.pcolormesh(exp_C2.da.time/3600,exp_C2.nr.zt,updates_C2,vmin=-0.03,vmax=0.03,cmap='RdBu_r')
-        c2 = ax2.pcolormesh(exp_C3.da.time/3600,exp_C2.nr.zt,updates_C3,vmin=-0.03,vmax=0.03,cmap='RdBu_r')
-        #
-        ax1.set_title(r'a) Exp. C$_2$, $\theta_a-\theta_b$')
-        ax2.set_title(r'b) Exp. C$_3$, $\theta_a-\theta_b$')
-        ax1.set_xlabel('Time (h)')
-        ax2.set_xlabel('Time (h)')
-        ax1.set_ylabel('Height (m)')
-        ax2.sharey(ax1)
-        p.setp(ax2.get_yticklabels(), visible=False)
-        #
-        p.colorbar(c1,orientation='horizontal')
-        p.colorbar(c2,orientation='horizontal')
-        #
         fig.savefig('fig07.png',format='png',dpi=300)
         p.close(fig)