Select Git revision
test_dart-rttov.py
verification.py 7.05 KiB
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as p
from PE_CBL_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)