Skip to content
Snippets Groups Projects
Commit 2b04fb8b authored by Stefano Serafin's avatar Stefano Serafin
Browse files

code cleanup, again

parent 7bdaa7c4
No related branches found
No related tags found
No related merge requests found
...@@ -2,70 +2,11 @@ ...@@ -2,70 +2,11 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import numpy as np import numpy as np
from matplotlib import pyplot as p from matplotlib import pyplot as p
import pickle
import os
import glob import glob
import pandas as pd import pandas as pd
from models import CBL
is_sorted = lambda a: np.all(a[:-1] <= a[1:]) 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): def decode_observations(observation_files):
# Read set of CSV files into a pandas dataframe. # Read set of CSV files into a pandas dataframe.
...@@ -140,9 +81,6 @@ if __name__ == '__main__': ...@@ -140,9 +81,6 @@ if __name__ == '__main__':
printout = False printout = False
CASE = 'energy_diagnostics'
if CASE == 'display_obs':
# Read the observations # Read the observations
observations,ocoords = decode_observations('./LES_data/Theta/*csv') observations,ocoords = decode_observations('./LES_data/Theta/*csv')
...@@ -192,80 +130,3 @@ if __name__ == '__main__': ...@@ -192,80 +130,3 @@ if __name__ == '__main__':
p.colorbar(c,orientation='horizontal') p.colorbar(c,orientation='horizontal')
fig.savefig('LES_obs.png',format='png',dpi=300) fig.savefig('LES_obs.png',format='png',dpi=300)
p.close(fig) 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*(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.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.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)
p.close(fig)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment