diff --git a/PE_CBL.py b/PE_CBL.py
index 9c7ef00fd5a6271e5563ffe3a1026eb122af09c9..788a9a0436bd15fca7431c693dd7485775f27ee8 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -55,6 +55,8 @@ if __name__ == '__main__':
         'perturbations_uv_amplitude' : 0.1,
         'perturbations_smooth_number' : 11,
         'perturbations_symmetric' : True,
+        'simulate_error_growth' : False,
+        'error_growth_perturbations_amplitude' : 0.5,
         # Part 2: perturbation of parameters
         'perturb_ensemble_parameters' : True,
         'parameter_number' : 1,
@@ -88,7 +90,7 @@ if __name__ == '__main__':
         "Assimilation interval must be an integer multiplier of model dt"
 
     # Decide what figures to plot
-    fig01 = False
+    fig01 = True
     fig02 = True
     fig03 = True
     fig04 = True
@@ -550,7 +552,7 @@ if __name__ == '__main__':
         cbl_settings_E['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
         cbl_settings_E['perturbations_theta_amplitude'] = sigma_b_init*10
         cbl_settings_E['Hmax'] = 0.15
-        cbl_settings_E['H0_perturbation_ampl_init'] = 0.15
+        cbl_settings_E['simulate_error_growth'] = True
         cbl_settings_E['is_cgrad'] = False
         da_settings_E['cbl_settings'] = cbl_settings_E
         da_settings_E['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*10
@@ -598,6 +600,38 @@ if __name__ == '__main__':
 
     if fig07:
 
+        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"))
@@ -616,50 +650,13 @@ if __name__ == '__main__':
         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]
-        experiments_2 = [exp_A_noPE,exp_B1_noPE,exp_B2_noPE,exp_B3_noPE,exp_C1_noPE,exp_C2_noPE,exp_C3_noPE]
-        labels = ["exp_A","exp_B1","exp_B2","exp_B3","exp_C1","exp_C2","exp_C3"]
-        #experiments_1 = [exp_A,exp_D]
-        #experiments_2 = [exp_A_noPE,exp_D_noPE]
-        #labels = ["exp_A","exp_D"]
+        experiments_1 = [exp_A,exp_B1,exp_B2,exp_B3,exp_C1,exp_C2,exp_C3,exp_D,exp_E]
+        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,exp_E_noPE]
+        labels = ["Exp. A","Exp. B1","Exp. B2","Exp. B3","Exp. C1","Exp. C2","Exp. C3", "Exp. D", "Exp. E"]
         linecolors = p.rcParams['axes.prop_cycle'].by_key()['color']
-        #["red","blue","green","gold","black","orange"]
 
-        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 range(len(experiments_1)):
-        for i in range(4):
-            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)')
-        #
-        #ax1.set_xscale('log')
-        #ax2.set_xscale('log')
-        #ax3.set_xscale('log')
-        #
-        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('fig07.png',format='png',dpi=300)
-        p.close(fig)
+        plotfig(range(4),  'fig07a.png')
+        plotfig(range(4,9),'fig07b.png')
 
     if fig08: