diff --git a/PE_CBL.py b/PE_CBL.py
index f348c9ca395331fd027ffc0c57e258904141cadc..750dc0a7783f5b75b54fdfeeae046cceec8fa278 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -97,12 +97,12 @@ if __name__ == '__main__':
 
     # Decide what figures to plot
     fig01 = False
-    fig02 = True
-    fig03 = True
-    fig04 = True
-    fig05 = True
-    fig06 = True
-    fig07 = True
+    fig02 = False
+    fig03 = False
+    fig04 = False
+    fig05 = False
+    fig06 = False
+    fig07 = False
     fig08 = True
     
     # Other possible, optional, plots
@@ -509,17 +509,19 @@ if __name__ == '__main__':
 
     if fig02:
 
+        exp_A = pickle.load(open("exp_A.pickle", "rb"))
+
         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_A,da_settings_A['obs_error_sdev_assimilate'][0],None,ax=[ax0,ax1,ax2])
+        [ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
         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, $\delta y\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
+        ax1.set_title(r'b) Exp. A, $\delta\overline{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, $\delta p\prime\prime$')
+        ax2.set_title(r'c) Exp. A, $\delta\overline{p}\prime\prime$')
         ax2.set_xlabel('Time (h)')
         ax2.set_ylabel('Height (m)')
         ax3 = plot_CBL_PE(exp_A,None,ax=ax3)
@@ -546,6 +548,11 @@ if __name__ == '__main__':
 
     if fig04:
 
+        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_B4 = pickle.load(open("exp_B4.pickle", "rb"))
+
         fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
         fig.set_size_inches(6,4)
         #
@@ -579,6 +586,9 @@ if __name__ == '__main__':
 
     if fig05:
 
+        exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
+        exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
+
         fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
         fig.set_size_inches(6,2)
         #
@@ -613,8 +623,10 @@ if __name__ == '__main__':
         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_ylim([0,2000])
+        ax2.set_ylim([0,2000])
+        ax1.set_title(r'a) Exp. C$_1$, $\overline{\theta}_a-\overline{\theta}_b$')
+        ax2.set_title(r'b) Exp. C$_2$, $\overline{\theta}_a-\overline{\theta}_b$')
         ax1.set_xlabel('Time (h)')
         ax2.set_xlabel('Time (h)')
         ax1.set_ylabel('Height (m)')
@@ -629,43 +641,54 @@ if __name__ == '__main__':
 
     if fig07:
 
+        exp_D = pickle.load(open("exp_D.pickle", "rb"))
+
+        fig, [ax0, ax1] = p.subplots(1,2,constrained_layout=True)
+        fig.set_size_inches(6,3)
+        #
+        ax0,c0 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],ax=ax0)
+        ax0.set_title(r'a) Exp. D, $\rho(p\prime\prime,y_b}$)')
+        ax0.set_xlabel('Time (h)')
+        ax0.set_ylabel('Height (m)')
+        ax1 = plot_CBL_PE(exp_D,None,ax=ax1)
+        ax1.set_title(r'b) Exp. D, evolution of $p$')
+        ax1.set_xlabel('Time (h)')
+        ax1.set_yticks([0,1,2,3,4,5])
+        p.colorbar(c0,orientation='horizontal')
+        #
+        fig.savefig('fig07.png',format='png',dpi=300)
+        p.close(fig)
+
+    if fig08:
+
+        exp_A = pickle.load(open("exp_A.pickle", "rb"))
+        exp_D = pickle.load(open("exp_D.pickle", "rb"))
+
         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, $\rho(p\prime\prime,y_b}$)')
+        [ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax1,ax0])
+        ax0.set_title(r'a) Exp. A, $\delta\overline{y}$')
         ax0.set_xlabel('Time (h)')
         ax0.set_ylabel('Height (m)')
-        ax1.set_title(r'b) Exp. D, $\delta y\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
+        ax1.set_title(r'b) Exp. A, $\sigma_{p\prime\prime}/\sigma_{y^b}$')
         ax1.set_xlabel('Time (h)')
         ax1.set_ylabel('Height (m)')
-        ax2.set_title(r'c) Exp. D, $\delta p\prime\prime$')
+        [ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],ax=[ax3,ax2])
+        ax2.set_title(r'c) Exp. D, $\delta\overline{y}$')
         ax2.set_xlabel('Time (h)')
         ax2.set_ylabel('Height (m)')
-        ax3 = plot_CBL_PE(exp_D,None,ax=ax3)
-        ax3.set_title(r'd) Exp. D, evolution of $p$')
+        ax3.set_title(r'd) Exp. D, $\sigma_{p\prime\prime}/\sigma_{y^b}$')
         ax3.set_xlabel('Time (h)')
-        ax3.set_yticks([0,1,2,3,4,5])
+        ax3.set_ylabel('Height (m)')
         p.colorbar(c0,orientation='horizontal')
         p.colorbar(c1,orientation='horizontal')
         p.colorbar(c2,orientation='horizontal')
+        p.colorbar(c3,orientation='horizontal')
         #
-
-        fig.savefig('fig07.png',format='png',dpi=300)
+        fig.savefig('fig08.png',format='png',dpi=300)
         p.close(fig)
 
-    if fig08:
-
-        experiments_pe = [exp_A]
-        experiments_nope = [exp_A_noPE]
-        labels = ["Exp. A"]
-        plot_diagnostics(experiments_pe,experiments_nope,labels,'fig08a.png')
-
-        experiments_pe = [exp_D]
-        experiments_nope = [exp_D_noPE]
-        labels = ["Exp. D"]
-        plot_diagnostics(experiments_pe,experiments_nope,labels,'fig08d.png')
-
     if opt01:
 
         da_settings = {'cbl_settings' : dict(default_cbl_settings),