diff --git a/observations.py b/observations.py
index dc8e6a71aad832e94754bb75fdf456cc3fef41cb..6d83d25d7ac4864a45f895714d208e1f3ba5a8da 100644
--- a/observations.py
+++ b/observations.py
@@ -7,8 +7,65 @@ 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.
@@ -83,52 +140,132 @@ if __name__ == '__main__':
printout = False
- # 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)
+ 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)
+ 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*(ocoords.shape[0]-1),ocoords.shape[0])[:,None]+np.zeros(ocoords.shape)
+ 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 = 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,
+
+ 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 (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')
- fig.savefig('LES_obs.png',format='png',dpi=300)
+
+ 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)
p.close(fig)