diff --git a/ENDA.py b/ENDA.py
index 4c573296f39f8167cf98000a3fc016642696a97a..5545032abd366f29a47468a48b80ce38a1cf5c8c 100644
--- a/ENDA.py
+++ b/ENDA.py
@@ -353,7 +353,6 @@ class cycle:
             if isinstance(nr,CBL):
                 if cbl_settings["simulate_error_growth"]:
                     x0[:-cbl_settings["parameter_number"],:] += RNG.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
-                    #x0[:-cbl_settings["parameter_number"],:] += np.random.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
 
             # Save initial conditions for next cycle
             da.update(x0)
@@ -447,7 +446,6 @@ class experiment:
                     truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate)
                     observations[i,j] = truths[i,j]+\
                         RNG.normal(0,self.obs_error_sdev_generate[j])
-                        #np.random.normal(0,self.obs_error_sdev_generate[j])
 
             # Store truths and observations
             self.truths = truths
diff --git a/models.py b/models.py
index e9136eceed8117aadce7689f9cc6b708139f2ae3..0e6abe7e072364325dd75523c91bffbe18fe04fe 100644
--- a/models.py
+++ b/models.py
@@ -102,12 +102,9 @@ class CBL:
             
             if self.perturbations_type == "random" or self.perturbations_type == "uniform":
                 ppt = RNG.standard_normal((randomsize,self.nens))
-                #ppt = np.random.randn(randomsize,self.nens)
                 if self.is_bwind:
                     ppu = RNG.standard_normal((self.nens,randomsize))
                     ppv = RNG.standard_normal((self.nens,randomsize))
-                    #ppu = np.random.randn(self.nens,randomsize)
-                    #ppv = np.random.randn(self.nens,randomsize)
 
             # Smooth perturbations are slightly more complicated
             if self.perturbations_type == "smooth":
@@ -121,7 +118,6 @@ class CBL:
                 randomsize = self.perturbations_smooth_number
                 ppt = np.zeros((self.nz,self.nens))+np.nan
                 pert_t = RNG.standard_normal((randomsize,self.nens))
-                #pert_t = np.random.randn(randomsize,self.nens)
                 for n in range(self.nens):
                     f = CubicSpline(ipert,pert_t[:,n])
                     ppt[:,n] = f(np.arange(self.nz))
@@ -130,8 +126,6 @@ class CBL:
                     ppv = np.zeros((self.nz,self.nens))+np.nan
                     pert_u = RNG.standard_normal((randomsize,self.nens))
                     pert_v = RNG.standard_normal((randomsize,self.nens))
-                    #pert_u = np.random.randn(randomsize,self.nens)
-                    #pert_v = np.random.randn(randomsize,self.nens)
                     for n in range(self.nens):
                         f = CubicSpline(ipert,pert_u[:,n])
                         ppu[:,n] = f(np.arange(self.nz))
@@ -222,7 +216,6 @@ class CBL:
 
         for k in range(-self.parameter_number,0):
             dum = RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
-            #dum = np.random.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
             pp[k,:] = self.parameter_transform[k](dum , kind='dir')
 
         return pp
@@ -318,7 +311,6 @@ class CBL:
             
             # Add random perturbations to the initial value
             H0 += RNG.normal(scale=H0_perturbation_ampl_init)
-            #H0 += np.random.normal(scale=H0_perturbation_ampl_init)
 
             # Set the surface momentum flux (ustar)
             ustar = 0
@@ -369,7 +361,6 @@ class CBL:
             # Add time-dependent surface perturbations
             # Then compute sensible heat flux and integrate T equation
             H[0,j]     = H0 + RNG.normal(scale=H0_perturbation_ampl_time)
-            #H[0,j]     = H0 + np.random.normal(scale=H0_perturbation_ampl_time)
             H[1:-1,j]  = -K[1:-1,j]*( (thetap[1:]-thetap[:-1])*rdz - gammac)
             H[nz,j]    = 2*H[nz-1,j]-H[nz-2,j]
             theta[:,j] = thetap[:]-dt*rdz*(H[1:,j]-H[:-1,j])