diff --git a/ENDA.py b/ENDA.py
index 3fd5c2343d2604b388d9b7fcc8678fc780385d5e..58d2aab3146375e2223937d8f64c5508d4a6297f 100644
--- a/ENDA.py
+++ b/ENDA.py
@@ -281,7 +281,7 @@ class cycle:
localization_cutoff,\
inflate_rtps,\
inflation_rtps_alpha,\
- RNG):
+ rnseed):
# How many cycles do you need to run?
ncycles = int(trun/assimilation_interval)
@@ -515,8 +515,13 @@ class experiment:
elif self.obs_kinds[j]=='v':
variable = nr.history['v'][:,time_index]
truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate)
+ # Change the seed every time you go through this,
+ # otherwise the perturbations are always the same;
+ # and you have bias instead of randomness
+ seed = self.rnseed+j*1000+i*100000
+ RNG = np.random.default_rng(seed=seed)
observations[i,j] = truths[i,j]+\
- self.RNG.normal(0,self.obs_error_sdev_generate[j])
+ RNG.normal(0,self.obs_error_sdev_generate[j])
# Store truths and observations
self.truths = truths
@@ -540,7 +545,7 @@ class experiment:
self.localization_cutoff,
self.inflate_rtps,
self.inflation_rtps_alpha,
- self.RNG)
+ self.rnseed)
# Compute diagnostics
if not hasattr(self,'dg'):
diff --git a/PE_CBL.py b/PE_CBL.py
index 941c268f332cb4e18d7882d5cbf9bb655b21e8e6..9cf0c58006b5fd242a0159d162954c2fb402b8d4 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -3,6 +3,7 @@
import numpy as np
from matplotlib import pyplot as p
import pickle
+import os
from copy import deepcopy
# Own modules
@@ -19,8 +20,8 @@ if __name__ == '__main__':
sigma_o_gn = 0.0 # set equal to sigma_o_as in experiments C and F
sigma_o_as = 0.1 # increased by a factor 10 in experiments C and F
- # Initiate a random number generator with preset seed, to ensure reproducibility
- RNG = np.random.default_rng(seed=181612)
+ # Define seed for random number generation, to ensure reproducibility
+ rnseed = 181612
# Default CBL model settings
default_cbl_settings ={
@@ -47,12 +48,11 @@ if __name__ == '__main__':
# Domain parameters (Kmax determines dt)
'dz' : 50,
'ztop' : 4000,
- 'maxtime' : 10800,
# 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
- 'RNG' : RNG,
+ 'rnseed' : rnseed,
'perturb_ensemble_state' : True,
'perturbations_type' : "smooth",
'perturbations_theta_amplitude' : sigma_b_init,
@@ -78,7 +78,7 @@ if __name__ == '__main__':
# Default assimilation settings
default_da_settings = {'cbl_settings' : default_cbl_settings,
'tspinup' : 3600,
- 'trun' : 10800,
+ 'trun' : 21600,
'assimilation_interval' : 300,
'obs_coordinates' : np.linspace(0,2000,nobs),
'obs_kinds' : np.array(['theta']*nobs),
@@ -88,7 +88,8 @@ if __name__ == '__main__':
'nens' : nens,
'FILTER' : 'EAKF',
'inflation_rtps_alpha' : 0.2,
- 'RNG' : RNG
+ 'rnseed' : rnseed,
+ 'path' : "./runs/seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
}
# Safety check on time steps
@@ -99,6 +100,10 @@ if __name__ == '__main__':
# Applies only to sets B&C; no-PE run is always computed for experiments A&D
noPE_runs = False
+ # Check if output path exists, if not, create it.
+ if not os.path.exists(default_da_settings['path']):
+ os.makedirs(default_da_settings['path'])
+
############################################################################
# Experiment A (control)
# Create a copy of the default settings
@@ -107,11 +112,11 @@ if __name__ == '__main__':
# Run it
try:
- exp_A = pickle.load(open("exp_A.pickle", "rb"))
+ exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
except:
exp_A = experiment(da_settings_A)
setattr(exp_A,'label','A')
- pickle.dump(exp_A, open('exp_A.pickle', 'wb'))
+ pickle.dump(exp_A, open(default_da_settings['path']+'exp_A.pickle', 'wb'))
# Store the nature run, because it is used also in other experiments/figures
exp_A.nr.has_parameters_in_state = True
@@ -127,14 +132,15 @@ if __name__ == '__main__':
da_settings_A_noPE['observations'] = observations
da_settings_A_noPE['nr'].do_parameter_estimation = False
cbl_settings_A_noPE['do_parameter_estimation'] = False
+ cbl_settings_A_noPE['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
da_settings_A_noPE['cbl_settings'] = cbl_settings_A_noPE
# Run it
try:
- exp_A_noPE = pickle.load(open("exp_A_noPE.pickle", "rb"))
+ exp_A_noPE = pickle.load(open(default_da_settings['path']+'exp_A_noPE.pickle', "rb"))
except:
exp_A_noPE = experiment(da_settings_A_noPE)
- pickle.dump(exp_A_noPE, open('exp_A_noPE.pickle', 'wb'))
+ pickle.dump(exp_A_noPE, open(default_da_settings['path']+'exp_A_noPE.pickle', 'wb'))
########################################################################
# Experiment B1
@@ -154,11 +160,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B1 = pickle.load(open("exp_B1.pickle", "rb"))
+ exp_B1 = pickle.load(open(default_da_settings['path']+'exp_B1.pickle', "rb"))
except:
exp_B1 = experiment(da_settings_B1)
setattr(exp_B1,'label','B1')
- pickle.dump(exp_B1, open('exp_B1.pickle', 'wb'))
+ pickle.dump(exp_B1, open(default_da_settings['path']+'exp_B1.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -171,10 +177,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B1_noPE = pickle.load(open("exp_B1_noPE.pickle", "rb"))
+ exp_B1_noPE = pickle.load(open(default_da_settings['path']+'exp_B1_noPE.pickle', "rb"))
except:
exp_B1_noPE = experiment(da_settings_B1_noPE)
- pickle.dump(exp_B1_noPE, open('exp_B1_noPE.pickle', 'wb'))
+ pickle.dump(exp_B1_noPE, open(default_da_settings['path']+'exp_B1_noPE.pickle', 'wb'))
########################################################################
# Experiment B2
@@ -194,11 +200,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B2 = pickle.load(open("exp_B2.pickle", "rb"))
+ exp_B2 = pickle.load(open(default_da_settings['path']+'exp_B2.pickle', "rb"))
except:
exp_B2 = experiment(da_settings_B2)
setattr(exp_B2,'label','B2')
- pickle.dump(exp_B2, open('exp_B2.pickle', 'wb'))
+ pickle.dump(exp_B2, open(default_da_settings['path']+'exp_B2.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -211,10 +217,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B2_noPE = pickle.load(open("exp_B2_noPE.pickle", "rb"))
+ exp_B2_noPE = pickle.load(open(default_da_settings['path']+'exp_B2_noPE.pickle', "rb"))
except:
exp_B2_noPE = experiment(da_settings_B2_noPE)
- pickle.dump(exp_B2_noPE, open('exp_B2_noPE.pickle', 'wb'))
+ pickle.dump(exp_B2_noPE, open(default_da_settings['path']+'exp_B2_noPE.pickle', 'wb'))
########################################################################
# Experiment B3
@@ -234,11 +240,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B3 = pickle.load(open("exp_B3.pickle", "rb"))
+ exp_B3 = pickle.load(open(default_da_settings['path']+'exp_B3.pickle', "rb"))
except:
exp_B3 = experiment(da_settings_B3)
setattr(exp_B3,'label','B3')
- pickle.dump(exp_B3, open('exp_B3.pickle', 'wb'))
+ pickle.dump(exp_B3, open(default_da_settings['path']+'exp_B3.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -251,10 +257,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B3_noPE = pickle.load(open("exp_B3_noPE.pickle", "rb"))
+ exp_B3_noPE = pickle.load(open(default_da_settings['path']+'exp_B3_noPE.pickle', "rb"))
except:
exp_B3_noPE = experiment(da_settings_B3_noPE)
- pickle.dump(exp_B3_noPE, open('exp_B3_noPE.pickle', 'wb'))
+ pickle.dump(exp_B3_noPE, open(default_da_settings['path']+'exp_B3_noPE.pickle', 'wb'))
########################################################################
# Experiment B4
@@ -272,11 +278,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B4 = pickle.load(open("exp_B4.pickle", "rb"))
+ exp_B4 = pickle.load(open(default_da_settings['path']+'exp_B4.pickle', "rb"))
except:
exp_B4 = experiment(da_settings_B4)
setattr(exp_B4,'label','B4')
- pickle.dump(exp_B4, open('exp_B4.pickle', 'wb'))
+ pickle.dump(exp_B4, open(default_da_settings['path']+'exp_B4.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -287,10 +293,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_B4_noPE = pickle.load(open("exp_B4_noPE.pickle", "rb"))
+ exp_B4_noPE = pickle.load(open(default_da_settings['path']+'exp_B4_noPE.pickle', "rb"))
except:
exp_B4_noPE = experiment(da_settings_B4_noPE)
- pickle.dump(exp_B4_noPE, open('exp_B4_noPE.pickle', 'wb'))
+ pickle.dump(exp_B4_noPE, open(default_da_settings['path']+'exp_B4_noPE.pickle', 'wb'))
########################################################################
# Experiment C1
@@ -310,11 +316,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
+ exp_C1 = pickle.load(open(default_da_settings['path']+'exp_C1.pickle', "rb"))
except:
exp_C1 = experiment(da_settings_C1)
setattr(exp_C1,'label','C1')
- pickle.dump(exp_C1, open('exp_C1.pickle', 'wb'))
+ pickle.dump(exp_C1, open(default_da_settings['path']+'exp_C1.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -325,10 +331,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_C1_noPE = pickle.load(open("exp_C1_noPE.pickle", "rb"))
+ exp_C1_noPE = pickle.load(open(default_da_settings['path']+'exp_C1_noPE.pickle', "rb"))
except:
exp_C1_noPE = experiment(da_settings_C1_noPE)
- pickle.dump(exp_C1_noPE, open('exp_C1_noPE.pickle', 'wb'))
+ pickle.dump(exp_C1_noPE, open(default_da_settings['path']+'exp_C1_noPE.pickle', 'wb'))
########################################################################
# Experiment C2
@@ -348,11 +354,11 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
+ exp_C2 = pickle.load(open(default_da_settings['path']+'exp_C2.pickle', "rb"))
except:
exp_C2 = experiment(da_settings_C2)
setattr(exp_C2,'label','C2')
- pickle.dump(exp_C2, open('exp_C2.pickle', 'wb'))
+ pickle.dump(exp_C2, open(default_da_settings['path']+'exp_C2.pickle', 'wb'))
if noPE_runs:
# Corresponding experiment without parameter estimation
@@ -363,10 +369,10 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_C2_noPE = pickle.load(open("exp_C2_noPE.pickle", "rb"))
+ exp_C2_noPE = pickle.load(open(default_da_settings['path']+'exp_C2_noPE.pickle', "rb"))
except:
exp_C2_noPE = experiment(da_settings_C2_noPE)
- pickle.dump(exp_C2_noPE, open('exp_C2_noPE.pickle', 'wb'))
+ pickle.dump(exp_C2_noPE, open(default_da_settings['path']+'exp_C2_noPE.pickle', 'wb'))
########################################################################
# Experiment D
@@ -389,23 +395,24 @@ if __name__ == '__main__':
# Run and save to disk
try:
- exp_D = pickle.load(open("exp_D.pickle", "rb"))
+ exp_D = pickle.load(open(default_da_settings['path']+'exp_D.pickle', "rb"))
except:
exp_D = experiment(da_settings_D)
setattr(exp_D,'label','D')
- pickle.dump(exp_D, open('exp_D.pickle', 'wb'))
+ pickle.dump(exp_D, open(default_da_settings['path']+'exp_D.pickle', 'wb'))
# Experiment matching D, but without parameter estimation
- cbl_settings_D_noPE = dict(cbl_settings_D)
- da_settings_D_noPE = dict(da_settings_D)
- cbl_settings_D_noPE['do_parameter_estimation'] = False
- da_settings_D_noPE['cbl_settings'] = cbl_settings_D_noPE
+ if noPE_runs:
+ cbl_settings_D_noPE = dict(cbl_settings_D)
+ da_settings_D_noPE = dict(da_settings_D)
+ cbl_settings_D_noPE['do_parameter_estimation'] = False
+ da_settings_D_noPE['cbl_settings'] = cbl_settings_D_noPE
- try:
- exp_D_noPE = pickle.load(open("exp_D_noPE.pickle", "rb"))
- except:
- exp_D_noPE = experiment(da_settings_D_noPE)
- pickle.dump(exp_A_noPE, open('exp_D_noPE.pickle', 'wb'))
+ try:
+ exp_D_noPE = pickle.load(open(default_da_settings['path']+'exp_D_noPE.pickle', "rb"))
+ except:
+ exp_D_noPE = experiment(da_settings_D_noPE)
+ pickle.dump(exp_A_noPE, open(default_da_settings['path']+'exp_D_noPE.pickle', 'wb'))
########################################################################
@@ -417,7 +424,8 @@ if __name__ == '__main__':
fig05 = False
fig06 = False
fig07 = False
- fig08 = True
+ fig08 = False
+ fig09 = True
# Other possible, optional, plots
opt01 = False # assimilation of a single observation
@@ -427,6 +435,8 @@ if __name__ == '__main__':
# Create a copy of the default settings
cbl_settings = dict(default_cbl_settings)
+ maxtime = 21600
+ cbl_settings['maxtime'] = maxtime
# Disable parameter estimation (not used here)
cbl_settings['do_parameter_estimation'] = False
@@ -454,6 +464,7 @@ if __name__ == '__main__':
# Panel c) spread induced by p
# Do a free ensemble run (ensemble size set expliclity)
cbl_settings_free = dict(default_cbl_settings)
+ cbl_settings_free['maxtime'] = maxtime
cbl_settings_free['perturb_ensemble_state'] = False
cbl_free = CBL(cbl_settings_free)
cbl_free.initialize(nens)
@@ -468,7 +479,7 @@ if __name__ == '__main__':
ax1.set_ylim([0,zmax])
ax1.set_ylabel(r'Height (m)')
ax1.set_xlabel(r'Time (h)')
- ax1.set_xticks(np.arange(4))
+ ax1.set_xticks(np.arange((maxtime+1)/3600))
ax1.set_title(r'c) $\overline{\theta}$ (K)')
p.colorbar(c1,orientation='horizontal')
ax1.contour(cbl_det.history['time']/3600,cbl_det.zt,cbl_det.history['theta'],
@@ -489,7 +500,7 @@ if __name__ == '__main__':
ax3.set_ylabel(r'Height (m)')
ax3.set_title(r'd) $\sigma_\theta$ (K)')
ax3.set_xlabel('Time (h)')
- ax3.set_xticks(np.arange(4))
+ ax3.set_xticks(np.arange((maxtime+1)/3600))
p.colorbar(c3,orientation='horizontal')
zoverh = np.linspace(0,1,101)
@@ -509,25 +520,34 @@ if __name__ == '__main__':
if fig02:
- exp_A = pickle.load(open("exp_A.pickle", "rb"))
+ exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
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_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. A, $\rho(p\prime\prime,y_b}$)')
+ ax0.set_title(r'a) Exp. A$_1$, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. A, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
+ ax1.set_title(r'b) Exp. A$_1$, $\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. A, $\delta\overline{p}\prime\prime$')
+ ax2.set_title(r'c) Exp. A$_1$, $\delta\overline{p}\prime\prime$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3 = plot_CBL_PE(exp_A,None,ax=ax3)
- ax3.set_title(r'd) Exp. A, evolution of $p$')
+ ax3.set_title(r'd) Exp. A$_1$, evolution of $p$')
ax3.set_xlabel('Time (h)')
ax3.set_yticks([0,1,2,3,4,5])
+ maxtime = exp_A.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')
@@ -537,8 +557,8 @@ if __name__ == '__main__':
if fig03:
- exp_A = pickle.load(open("exp_A.pickle", "rb"))
- exp_A_noPE = pickle.load(open("exp_A_noPE.pickle", "rb"))
+ 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"))
experiments_pe = [exp_A]
experiments_nope = [exp_A_noPE]
@@ -548,10 +568,10 @@ if __name__ == '__main__':
if fig04:
- exp_B1 = pickle.load(open("exp_B1.pickle", "rb"))
- exp_B2 = pickle.load(open("exp_B2.pickle", "rb"))
- exp_B3 = pickle.load(open("exp_B3.pickle", "rb"))
- exp_B4 = pickle.load(open("exp_B4.pickle", "rb"))
+ 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"))
+ exp_B3 = pickle.load(open(default_da_settings['path']+"exp_B3.pickle", "rb"))
+ exp_B4 = pickle.load(open(default_da_settings['path']+"exp_B4.pickle", "rb"))
fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,4)
@@ -586,8 +606,8 @@ if __name__ == '__main__':
if fig05:
- exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
- exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
+ 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"))
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,2)
@@ -610,8 +630,8 @@ if __name__ == '__main__':
if fig06:
- exp_C1 = pickle.load(open("exp_C1.pickle", "rb"))
- exp_C2 = pickle.load(open("exp_C2.pickle", "rb"))
+ 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"))
npar2 = exp_C1.nr.parameter_number
npar3 = exp_C2.nr.parameter_number
@@ -631,6 +651,11 @@ if __name__ == '__main__':
ax2.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
ax2.sharey(ax1)
+ maxtime = exp_C1.trun
+ ax1.set_xlim([0,maxtime/3600])
+ ax2.set_xlim([0,maxtime/3600])
+ ax1.set_xticks(np.arange((maxtime+1)/3600))
+ ax2.set_xticks(np.arange((maxtime+1)/3600))
p.setp(ax2.get_yticklabels(), visible=False)
#
p.colorbar(c1,orientation='horizontal')
@@ -641,7 +666,7 @@ if __name__ == '__main__':
if fig07:
- exp_D = pickle.load(open("exp_D.pickle", "rb"))
+ exp_D = pickle.load(open(default_da_settings['path']+"exp_D.pickle", "rb"))
fig, [ax0, ax1] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
@@ -654,6 +679,11 @@ if __name__ == '__main__':
ax1.set_title(r'b) Exp. D, evolution of $p$')
ax1.set_xlabel('Time (h)')
ax1.set_yticks([0,1,2,3,4,5])
+ maxtime = exp_D.trun
+ ax0.set_xlim([0,maxtime/3600])
+ ax1.set_xlim([0,maxtime/3600])
+ ax0.set_xticks(np.arange((maxtime+1)/3600))
+ ax1.set_xticks(np.arange((maxtime+1)/3600))
p.colorbar(c0,orientation='horizontal')
#
fig.savefig('fig07.png',format='png',dpi=300)
@@ -661,19 +691,24 @@ if __name__ == '__main__':
if fig08:
- exp_A = pickle.load(open("exp_A.pickle", "rb"))
- exp_D = pickle.load(open("exp_D.pickle", "rb"))
+ 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"))
fig, [[ax0, ax1],[ax2,ax3]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,6)
#
[ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax1,ax0])
- ax0.set_title(r'a) Exp. A, $\delta\overline{y}$ (K)')
+ ax0.set_title(r'a) Exp. A$_1$, $\delta\overline{y}$ (K)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. A, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
+ ax1.set_title(r'b) Exp. A$_1$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
+ maxtime = exp_A.trun
+ ax0.set_xlim([0,maxtime/3600])
+ ax1.set_xlim([0,maxtime/3600])
+ ax0.set_xticks(np.arange((maxtime+1)/3600))
+ ax1.set_xticks(np.arange((maxtime+1)/3600))
[ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],ax=[ax3,ax2])
ax2.set_title(r'c) Exp. D, $\delta\overline{y}$ (K)')
ax2.set_xlabel('Time (h)')
@@ -681,6 +716,11 @@ if __name__ == '__main__':
ax3.set_title(r'd) Exp. D, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
ax3.set_xlabel('Time (h)')
ax3.set_ylabel('Height (m)')
+ maxtime = exp_D.trun
+ ax2.set_xlim([0,maxtime/3600])
+ ax3.set_xlim([0,maxtime/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')
@@ -689,6 +729,13 @@ if __name__ == '__main__':
fig.savefig('fig08.png',format='png',dpi=300)
p.close(fig)
+ if fig09:
+
+ path_EAKF = "./runs/seed=181612_enssize=200_EAKF_6hrs"
+ path_LETKF = "./runs/seed=181612_enssize=200_LETKF_6hrs"
+
+
+
if opt01:
da_settings = {'cbl_settings' : dict(default_cbl_settings),
@@ -732,3 +779,22 @@ if __name__ == '__main__':
fig.savefig('profiles.png',format='png',dpi=300)
p.close(fig)
+
+ check_reproducibility = False
+ if check_reproducibility:
+ path = default_da_settings['path']
+ for i in ["exp_A.pickle",
+ "exp_A_noPE.pickle",
+ "exp_B1.pickle",
+ "exp_B2.pickle",
+ "exp_B3.pickle",
+ "exp_B4.pickle",
+ "exp_C1.pickle",
+ "exp_C2.pickle",
+ "exp_D.pickle"]:
+ experiment = pickle.load(open(path+i, "rb"))
+ print(i)
+ print("OBS :",experiment.da.observations[:,0])
+ print("PARAMETERS:",experiment.da.initial_perturbed_parameters[0,:5])
+ print("STATE :",experiment.da.initial_state[0,0,:5])
+
diff --git a/graphics.py b/graphics.py
index 56f6fda759478d5027a20dd8563057fb7154fa7e..7afdc5bbfaeec983087e8c180fd628773336a35d 100644
--- a/graphics.py
+++ b/graphics.py
@@ -421,7 +421,7 @@ def plot_diagnostics(experiments_pe,experiments_nope,labels,filename):
fig.set_size_inches(6,4)
z = experiments_pe[0].obs_coordinates
z_pbl = z*1.
- z_pbl[z>1000] = np.nan
+ z_pbl[z>1500] = np.nan
for i in range(len(experiments_pe)):
i1 = experiments_pe[i].dg
i2 = experiments_nope[i].dg
diff --git a/models.py b/models.py
index 8e2915a63b892858e2860a35cca2f9271bea5463..62c584d910e7403c826be3f85dc28cb8466bce74 100644
--- a/models.py
+++ b/models.py
@@ -97,7 +97,8 @@ class CBL:
# Overwrite the pre-existing unperturbed values
# This is needed for safety, because parameter transformations may be nonlinear
if self.nens>1:
- if self.do_parameter_estimation and hasattr(self,"initial_perturbed_parameters"):
+ #if self.do_parameter_estimation and hasattr(self,"initial_perturbed_parameters"):
+ if hasattr(self,"initial_perturbed_parameters"):
x0[-self.parameter_number:,:] = self.initial_perturbed_parameters
else:
if self.perturb_ensemble_parameters or self.do_parameter_estimation:
@@ -158,7 +159,8 @@ class CBL:
pp = np.zeros((self.parameter_number,self.nens))
for k in range(-self.parameter_number,0):
- dum = self.RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
+ RNG = np.random.default_rng(seed=self.rnseed)
+ dum = RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
pp[k,:] = self.parameter_transform[k](dum , kind='dir')
return pp
@@ -172,10 +174,12 @@ class CBL:
randomsize = 1
if self.perturbations_type == "random" or self.perturbations_type == "uniform":
- ppt = self.RNG.standard_normal((randomsize,self.nens))
+ RNG = np.random.default_rng(seed=self.rnseed)
+ ppt = RNG.standard_normal((randomsize,self.nens))
if self.is_bwind:
- ppu = self.RNG.standard_normal((self.nens,randomsize))
- ppv = self.RNG.standard_normal((self.nens,randomsize))
+ RNG = np.random.default_rng(seed=self.rnseed)
+ ppu = RNG.standard_normal((self.nens,randomsize))
+ ppv = RNG.standard_normal((self.nens,randomsize))
# Smooth perturbations are slightly more complicated
if self.perturbations_type == "smooth":
@@ -188,15 +192,17 @@ class CBL:
# Draw random perturbations, then interpolate
randomsize = self.perturbations_smooth_number
ppt = np.zeros((self.nz,self.nens))+np.nan
- pert_t = self.RNG.standard_normal((randomsize,self.nens))
+ RNG = np.random.default_rng(seed=self.rnseed)
+ pert_t = RNG.standard_normal((randomsize,self.nens))
for n in range(self.nens):
f = CubicSpline(ipert,pert_t[:,n])
ppt[:,n] = f(np.arange(self.nz))
if self.is_bwind:
ppu = np.zeros((self.nz,self.nens))+np.nan
ppv = np.zeros((self.nz,self.nens))+np.nan
- pert_u = self.RNG.standard_normal((randomsize,self.nens))
- pert_v = self.RNG.standard_normal((randomsize,self.nens))
+ RNG = np.random.default_rng(seed=self.rnseed)
+ pert_u = RNG.standard_normal((randomsize,self.nens))
+ pert_v = RNG.standard_normal((randomsize,self.nens))
for n in range(self.nens):
f = CubicSpline(ipert,pert_u[:,n])
ppu[:,n] = f(np.arange(self.nz))
@@ -311,7 +317,11 @@ class CBL:
H0 = Hmax
# Add random perturbations to the initial value
- H0 += self.RNG.normal(scale=H0_perturbation_ampl_init)
+ # Change the seed every time you go through this,
+ # otherwise the perturbations are always the same
+ seed = self.rnseed+j*200000
+ RNG = np.random.default_rng(seed=seed)
+ H0 += RNG.normal(scale=H0_perturbation_ampl_init)
# Set the surface momentum flux (ustar)
ustar = 0
@@ -360,8 +370,13 @@ class CBL:
rdz = 1./self.dz
# Add time-dependent surface perturbations
+ # Change the seed every time you go through this,
+ # otherwise the perturbations are always the same
+ seed = self.rnseed+j*200000
+ RNG = np.random.default_rng(seed=seed)
+
# Then compute sensible heat flux and integrate T equation
- H[0,j] = H0 + self.RNG.normal(scale=H0_perturbation_ampl_time)
+ H[0,j] = H0 + RNG.normal(scale=H0_perturbation_ampl_time)
H[1:-1,j] = -K[1:-1,j]*( (thetap[1:]-thetap[:-1])*rdz - gammac)
H[nz,j] = 2*H[nz-1,j]-H[nz-2,j]
theta[:,j] = thetap[:]-dt*rdz*(H[1:,j]-H[:-1,j])