diff --git a/models.py b/models.py
index 27e66643a80d55612af0e300767001768beb7301..8e2915a63b892858e2860a35cca2f9271bea5463 100644
--- a/models.py
+++ b/models.py
@@ -92,64 +92,6 @@ class CBL:
         if self.nens>1 and self.perturb_ensemble_state:
             x0 = self.perturb_state(x0)
 
-            """
-            # Random or uniform perturbations perturbations
-            # differ only by magnitude of random sample
-            if self.perturbations_type == "random":
-                randomsize = self.nz
-            elif self.perturbations_type == "uniform":
-                randomsize = 1
-            
-            if self.perturbations_type == "random" or self.perturbations_type == "uniform":
-                ppt = self.RNG.standard_normal((randomsize,self.nens))
-                if self.is_bwind:
-                    ppu = self.RNG.standard_normal((self.nens,randomsize))
-                    ppv = self.RNG.standard_normal((self.nens,randomsize))
-
-            # Smooth perturbations are slightly more complicated
-            if self.perturbations_type == "smooth":
-
-                # Define positions of perturbations
-                # enforce two at boundaries to avoid extrapolation
-                ipert = np.arange(0,self.perturbations_smooth_number-1)*(self.nz//(self.perturbations_smooth_number))
-                ipert = np.append(ipert,self.nz-1)
-
-                # Draw random perturbations, then interpolate
-                randomsize = self.perturbations_smooth_number
-                ppt = np.zeros((self.nz,self.nens))+np.nan
-                pert_t = self.RNG.standard_normal((randomsize,self.nens))
-                for n in range(self.nens):
-                    f = CubicSpline(ipert,pert_t[:,n])
-                    ppt[:,n] = f(np.arange(self.nz))
-                if self.is_bwind:
-                    ppu = np.zeros((self.nz,self.nens))+np.nan
-                    ppv = np.zeros((self.nz,self.nens))+np.nan
-                    pert_u = self.RNG.standard_normal((randomsize,self.nens))
-                    pert_v = self.RNG.standard_normal((randomsize,self.nens))
-                    for n in range(self.nens):
-                        f = CubicSpline(ipert,pert_u[:,n])
-                        ppu[:,n] = f(np.arange(self.nz))
-                        f = CubicSpline(ipert,pert_v[:,n])
-                        ppv[:,n] = f(np.arange(self.nz)) 
-
-            if self.perturbations_symmetric:
-                # if n is odd, the last member remains unperturbed
-                for n in range(self.nens//2):
-                    x0[:self.nz,n] += self.perturbations_theta_amplitude*ppt[:,n]
-                    x0[:self.nz,self.nens//2+n] -= self.perturbations_theta_amplitude*ppt[:,n]
-                    if self.is_bwind:
-                        x0[self.nz:self.nz*2,n] += self.perturbations_uv_amplitude*ppu[:,n]
-                        x0[self.nz:self.nz*2,self.nens//2+n] -= self.perturbations_uv_amplitude*ppu[:,n]
-                        x0[self.nz*2:self.nz*3,n] += self.perturbations_uv_amplitude*ppv[:,n]
-                        x0[self.nz*2:self.nz*3,self.nens//2+n] -= self.perturbations_uv_amplitude*ppv[:,n]
-            else:
-                for n in range(self.nens):
-                    x0[:self.nz,n] += self.perturbations_theta_amplitude*ppt[:,n]
-                    if self.is_bwind:
-                        x0[self.nz:self.nz*2,n] += self.perturbations_uv_amplitude*ppu[:,n]
-                        x0[self.nz*2:self.nz*3,n] += self.perturbations_uv_amplitude*ppv[:,n]
-            """
-
         # If desired, perturb parameters too
         # (last parameter_number elements of state vector)
         # Overwrite the pre-existing unperturbed values