diff --git a/PE_CBL.py b/PE_CBL.py
index 9cf0c58006b5fd242a0159d162954c2fb402b8d4..6b50f42731a3f070fa309053e81e69e9dbf156c1 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -731,10 +731,91 @@ if __name__ == '__main__':
 
     if fig09:
 
-        path_EAKF = "./runs/seed=181612_enssize=200_EAKF_6hrs"
-        path_LETKF = "./runs/seed=181612_enssize=200_LETKF_6hrs"
-
-
+        path_EAKF = "./runs/seed=181612_enssize=200_EAKF_6hrs/"
+        path_LETKF = "./runs/seed=181612_enssize=200_LETKF_6hrs/"
+
+        exp_Aa = pickle.load(open(path_EAKF+"exp_A.pickle", "rb"))
+        exp_B1a = pickle.load(open(path_EAKF+"exp_B1.pickle", "rb"))
+        exp_B2a = pickle.load(open(path_EAKF+"exp_B2.pickle", "rb"))
+        exp_B3a = pickle.load(open(path_EAKF+"exp_B3.pickle", "rb"))
+        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_Da = pickle.load(open(path_EAKF+"exp_D.pickle", "rb"))
+
+        exp_Ab = pickle.load(open(path_LETKF+"exp_A.pickle", "rb"))
+        exp_B1b = pickle.load(open(path_LETKF+"exp_B1.pickle", "rb"))
+        exp_B2b= pickle.load(open(path_LETKF+"exp_B2.pickle", "rb"))
+        exp_B3b = pickle.load(open(path_LETKF+"exp_B3.pickle", "rb"))
+        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_Db = pickle.load(open(path_LETKF+"exp_D.pickle", "rb"))
+
+        fig, [ax1, ax2] = p.subplots(1,2)
+        fig.set_size_inches(6,3)
+        #
+        cc = 0
+        linecolors = ['darkviolet','chocolate','darkorange','orange','gold','darkgreen','seagreen','royalblue']
+        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']
+        for exp in [exp_Aa,exp_B1a,exp_B2a,exp_B3a,exp_B4a,exp_C1a,exp_C2a,exp_Da]:
+            # Experiment settings
+            true_value = exp.nr.pfac
+            parameter_number = exp.nr.parameter_number
+            parameter_id = 0
+            pt = exp.nr.parameter_transform
+            par_tran = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
+            par_phys = pt[parameter_id](par_tran,kind='inv')
+            initpar_phys = pt[parameter_id](exp.da.initial_perturbed_parameters,kind='inv')
+            t = exp.da.time/3600.
+
+            # Initial parameter distribution
+            ax1.step([0,t[0]],np.median(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc],label=labels[cc])
+
+            # Later times
+            ax1.step(t,np.median(par_phys,axis=1),color=linecolors[cc])
+            cc += 1
+
+        cc = 0
+        for exp in [exp_Ab,exp_B1b,exp_B2b,exp_B3b,exp_B4b,exp_C1b,exp_C2b,exp_Db]:
+            # Experiment settings
+            true_value = exp.nr.pfac
+            parameter_number = exp.nr.parameter_number
+            parameter_id = 0
+            pt = exp.nr.parameter_transform
+            par_tran = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
+            par_phys = pt[parameter_id](par_tran,kind='inv')
+            initpar_phys = pt[parameter_id](exp.da.initial_perturbed_parameters,kind='inv')
+            t = exp.da.time/3600.
+
+            # Initial parameter distribution
+            ax2.step([0,t[0]],np.median(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc])
+
+            # Later times
+            ax2.step(t,np.median(par_phys,axis=1),color=linecolors[cc])
+            cc += 1
+
+        # Other stuff
+        ax1.axhline(y=true_value,linestyle='--',color='black')
+        ax1.set_title(r'a) EAKF experiments')
+        ax1.set_xlim([0,t.max()])
+        ax1.set_xlabel('Time (h)')
+        ax1.set_ylabel(r'$p$')
+        ax1.set_yticks([1,2,3])
+        #
+        ax2.axhline(y=true_value,linestyle='--',color='black')
+        ax2.set_title(r'b) LETKF experiments')
+        ax2.set_xlim([0,t.max()])
+        ax2.set_xlabel('Time (h)')
+        ax2.set_ylabel(r'$p$')
+        ax2.set_yticks([1,2,3])
+        #
+        ax2.sharey(ax1)
+        fig.legend(loc=8,ncol=4,frameon=False)
+        fig.subplots_adjust(bottom=0.35, wspace=0.38)
+        #
+        fig.savefig('fig09.png',format='png',dpi=300)
+        p.close(fig)
 
     if opt01: