diff --git a/verification.py b/verification.py
new file mode 100644
index 0000000000000000000000000000000000000000..e1e9646f32aaeafabadd181f96e0ab0ce01f0f1c
--- /dev/null
+++ b/verification.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+from matplotlib import pyplot as p
+
+from models import CBL
+from observations import decode_observations
+
+# 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.12,
+ 'H0_perturbation_ampl_init' : 0.0,
+ 'H0_perturbation_ampl_time' : 0.0,
+ 'exct' : 0.3,
+ 'pfac' : 1.5,
+ '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
+
+if __name__ == '__main__':
+
+ # Settings
+ energy_diagnostics = False
+ sensitivity_to_p = True
+
+ # 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)
+
+ if energy_diagnostics:
+
+ # 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 = 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],[ax3,ax4]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,6)
+
+ # Time coordinate is hours
+ t = t/3600.
+
+ # Observations
+ c = ax1.pcolormesh(t,z,delta_theta,vmin=-3,vmax=3)
+ ax1.contour(t,z,theta,
+ 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 (h)')
+ ax1.set_title(r'$\Delta\overline{\theta}$ (K) observations')
+ p.colorbar(c,orientation='horizontal')
+
+ # Observations
+ d = ax2.pcolormesh(t,z,model_e_delta_theta,vmin=-3,vmax=3)
+ ax2.contour(t,z,model_equivalents,
+ np.linspace(290,290+0.003*2000,13),
+ colors='white',
+ linestyles='--',
+ linewidths=0.75)
+ ax2.set_ylim([0,2000])
+ ax2.set_ylabel(r'Height (m)')
+ ax2.set_xlabel(r'Time (h)')
+ ax2.set_title(r'$\Delta\overline{\theta}$ (K) model equivalents')
+ p.colorbar(d,orientation='horizontal')
+
+ # Difference
+ d = ax3.pcolormesh(t,z,model_equivalents-theta,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
+ ax3.set_ylim([0,2000])
+ ax3.set_ylabel(r'Height (m)')
+ ax3.set_xlabel(r'Time (h)')
+ ax3.set_title(r'$\overline{\theta}_{fc}-\overline{\theta}_{obs}$ (K) differences')
+ ax3.text(3.5,100,'RMSD=%6.4f K'%np.mean((theta-model_equivalents)**2))
+ p.colorbar(d,orientation='horizontal')
+
+ # Energy diagnostics
+ ax4.plot(t[:,0],t[:,0]*3600*hfx,color='k',dashes=[3,1],zorder=10,label=r'$\int H dt$')
+ ax4.plot(model_times/3600,np.sum(model_delta_theta*model_thicknesses,axis=1),color='r',label=r'$\int\Delta\theta dz$, model state')
+ ax4.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')
+ ax4.plot(t[:,0],np.sum(delta_theta*thicknesses,axis=1),label=r'$\int\Delta\theta dz$, observations')
+ ax4.set_ylabel(r'$E/(\rho c_p)$ (K m)')
+ ax4.set_xlabel(r'Time (h)')
+ ax4.legend(prop={'size': 6})
+
+ fig.savefig('verification_1.png',format='png',dpi=300)
+ p.close(fig)
+
+ if sensitivity_to_p:
+
+ p_values = np.linspace(0.5,2,31)
+ rmsd = np.zeros(p_values.shape)
+
+ for k in range(p_values.size):
+ print('Run %02u'%k)
+
+ # Run model
+ run_settings = dict(default_cbl_settings)
+ run_settings['pfac'] = p_values[k]
+ run = CBL(run_settings)
+ run.maxtime = 25200
+ run.initialize(1)
+ run.run(output_full_history=True)
+
+ # Model equivalents
+ model_theta = run.history['theta']
+ 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],run.zt,z[i,j])
+
+ rmsd[k] = np.sqrt(np.mean((theta-model_equivalents)**2))
+
+ # Make plots
+ fig, ax1 = p.subplots(1,1,constrained_layout=True)
+ fig.set_size_inches(4,3)
+
+ ax1.plot(p_values,rmsd)
+ ax1.set_xlabel('$p$')
+ ax1.set_ylabel('RMSD (K)')
+
+ fig.savefig('verification_2.png',format='png',dpi=300)
+ p.close(fig)