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

added an experiment that assimilates 'real' observations (from an LES)

parent c48a6483
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 223 additions and 23 deletions
PE_CBL.py
+ 223
− 23
View file @ 1370b3ab
......@@ -8,7 +8,7 @@ from copy import deepcopy
# Own modules
from models import CBL
from ENDA import experiment, transform_log, transform_none, transform_pit
from ENDA import experiment, experiment_real, transform_log, transform_none, transform_pit
from graphics import *
if __name__ == '__main__':
......@@ -417,21 +417,23 @@ if __name__ == '__main__':
########################################################################
# Decide what figures to plot
fig01 = False
fig02 = False
fig03 = False
fig04 = False
fig05 = False
fig06 = False
fig07 = False
fig08 = False
fig09 = True
fig101 = False
fig102 = False
fig103 = False
fig104 = False
fig105 = False
fig106 = False
fig107 = False
fig108 = False
fig109 = False
# Other possible, optional, plots
opt01 = False # assimilation of a single observation
opt02 = False # assimilation of profiles at two times
opt03 = False # assimilation of real observations
opt04 = True # assimilation of real observations, different spinup times
if fig01:
if fig101:
# Create a copy of the default settings
cbl_settings = dict(default_cbl_settings)
......@@ -518,7 +520,7 @@ if __name__ == '__main__':
fig.savefig('fig01.png',format='png',dpi=300)
p.close(fig)
if fig02:
if fig102:
exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
......@@ -555,7 +557,7 @@ if __name__ == '__main__':
fig.savefig('fig02.png',format='png',dpi=300)
p.close(fig)
if fig03:
if fig103:
exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
exp_A_noPE = pickle.load(open(default_da_settings['path']+"exp_A_noPE.pickle", "rb"))
......@@ -566,7 +568,7 @@ if __name__ == '__main__':
plot_diagnostics(experiments_pe,experiments_nope,labels,'fig03.png')
if fig04:
if fig104:
exp_B1 = pickle.load(open(default_da_settings['path']+"exp_B1.pickle", "rb"))
exp_B2 = pickle.load(open(default_da_settings['path']+"exp_B2.pickle", "rb"))
......@@ -604,7 +606,7 @@ if __name__ == '__main__':
fig.savefig('fig04.png',format='png',dpi=300)
p.close(fig)
if fig05:
if fig105:
exp_C1 = pickle.load(open(default_da_settings['path']+"exp_C1.pickle", "rb"))
exp_C2 = pickle.load(open(default_da_settings['path']+"exp_C2.pickle", "rb"))
......@@ -628,7 +630,7 @@ if __name__ == '__main__':
fig.savefig('fig05.png',format='png',dpi=300)
p.close(fig)
if fig06:
if fig106:
exp_C1 = pickle.load(open(default_da_settings['path']+"exp_C1.pickle", "rb"))
exp_C2 = pickle.load(open(default_da_settings['path']+"exp_C2.pickle", "rb"))
......@@ -664,7 +666,7 @@ if __name__ == '__main__':
fig.savefig('fig06.png',format='png',dpi=300)
p.close(fig)
if fig07:
if fig107:
exp_D = pickle.load(open(default_da_settings['path']+"exp_D.pickle", "rb"))
......@@ -689,7 +691,7 @@ if __name__ == '__main__':
fig.savefig('fig07.png',format='png',dpi=300)
p.close(fig)
if fig08:
if fig108:
exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
exp_D = pickle.load(open(default_da_settings['path']+"exp_D.pickle", "rb"))
......@@ -729,7 +731,7 @@ if __name__ == '__main__':
fig.savefig('fig08.png',format='png',dpi=300)
p.close(fig)
if fig09:
if fig109:
path_EAKF = "./runs/seed=181612_enssize=200_EAKF_6hrs/"
path_LETKF = "./runs/seed=181612_enssize=200_LETKF_6hrs/"
......@@ -770,10 +772,10 @@ if __name__ == '__main__':
t = exp.da.time/3600.
# Initial parameter distribution
ax1.step([0,t[0]],np.median(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc],label=labels[cc])
ax1.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc],label=labels[cc])
# Later times
ax1.step(t,np.median(par_phys,axis=1),color=linecolors[cc])
ax1.step(t,np.mean(par_phys,axis=1),color=linecolors[cc])
cc += 1
cc = 0
......@@ -789,10 +791,10 @@ if __name__ == '__main__':
t = exp.da.time/3600.
# Initial parameter distribution
ax2.step([0,t[0]],np.median(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc])
ax2.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc])
# Later times
ax2.step(t,np.median(par_phys,axis=1),color=linecolors[cc])
ax2.step(t,np.mean(par_phys,axis=1),color=linecolors[cc])
cc += 1
# Other stuff
......@@ -861,6 +863,204 @@ if __name__ == '__main__':
fig.savefig('profiles.png',format='png',dpi=300)
p.close(fig)
# Test with "real" observations (no nature run)
if opt03:
nens = 20
############################################################################
# Control experiment (perfect-model)
# Create a copy of the default settings
# and change some of the settings as needed
cbl_settings_RC = dict(default_cbl_settings)
cbl_settings_RC['theta_0'] = 290
cbl_settings_RC['gamma'] = 3e-3
cbl_settings_RC['Hmax'] = 0.12
#cbl_settings_RC['parameter_ensemble_min'] = np.array([0.5])
#cbl_settings_RC['parameter_ensemble_max'] = np.array([2.5])
da_settings_RC = dict(default_da_settings)
da_settings_RC['cbl_settings'] = dict(cbl_settings_RC)
da_settings_RC['nens'] = nens
da_settings_RC['type'] = 'OSSE'
da_settings_RC['path'] = "./runs/real_init_exp_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
# Run and save to disk
try:
exp_RC = pickle.load(open(da_settings_RC['path']+"exp_RC.pickle", "rb"))
except:
exp_RC = experiment_real(da_settings_RC)
setattr(exp_RC,'label','RC')
pickle.dump(exp_RC, open(da_settings_RC['path']+'exp_RC.pickle', 'wb'))
############################################################################
# Imperfect-model experiment
# Create a copy of the default settings
# and change some of the settings as needed
cbl_settings_R = dict(default_cbl_settings)
cbl_settings_R['theta_0'] = 290
cbl_settings_R['gamma'] = 3e-3
cbl_settings_R['Hmax'] = 0.12
#cbl_settings_R['parameter_ensemble_min'] = np.array([0.5])
#cbl_settings_R['parameter_ensemble_max'] = np.array([2.5])
da_settings_R = dict(default_da_settings)
da_settings_R['cbl_settings'] = dict(cbl_settings_R)
da_settings_R['nens'] = nens
da_settings_R['type'] = 'real'
da_settings_R['observation_files'] = './LES_data/Theta/*csv'
da_settings_R['obs_coordinates'] = None
da_settings_R['obs_error_sdev_generate'] = None
da_settings_R['path'] = "./runs/real_init_exp_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
# Run and save to disk
try:
exp_R = pickle.load(open(da_settings_R['path']+'exp_R.pickle', "rb"))
except:
exp_R = experiment_real(da_settings_R)
setattr(exp_R,'label','R')
pickle.dump(exp_R, open(da_settings_R['path']+'exp_R.pickle', 'wb'))
############################################################################
# Run a free ensemble if desired
freens = False
if freens:
maxtime = 3600+21600 # spinup + analysis run
# Do a free ensemble run (ensemble size set expliclity)
cbl_settings_free = dict(cbl_settings_R)
cbl_settings_free['maxtime'] = maxtime
cbl_settings_free['perturb_ensemble_state'] = False
cbl_settings_free['parameter_ensemble_min'] = np.array([0.5])
cbl_settings_free['parameter_ensemble_max'] = np.array([2.5])
cbl_free = CBL(cbl_settings_free)
cbl_free.initialize(nens)
cbl_free.run(output_full_history=True)
# Make a plot about the free ensemble
fig, ax1 = p.subplots(1,1,constrained_layout=True)
fig.set_size_inches(4,3)
ax1,c1 = plot_spread(cbl_free,plot='mean',ax=ax1)
ax1.set_ylabel(r'Height (m)')
ax1.set_title(r'$\overline{\theta}$ (K)')
ax1.set_xlabel('Time (h)')
ax1.set_xticks(np.arange((maxtime+1)/3600))
c1.set_clim([290,296])
p.colorbar(c1,orientation='horizontal')
fig.savefig('expR_free_ens.png',format='png',dpi=300)
p.close(fig)
#########################################################################################à
# Make plots: perfect-model experiment (RC)
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
which_cycle = 0
tt1 = (which_cycle+1)*da_settings_RC['assimilation_interval']
ax1 = plot_CBL_assimilation(exp_RC,None,which_cycle=which_cycle,ax=ax1)
ax1.set_ylim([0,2000])
ax1.set_xlim([288,298])
ax1.set_ylabel(r'height (m)')
ax1.set_title(r'a) Analysis at $t=$%u s'%tt1)
ax1.legend(frameon=False)
which_cycle = 31
tt2 = (which_cycle+1)*da_settings_RC['assimilation_interval']
ax2 = plot_CBL_assimilation(exp_RC,None,which_cycle=which_cycle,ax=ax2)
ax2.set_ylim([0,2000])
ax2.set_xlim([288,298])
ax2.set_ylabel(r'height (m)')
ax2.set_title(r'b) Analysis at $t=$%u s'%tt2)
fig.savefig('expRC_profiles.png',format='png',dpi=300)
p.close(fig)
fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,6)
#
[ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_RC,da_settings_RC['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
ax0.set_title(r'a) Exp. RC, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
ax1.set_title(r'b) Exp. RC, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
ax2.set_title(r'c) Exp. RC, $\delta\overline{p}\prime\prime$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3 = plot_CBL_PE(exp_RC,None,ax=ax3)
ax3.set_title(r'd) Exp. RC, evolution of $p$')
ax3.set_xlabel('Time (h)')
ax3.set_yticks([0,1,2,3,4,5])
maxtime = exp_R.trun
ax0.set_xlim([0,maxtime/3600])
ax1.set_xlim([0,maxtime/3600])
ax2.set_xlim([0,maxtime/3600])
ax3.set_xlim([0,maxtime/3600])
ax0.set_xticks(np.arange((maxtime+1)/3600))
ax1.set_xticks(np.arange((maxtime+1)/3600))
ax2.set_xticks(np.arange((maxtime+1)/3600))
ax3.set_xticks(np.arange((maxtime+1)/3600))
p.colorbar(c0,orientation='horizontal')
p.colorbar(c1,orientation='horizontal')
p.colorbar(c2,orientation='horizontal')
#
fig.savefig('expRC_PE.png',format='png',dpi=300)
p.close(fig)
#########################################################################################à
# Make plots: perfect-model experiment (R)
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
which_cycle = 0
tt1 = (which_cycle+1)*da_settings_R['assimilation_interval']
ax1 = plot_CBL_assimilation(exp_R,None,which_cycle=which_cycle,ax=ax1)
ax1.set_ylim([0,2000])
ax1.set_xlim([288,298])
ax1.set_ylabel(r'height (m)')
ax1.set_title(r'a) Analysis at $t=$%u s'%tt1)
ax1.legend(frameon=False)
which_cycle = 31
tt2 = (which_cycle+1)*da_settings_R['assimilation_interval']
ax2 = plot_CBL_assimilation(exp_R,None,which_cycle=which_cycle,ax=ax2)
ax2.set_ylim([0,2000])
ax2.set_xlim([288,298])
ax2.set_ylabel(r'height (m)')
ax2.set_title(r'b) Analysis at $t=$%u s'%tt2)
fig.savefig('expR_profiles.png',format='png',dpi=300)
p.close(fig)
fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,6)
#
[ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_R,da_settings_R['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
ax0.set_title(r'a) Exp. R, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
ax1.set_title(r'b) Exp. R, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
ax2.set_title(r'c) Exp. R, $\delta\overline{p}\prime\prime$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3 = plot_CBL_PE(exp_R,None,ax=ax3)
ax3.set_title(r'd) Exp. R, evolution of $p$')
ax3.set_xlabel('Time (h)')
ax3.set_yticks([0,1,2,3,4,5])
maxtime = exp_R.trun
ax0.set_xlim([0,maxtime/3600])
ax1.set_xlim([0,maxtime/3600])
ax2.set_xlim([0,maxtime/3600])
ax3.set_xlim([0,maxtime/3600])
ax0.set_xticks(np.arange((maxtime+1)/3600))
ax1.set_xticks(np.arange((maxtime+1)/3600))
ax2.set_xticks(np.arange((maxtime+1)/3600))
ax3.set_xticks(np.arange((maxtime+1)/3600))
p.colorbar(c0,orientation='horizontal')
p.colorbar(c1,orientation='horizontal')
p.colorbar(c2,orientation='horizontal')
#
fig.savefig('expR_PE.png',format='png',dpi=300)
p.close(fig)
check_reproducibility = False
if check_reproducibility:
path = default_da_settings['path']
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment