diff --git a/PE_CBL.py b/PE_CBL.py
index 6b50f42731a3f070fa309053e81e69e9dbf156c1..56e9a99fd483d33f2836f66823b88e7c20862034 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -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']