old Fig 6 deleted

PE_CBL.py

@@ -19,6 +19,9 @@ if __name__ == '__main__':
     sigma_o_gn = 0.0 # set equal to sigma_o_as in experiments C and F
     sigma_o_as = 0.1 # increased by a factor 10 in experiments C and F
 
+    # Initiate a random number generator with preset seed, to ensure reproducibility
+    RNG = np.random.default_rng(seed=181612)
+
     # Default CBL model settings
     default_cbl_settings ={
         # Physical parameters
@@ -49,14 +52,15 @@ if __name__ == '__main__':
         # Part 1: perturbation of initial state
         # perturbation_type can be "smooth", "random", "uniform"
         # the adjective refers to vertical variability
+        'RNG' : RNG,
         'perturb_ensemble_state' : True,
-        'perturbations_type' : "uniform",
+        'perturbations_type' : "smooth",
         'perturbations_theta_amplitude' : sigma_b_init,
         'perturbations_uv_amplitude' : 0.1,
         'perturbations_smooth_number' : 11,
         'perturbations_symmetric' : True,
         'simulate_error_growth' : False,
-        'error_growth_perturbations_amplitude' : 0.5,
+        'error_growth_perturbations_amplitude' : 0.0,
         # Part 2: perturbation of parameters
         'perturb_ensemble_parameters' : True,
         'parameter_number' : 1,
@@ -83,6 +87,7 @@ if __name__ == '__main__':
             'nens' : nens,
             'FILTER' : 'EAKF',
             'inflation_rtps_alpha' : 0.2,
+            'RNG' : RNG
             }
 
     # Safety check on time steps
@@ -497,8 +502,10 @@ if __name__ == '__main__':
         cbl_settings_D['perturbations_theta_amplitude'] = sigma_b_init*10
         cbl_settings_D['Hmax'] = 0.15
         cbl_settings_D['is_cgrad'] = False
+        cbl_settings_D['simulate_error_growth'] = True
+        cbl_settings_D['error_growth_perturbations_amplitude'] = sigma_b_init*5
         da_settings_D['cbl_settings'] = cbl_settings_D
-        da_settings_D['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*10
+        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
@@ -544,65 +551,6 @@ if __name__ == '__main__':
 
     if fig06:
 
-        # Create a copy of the default settings
-        cbl_settings_E = dict(default_cbl_settings)
-        da_settings_E = dict(default_da_settings)
-
-        # Change settings as necessary
-        # Changes include generation of observations, so the existing nature run
-        # can't be reused.
-        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['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
-        da_settings_E['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
-
-        # Run
-        exp_E = experiment(da_settings_E)
-        setattr(exp_E,'label','E')
-
-        # Run the corresponding experiment without parameter estimation
-        cbl_settings_E_noPE = dict(cbl_settings_E)
-        da_settings_E_noPE = dict(da_settings_E)
-        cbl_settings_E_noPE['do_parameter_estimation'] = False
-        da_settings_E_noPE['cbl_settings'] = cbl_settings_E_noPE
-        exp_E_noPE = experiment(da_settings_E_noPE)
-
-        # Save experiments to disk
-        pickle.dump(exp_E, open('exp_E.pickle', 'wb')) 
-        pickle.dump(exp_E_noPE, open('exp_E_noPE.pickle', 'wb')) 
-
-        # Make plots
-        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_E,da_settings_E['obs_error_sdev_assimilate'][0],None,ax=[ax0,ax1,ax2])
-        ax0.set_title(r'a) Exp. E, cov($p,y_b}$) (K)')
-        ax0.set_xlabel('Time (h)')
-        ax0.set_ylabel('Height (m)')
-        ax1.set_title(r'b) Exp. E, $\sigma^2_{y^b}}$ (K)')
-        ax1.set_xlabel('Time (h)')
-        ax1.set_ylabel('Height (m)')
-        ax2.set_title(r'c) Exp. E, $K_{p,y_b}$ (K$^{-1}$)')
-        ax2.set_xlabel('Time (h)')
-        ax2.set_ylabel('Height (m)')
-        ax3 = plot_CBL_PE(exp_E,None,ax=ax3)
-        ax3.set_title(r'd) Exp. E, evolution of $p$')
-        ax3.set_xlabel('Time (h)')
-        ax3.set_yticks([0,1,2,3,4,5])
-        p.colorbar(c0,orientation='horizontal')
-        p.colorbar(c1,orientation='horizontal')
-        p.colorbar(c2,orientation='horizontal')
-        #
-
-        fig.savefig('fig06.png',format='png',dpi=300)
-        p.close(fig)
-
-    if fig07:
-
         def plotfig(exprange,filename):
             fig, [ax1, ax2, ax3] = p.subplots(1,3,constrained_layout=True)
             fig.set_size_inches(8,3)
@@ -653,15 +601,14 @@ 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,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"]
+        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),  'fig07a.png')
-        plotfig(range(4,9),'fig07b.png')
+        plotfig(range(4),  'fig06.png')
 
-    if fig08:
+    if fig07:
 
         exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
         exp_C3 = pickle.load(open("exp_C3.pickle", "rb"))
@@ -687,7 +634,7 @@ if __name__ == '__main__':
         p.colorbar(c1,orientation='horizontal')
         p.colorbar(c2,orientation='horizontal')
         #
-        fig.savefig('fig08.png',format='png',dpi=300)
+        fig.savefig('fig07.png',format='png',dpi=300)
         p.close(fig)
 
     if opt01: