diff --git a/observations.py b/observations.py
index 6d83d25d7ac4864a45f895714d208e1f3ba5a8da..c5582c13c133f6964730b46147d8d1ef8c7f4e86 100644
--- a/observations.py
+++ b/observations.py
@@ -2,70 +2,11 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 from matplotlib import pyplot as p
-import pickle
-import os
 import glob
 import pandas as pd
 
-from models import CBL
-
 is_sorted = lambda a: np.all(a[:-1] <= a[1:])
 
-# Default CBL model settings
-default_cbl_settings ={
-    # Physical parameters
-    'g' : 9.806,
-    'f' : 1e-4,
-    'kvonk' : 0.4,
-    'z0' : 0.1,
-    'ug' : 0,
-    'vg' : 10,
-    'theta_0' : 290,
-    'gamma' : 3e-3,
-    'Hmax' : 0.1,
-    'H0_perturbation_ampl_init' : 0.0,
-    'H0_perturbation_ampl_time' : 0.0,
-    'exct' : 0.3,
-    'pfac' : 0.8,
-    'Kmin' : 0.1,
-    'Kmax' : 200,
-    # Model formulation
-    'is_bwind' : False,
-    'is_cgrad' : True,
-    'is_cycle' : False,
-    # Domain parameters (Kmax determines dt)
-    'dz' : 50,
-    'ztop' : 4000,
-    # Ensemble generation parameters (only effective for ensemble runs)
-    # Part 1: perturbation of initial state
-    # perturbation_type can be "smooth", "random", "uniform"
-    # the adjective refers to vertical variability
-    'rnseed' : 181612,
-    'perturb_ensemble_state' : True,
-    'perturbations_type' : "smooth",
-    'perturbations_theta_amplitude' : 0.1,
-    'perturbations_uv_amplitude' : 0.1,
-    'perturbations_smooth_number' : 11,
-    'perturbations_symmetric' : True,
-    'simulate_error_growth' : False,
-    'error_growth_perturbations_amplitude' : 0.0,
-    # Part 2: perturbation of parameters
-    'perturb_ensemble_parameters' : True,
-    'parameter_number' : 1,
-    'parameter_transform' : None,
-    'parameter_ensemble_min' : np.array([0.5]),
-    'parameter_ensemble_max' : np.array([4.5]),
-    'parameter_true' : np.array([1.5]),
-    # In case of parameter estimation
-    'do_parameter_estimation' : False,
-    'parameter_inflation_rtps_alpha' : np.array([0.8]),
-    'return_covariances_increments_and_innovations' : True
-    }
-
-def observation_operator(fp,xp,x):
-    f = np.interp(x, xp, fp)
-    return f
-
 def decode_observations(observation_files):
 
     # Read set of CSV files into a pandas dataframe.
@@ -140,132 +81,52 @@ if __name__ == '__main__':
 
     printout = False
 
-    CASE = 'energy_diagnostics'
-
-    if CASE == 'display_obs':
-        # Read the observations
-        observations,ocoords = decode_observations('./LES_data/Theta/*csv')
-
-        if printout:
-            ntim = 0
-            for i in range(observations.shape[1]):
-                print(ocoords[ntim,i],observations[ntim,i])
-
-        # Sort the observations randomly
-        observations,ocoords = random_sort_observations(observations,ocoords)
-        if printout:
-            for i in range(observations.shape[1]):
-                print(ocoords[ntim,i],observations[ntim,i])
-
-        makeplots = True
-        # Make plots
-        if makeplots:
-
-            # If observations are shuffled, sort them again
-            ntimes = ocoords.shape[0]
-            for k in range(ntimes):
-                if is_sorted(ocoords[k,:]):
-                    pass
-                else:
-                    indices = np.argsort(ocoords[k,:])
-                    ocoords[k,:] = ocoords[k,indices]
-                    observations[k,:] = observations[k,indices]
-
-            fig, ax1 = p.subplots(1,1,constrained_layout=True)
-            fig.set_size_inches(4,3)
-
-            tint = 300
-            t = np.linspace(0,tint*(ocoords.shape[0]-1),ocoords.shape[0])[:,None]+np.zeros(ocoords.shape)
-            t = t/3600.
-            #c = ax1.scatter(t,ocoords,s=5,c=observations,vmin=290,vmax=296)
-            c = ax1.pcolormesh(t,ocoords,observations,vmin=290,vmax=296)
-            ax1.contour(t,ocoords,observations,
-                    np.linspace(290,290+0.003*2000,13),
-                    colors='white',
-                    linestyles='--',
-                    linewidths=0.75)
-            ax1.set_ylim([0,2000])
-            ax1.set_ylabel(r'Height (m)')
-            ax1.set_xlabel(r'Time (s)')
-            ax1.set_title(r'$\overline{\theta}$ (K) obs in t=(%u-%u) h'%(t.min(),t.max()))
-            ax1.legend(frameon=False)
-            p.colorbar(c,orientation='horizontal')
-            fig.savefig('LES_obs.png',format='png',dpi=300)
-            p.close(fig)
+    # Read the observations
+    observations,ocoords = decode_observations('./LES_data/Theta/*csv')
+
+    if printout:
+        ntim = 0
+        for i in range(observations.shape[1]):
+            print(ocoords[ntim,i],observations[ntim,i])
+
+    # Sort the observations randomly
+    observations,ocoords = random_sort_observations(observations,ocoords)
+    if printout:
+        for i in range(observations.shape[1]):
+            print(ocoords[ntim,i],observations[ntim,i])
+
+    makeplots = True
+    # Make plots
+    if makeplots:
+
+        # If observations are shuffled, sort them again
+        ntimes = ocoords.shape[0]
+        for k in range(ntimes):
+            if is_sorted(ocoords[k,:]):
+                pass
+            else:
+                indices = np.argsort(ocoords[k,:])
+                ocoords[k,:] = ocoords[k,indices]
+                observations[k,:] = observations[k,indices]
+
+        fig, ax1 = p.subplots(1,1,constrained_layout=True)
+        fig.set_size_inches(4,3)
 
-    if CASE == 'energy_diagnostics':
-
-        # Read the observations
-        theta,z = decode_observations('./LES_data/Theta/*csv')
-        nassim,nobs = theta.shape
-        thicknesses = np.zeros(z.shape)+np.nan
-        for k in range(nassim):
-            thicknesses[k,:]=np.hstack((z[k,0],np.diff(z[k,:])))
         tint = 300
-        t = np.linspace(0,tint*(z.shape[0]-1),z.shape[0])[:,None]+np.zeros(z.shape)
-
-        # Compute delta theta wrt IC
-        theta_init = 290+z*0.003
-        delta_theta = theta-theta_init
-
-        # Set the sensible heat flux (kinematic units)
-        hfx = 0.12
-
-        # Run model
-        run_settings = dict(default_cbl_settings)
-        run_settings['theta_0'] = 290
-        run_settings['gamma'] = 3e-3
-        run_settings['Hmax'] = 0.12
-        run = CBL(run_settings)
-        run.maxtime = 25200
-        run.initialize(1) 
-        run.run(output_full_history=True)
-        zt = run.zt
-
-        # Model state
-        model_theta = run.history['theta']
-        model_delta_theta = model_theta-model_theta[:,0][:,None]
-        model_thicknesses = run.dz*(np.ones(zt.size))[:,None]+np.zeros(model_theta.shape)
-        model_delta_theta = model_delta_theta.T
-        model_thicknesses = model_thicknesses.T
-        model_times = run.history['time']
-
-        # Model equivalents
-        model_equivalents = np.zeros(theta.shape)+np.nan
-        for i in range(int(nassim)):
-            valid_time = i*tint
-            time_index = np.argwhere(run.history['time'] == valid_time)[0][0]
-            for j in range(nobs):
-                model_equivalents[i,j] = observation_operator(model_theta[:,time_index],zt,z[i,j])
-        model_e_delta_theta = model_equivalents-theta_init
-        model_e_thicknesses = thicknesses
-        model_e_times = t
-
-        # Make plots
-        fig, [ax1,ax2] = p.subplots(1,2,constrained_layout=True)
-        fig.set_size_inches(6,3)
-
-        # Time coordinate is hours
+        t = np.linspace(0,tint*(ocoords.shape[0]-1),ocoords.shape[0])[:,None]+np.zeros(ocoords.shape)
         t = t/3600.
-
-        c = ax1.pcolormesh(t,z,delta_theta,vmin=-3,vmax=3)
-        ax1.contour(t,z,theta,
+        #c = ax1.scatter(t,ocoords,s=5,c=observations,vmin=290,vmax=296)
+        c = ax1.pcolormesh(t,ocoords,observations,vmin=290,vmax=296)
+        ax1.contour(t,ocoords,observations,
                 np.linspace(290,290+0.003*2000,13),
                 colors='white',
                 linestyles='--',
                 linewidths=0.75)
         ax1.set_ylim([0,2000])
-        #ax1.set_ylabel(r'Height (m)')
-        #ax1.set_xlabel(r'Time (s)')
-        #ax1.set_title(r'$\overline{\theta}$ (K) obs in t=(%u-%u) h'%(t.min(),t.max()))
+        ax1.set_ylabel(r'Height (m)')
+        ax1.set_xlabel(r'Time (s)')
+        ax1.set_title(r'$\overline{\theta}$ (K) obs in t=(%u-%u) h'%(t.min(),t.max()))
         ax1.legend(frameon=False)
         p.colorbar(c,orientation='horizontal')
-
-        ax2.plot(t[:,0],t[:,0]*3600*hfx,color='k',dashes=[3,1],zorder=10,label=r'$\int H dt$')
-        ax2.plot(model_times/3600,np.sum(model_delta_theta*model_thicknesses,axis=1),color='r',label=r'$\int\Delta\theta dz$, model state')
-        ax2.plot(model_e_times[:,0]/3600,np.sum(model_e_delta_theta*model_e_thicknesses,axis=1),color='orange',label=r'$\int\Delta\theta dz$, model equivalents')
-        ax2.plot(t[:,0],np.sum(delta_theta*thicknesses,axis=1),label=r'$\int\Delta\theta dz$, observations')
-        ax2.text(3.5,0,'RMSD=%6.4f K'%np.mean((theta-model_equivalents)**2))
-        ax2.legend(prop={'size': 6})
-        fig.savefig('LES_delta_theta.png',format='png',dpi=300)
+        fig.savefig('LES_obs.png',format='png',dpi=300)
         p.close(fig)