diff --git a/PE_CBL.py b/PE_CBL.py
index 56e9a99fd483d33f2836f66823b88e7c20862034..72a584b0d782df10a44cd130c67f37e186b84dfa 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -8,7 +8,7 @@ from copy import deepcopy
 
 # Own modules
 from models import CBL
-from ENDA import experiment, experiment_real, transform_log, transform_none, transform_pit
+from ENDA import experiment, transform_log, transform_none, transform_pit
 from graphics import *
 
 if __name__ == '__main__':
@@ -80,6 +80,7 @@ if __name__ == '__main__':
             'tspinup' : 3600,
             'trun' : 21600,
             'assimilation_interval' : 300,
+            'randomize_obs': True,
             'obs_coordinates' : np.linspace(0,2000,nobs),
             'obs_kinds' : np.array(['theta']*nobs),
             'obs_error_sdev_generate' : np.ones(nobs)*0,
@@ -863,8 +864,8 @@ if __name__ == '__main__':
         fig.savefig('profiles.png',format='png',dpi=300)
         p.close(fig)
 
-    # Test with "real" observations (no nature run)
     if opt03:
+        # Test with "real" observations (no nature run)
 
         nens = 20
         ############################################################################
@@ -887,7 +888,7 @@ if __name__ == '__main__':
         try:
             exp_RC = pickle.load(open(da_settings_RC['path']+"exp_RC.pickle", "rb"))
         except:
-            exp_RC = experiment_real(da_settings_RC)
+            exp_RC = experiment(da_settings_RC)
             setattr(exp_RC,'label','RC')
             pickle.dump(exp_RC, open(da_settings_RC['path']+'exp_RC.pickle', 'wb')) 
 
@@ -914,7 +915,7 @@ if __name__ == '__main__':
         try:
             exp_R = pickle.load(open(da_settings_R['path']+'exp_R.pickle', "rb"))
         except:
-            exp_R = experiment_real(da_settings_R)
+            exp_R = experiment(da_settings_R)
             setattr(exp_R,'label','R')
             pickle.dump(exp_R, open(da_settings_R['path']+'exp_R.pickle', 'wb')) 
 
@@ -1004,7 +1005,7 @@ if __name__ == '__main__':
         p.close(fig)
         
         #########################################################################################à
-        # Make plots: perfect-model experiment (R)
+        # Make plots: imperfect-model experiment (R)
         fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
         fig.set_size_inches(6,3)
 
@@ -1061,6 +1062,96 @@ if __name__ == '__main__':
         fig.savefig('expR_PE.png',format='png',dpi=300)
         p.close(fig)
 
+    if opt04:
+
+        # Use only 20 ensemble members for these experiments
+        nens = 20
+
+        # Base settings
+        cbl_settings_base = dict(default_cbl_settings)
+        cbl_settings_base['theta_0'] = 290
+        cbl_settings_base['gamma'] = 3e-3
+        cbl_settings_base['Hmax'] = 0.12
+        da_settings_base = dict(default_da_settings)
+        da_settings_base['cbl_settings'] = dict(cbl_settings_base)
+        da_settings_base['nens'] = nens
+        da_settings_base['type'] = 'real'
+        da_settings_base['observation_files'] = './LES_data/Theta/*csv'
+        da_settings_base['obs_coordinates'] = None
+        da_settings_base['obs_error_sdev_generate'] = None
+
+        # Base experiment (in opt03 has tspinup = 1 hour). Try other lengths.
+        da_settings_spinup_1 = dict(da_settings_base)
+        da_settings_spinup_1['tspinup'] = 7200
+        da_settings_spinup_1['trun'] = 25200-7200
+        da_settings_spinup_1['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
+
+        da_settings_spinup_2 = dict(da_settings_base)
+        da_settings_spinup_2['tspinup'] = 10800
+        da_settings_spinup_2['trun'] = 25200-10800
+        da_settings_spinup_2['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
+
+        da_settings_spinup_3 = dict(da_settings_base)
+        da_settings_spinup_3['tspinup'] = 14400
+        da_settings_spinup_3['trun'] = 25200-14400
+        da_settings_spinup_3['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
+
+        # Run and save to disk
+        try:
+            exp_spinup_1 = pickle.load(open(da_settings_spinup_1['path']+'exp_spinup_1.pickle', "rb"))
+        except:
+            exp_spinup_1 = experiment(da_settings_spinup_1)
+            setattr(exp_spinup_1,'label','spinup_1')
+            pickle.dump(exp_spinup_1, open(da_settings_spinup_1['path']+'exp_spinup_1.pickle', 'wb')) 
+        try:
+            exp_spinup_2 = pickle.load(open(da_settings_spinup_2['path']+'exp_spinup_2.pickle', "rb"))
+        except:
+            exp_spinup_2 = experiment(da_settings_spinup_2)
+            setattr(exp_spinup_2,'label','spinup_2')
+            pickle.dump(exp_spinup_2, open(da_settings_spinup_2['path']+'exp_spinup_2.pickle', 'wb')) 
+        try:
+            exp_spinup_3 = pickle.load(open(da_settings_spinup_3['path']+'exp_spinup_3.pickle', "rb"))
+        except:
+            exp_spinup_3 = experiment(da_settings_spinup_3)
+            setattr(exp_spinup_3,'label','spinup_3')
+            pickle.dump(exp_spinup_3, open(da_settings_spinup_3['path']+'exp_spinup_3.pickle', 'wb')) 
+
+        for exp,settings in zip([exp_spinup_1,exp_spinup_2,exp_spinup_3],[da_settings_spinup_1,da_settings_spinup_2,da_settings_spinup_3]):
+            label = exp.label
+            maxtime = exp.trun
+
+            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,settings['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
+            ax0.set_title(r'a) Exp. %s, $\rho(p\prime\prime,y_b}$)'%label)
+            ax0.set_xlabel('Time (h)')
+            ax0.set_ylabel('Height (m)')
+            ax1.set_title(r'b) Exp. %s, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$'%label)
+            ax1.set_xlabel('Time (h)')
+            ax1.set_ylabel('Height (m)')
+            ax2.set_title(r'c) Exp. %s, $\delta\overline{p}\prime\prime$'%label)
+            ax2.set_xlabel('Time (h)')
+            ax2.set_ylabel('Height (m)')
+            ax3 = plot_CBL_PE(exp,None,ax=ax3)
+            ax3.set_title(r'd) Exp. %s, evolution of $p$'%label)
+            ax3.set_xlabel('Time (h)')
+            ax3.set_yticks([0,1,2,3,4,5])
+            ax0.set_xlim([0,maxtime/3600])
+            ax1.set_xlim([0,maxtime/3600])
+            ax2.set_xlim([0,maxtime/3600])
+            ax3.set_xlim([0,maxtime/3600])
+            ax0.set_xticks(np.arange((maxtime+1)/3600))
+            ax1.set_xticks(np.arange((maxtime+1)/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')
+            #
+            fig.savefig('exp_%s_PE.png'%label,format='png',dpi=300)
+            p.close(fig)
+
     check_reproducibility = False
     if check_reproducibility:
         path = default_da_settings['path']