diff --git a/config/cfg.py b/config/cfg.py
index 005b1e14fa39d900fb3210ae9db605ee41e29cc0..0636712223129221510872e7ce30b960dcd371e9 100755
--- a/config/cfg.py
+++ b/config/cfg.py
@@ -9,22 +9,27 @@ class ExperimentConfiguration(object):
exp = ExperimentConfiguration()
-exp.expname = "exp_v1.16_Pwbub-1_40mem"
+exp.expname = "exp_v1.17_Pwbub-1_WV_obs20_loc5"
exp.model_dx = 2000
exp.n_ens = 40
exp.n_nodes = 10
-n_obs = 121 #961 900: 10km resoltn # radar: n_obs for each observation height level
+# localize vertically, if it has a vertical position
+# needs a horizontal scale too, to calculate the vertical normalization
+# since you can not specify different vertical localizations for diff. variables
+exp.cov_loc_vert_km_horiz_km = (3, 5)
+
+n_obs = 256 #121: 30km, 256:16x16 (20km); 961: 10km resoltn # radar: n_obs for each observation height level
vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs,
error_generate=0.03, error_assimilate=0.06,
- cov_loc_radius_km=32)
+ cov_loc_radius_km=7.5)
wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs,
- error_generate=1., error_assimilate=False,
- cov_loc_radius_km=32)
+ error_generate=1., error_assimilate=False, #
+ cov_loc_radius_km=5)
ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs,
@@ -34,21 +39,21 @@ ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
kind='RADAR_REFLECTIVITY', n_obs=n_obs,
error_generate=2.5, error_assimilate=5.,
- heights=np.arange(1000, 15001, 1000),
- cov_loc_radius_km=32, cov_loc_vert_km=4)
+ heights=[7000,], #np.arange(1000, 15001, 1000),
+ cov_loc_radius_km=5)
t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
kind='SYNOP_TEMPERATURE', n_obs=n_obs,
error_generate=0.1, error_assimilate=1.,
- cov_loc_radius_km=20, cov_loc_vert_km=3)
+ cov_loc_radius_km=20)
psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs,
error_generate=50., error_assimilate=100.,
- cov_loc_radius_km=32, cov_loc_vert_km=5)
+ cov_loc_radius_km=32)
-exp.observations = [] #wv73, vis] # 108, wv73, vis]
+exp.observations = [radar] # 108, wv73, vis]
exp.update_vars = ['T', 'QVAPOR', 'QCLOUD', 'QICE','CLDFRA']
#exp.update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'TSK', 'CLDFRA']
diff --git a/scripts/apply_obs_op_dart.py b/scripts/apply_obs_op_dart.py
index 8db2c2e40afadd10e809cd619f076eb5f65cf4c6..6939ad6e1ccfb4758f5659b0e543fcce722e197b 100755
--- a/scripts/apply_obs_op_dart.py
+++ b/scripts/apply_obs_op_dart.py
@@ -48,37 +48,34 @@ if __name__ == '__main__':
# only the prior state values are of interest in this file
# observation and truth is wrong in this file (dummy)
- istage = 0
- archive_stage = savedir+'/assim_stage'+str(istage)
- aso.archive_diagnostics(archive_stage, time)
+ aso.archive_osq_final(time, posterior_1min=True)
- sys.exit() # multi stage below
+ # sys.exit() # multi stage below
+ # n_stages = len(exp.observations)
+ # for istage, obscfg in enumerate(exp.observations):
- n_stages = len(exp.observations)
- for istage, obscfg in enumerate(exp.observations):
+ # n_obs = obscfg['n_obs']
+ # sat_channel = obscfg.get('sat_channel', False)
+ # obscfg['folder_obs_coords'] = False
- n_obs = obscfg['n_obs']
- sat_channel = obscfg.get('sat_channel', False)
- obscfg['folder_obs_coords'] = False
+ # aso.set_DART_nml(cov_loc_radius_km=obscfg['cov_loc_radius_km'],
+ # cov_loc_vert_km=obscfg.get('cov_loc_vert_km', False),
+ # just_prior_values=True)
- aso.set_DART_nml(cov_loc_radius_km=obscfg['cov_loc_radius_km'],
- cov_loc_vert_km=obscfg.get('cov_loc_vert_km', False),
- just_prior_values=True)
+ # osq.create_obsseq_in(time, obscfg)
- osq.create_obsseq_in(time, obscfg)
+ # # prepare dummy nature (this Hx is irrelevant)
+ # os.chdir(cluster.dartrundir)
+ # os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
+ # wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
+ # cluster.dartrundir+'/wrfout_d01')
+ # aso.run_perfect_model_obs()
+ # aso.assimilate(nproc=96)
- # prepare dummy nature (this Hx is irrelevant)
- os.chdir(cluster.dartrundir)
- os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
- wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
- cluster.dartrundir+'/wrfout_d01')
- aso.run_perfect_model_obs()
- aso.assimilate(nproc=96)
-
- # only the prior state values are of interest in this file
- # observation and truth is wrong in this file (dummy)
- archive_stage = savedir+'/assim_stage'+str(istage)
- aso.archive_diagnostics(archive_stage, time)
+ # # only the prior state values are of interest in this file
+ # # observation and truth is wrong in this file (dummy)
+ # archive_stage = savedir+'/assim_stage'+str(istage)
+ # aso.archive_diagnostics(archive_stage, time)
diff --git a/scripts/assim_synth_obs.py b/scripts/assim_synth_obs.py
index f80ddccc17c747f4ab3f53149498a0a513bc4921..e39ad55526388ecde8d7e634984f1fd37bcc65eb 100755
--- a/scripts/assim_synth_obs.py
+++ b/scripts/assim_synth_obs.py
@@ -97,7 +97,7 @@ def replace_errors_obsseqout(f, new_errors):
previous_err_var = obsseq[line_error_obs_i]
new_err_obs_i = new_errors[i_obs]**2 # variance in obs_seq.out
- if debug: print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
+ if debug: print('previous err var ', float(previous_err_var.strip()), 'new error', new_err_obs_i)
obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n'
i_obs += 1 # next iteration
@@ -110,6 +110,7 @@ def replace_errors_obsseqout(f, new_errors):
def set_DART_nml_singleobstype(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
just_prior_values=False):
+ """this function is outdated and will not work"""
cov_loc_radian = cov_loc_radius_km/earth_radius_km
if just_prior_values:
@@ -145,8 +146,19 @@ def set_DART_nml_singleobstype(sat_channel=False, cov_loc_radius_km=32, cov_loc_
rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
append_file(cluster.dartrundir+'/input.nml', rttov_nml)
-def set_DART_nml(cov_loc_radius_km=32, cov_loc_vert_km=False, just_prior_values=False):
- cov_loc_radian = cov_loc_radius_km/earth_radius_km
+def set_DART_nml(just_prior_values=False):
+
+ def to_radian_horizontal(cov_loc_horiz_km):
+ cov_loc_radian = cov_loc_horiz_km/earth_radius_km
+ return cov_loc_radian
+
+ def to_vertical_normalization(cov_loc_vert_km, cov_loc_horiz_km):
+ vert_norm_rad = earth_radius_km*cov_loc_vert_km/cov_loc_horiz_km*1000
+ return vert_norm_rad
+
+ list_obstypes = [obscfg['kind'] for obscfg in exp.observations]
+ list_cov_loc_radius_km = [obscfg['cov_loc_radius_km'] for obscfg in exp.observations]
+ list_cov_loc_radian = [str(to_radian_horizontal(a)) for a in list_cov_loc_radius_km]
if just_prior_values:
template = cluster.scriptsdir+'/../templates/input.eval.nml'
@@ -156,15 +168,20 @@ def set_DART_nml(cov_loc_radius_km=32, cov_loc_vert_km=False, just_prior_values=
# options keys are replaced in input.nml with values
options = {'<n_ens>': str(int(exp.n_ens)),
- '<cov_loc_radian>': str(cov_loc_radian)}
+ '<cov_loc_radian>': '0.00000001', # dummy value, used for types not mentioned below
+ '<list_obstypes>': "'" + "','".join(list_obstypes) + "'",
+ '<list_cutoffs>': ", ".join(list_cov_loc_radian),
+ }
- if cov_loc_vert_km:
- vert_norm_rad = earth_radius_km*cov_loc_vert_km/cov_loc_radius_km*1000
- options['<horiz_dist_only>'] = '.false.'
- options['<vert_norm_hgt>'] = str(vert_norm_rad)
- else:
- options['<horiz_dist_only>'] = '.true.'
- options['<vert_norm_hgt>'] = '50000.0' # dummy value
+ # if cov_loc_vert_km:
+ options['<horiz_dist_only>'] = '.true.'
+
+ cov_loc_vert_km, cov_loc_horiz_km = exp.cov_loc_vert_km_horiz_km
+ vert_norm_hgt = to_vertical_normalization(cov_loc_vert_km, cov_loc_horiz_km)
+ options['<vert_norm_hgt>'] = str(vert_norm_hgt)
+ # else:
+ # options['<horiz_dist_only>'] = '.true.'
+ # options['<vert_norm_hgt>'] = '50000.0' # dummy value
for key, value in options.items():
sed_inplace(cluster.dartrundir+'/input.nml', key, value)
@@ -295,42 +312,56 @@ def recycle_output():
############### archiving
-def archive_assimilation(time)
+def archive_osq_final(time, posterior_1min=False):
+ """Save obs_seq.final file for later.
+ time (dt.datetime) : time of sampling values from files
+ posterior_1min (bool) : False if usual assimilation
+ """
- print('archive obs space diagnostics')
- archive_dir = cluster.archivedir+'/obs_seq_final/'
+ if posterior_1min:
+ archive_dir = cluster.archivedir+'/obs_seq_final_1min/'
+ else:
+ archive_dir = cluster.archivedir+'/obs_seq_final/'
mkdir(archive_dir)
fout = archive_dir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final')
copy(cluster.dartrundir+'/obs_seq.final', fout)
print(fout, 'saved.')
- # try: # what are regression diagnostics?!
- # print('archive regression diagnostics')
+ # try:
+ # print('archive regression diagnostics') # what are regression diagnostics?!
# copy(cluster.dartrundir+'/reg_diagnostics', archive_dir+'/reg_diagnostics')
# except Exception as e:
# warnings.warn(str(e))
+def archive_filteroutput(time):
print('archiving output')
- archive_assim = cluster.archivedir + '/assim_stage0/'
+ archive_assim = cluster.archivedir + time.strftime('/%Y-%m-%d_%H:%M/assim_stage0/')
mkdir(archive_assim)
copy(cluster.dartrundir+'/input.nml', archive_assim+'/input.nml')
for iens in range(1, exp.n_ens+1): # single members
- filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
- copy(filter_out, archive_assim+'/filter_restart_d01.'+str(iens).zfill(4))
+ copy(cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4),
+ archive_assim+'/filter_restart_d01.'+str(iens).zfill(4))
- for f in ['output_mean.nc', 'output_sd.nc']: # copy mean and sd to archive
- copy(cluster.dartrundir+'/'+f, archive_assim+'/'+f)
+ try: # not necessary for next forecast run
+ for iens in range(1, exp.n_ens+1):
+ copy(cluster.dartrundir+'/postassim_member_'+str(iens).zfill(4)+'.nc',
+ archive_assim+'/postassim_member_'+str(iens).zfill(4)+'.nc')
-def archive_obs_generation(time):
+ for f in ['output_mean.nc', 'output_sd.nc']: # copy mean and sd to archive
+ copy(cluster.dartrundir+'/'+f, archive_assim+'/'+f)
+ except Exception as e:
+ warnings.warn(str(e))
+
+
+def archive_osq_out(time):
dir_obsseq = cluster.archivedir+'/obs_seq_out/'
os.makedirs(dir_obsseq, exist_ok=True)
copy(cluster.dartrundir+'/obs_seq.out', dir_obsseq+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.out'))
if __name__ == "__main__":
-
"""Assimilate observations (different obs types)
as defined in config/cfg.py
for a certain timestamp (argument) of the nature run (defined in config/clusters.py)
@@ -339,7 +370,7 @@ if __name__ == "__main__":
for each assimilation stage (one obs_seq.in and e.g. one observation type):
1) create obs_seq.in with obs-errors
2) prepare nature run for DART
- 3) create obs from nature (obs_seq.out)
+ 3) create obs from nature (obs_seq.out) using pre-defined obs-errors
4) Assimilate
- adapt obs errors for assimilation
- calculate assim obs error from parametrization
@@ -395,7 +426,7 @@ if __name__ == "__main__":
Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
err_assim = calc_obserr_WV73(Hx_nat, Hx_prior)
-
+
error_assimilate.extend(err_assim) # the obs-error we assume for assimilating observations
################################################
@@ -408,12 +439,10 @@ if __name__ == "__main__":
osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=error_generate)
- first_obstype = exp.observations[0] # TODO: different for each observation type
- set_DART_nml(cov_loc_radius_km=first_obstype['cov_loc_radius_km'],
- cov_loc_vert_km=first_obstype.get('cov_loc_vert_km', False))
+ set_DART_nml()
run_perfect_model_obs() # actually create observations that are used to assimilate
- archive_obs_generation(time)
+ archive_osq_out(time)
################################################
print(' 3) assimilate with observation-errors for assimilation')
@@ -423,4 +452,5 @@ if __name__ == "__main__":
assimilate()
print('filter took', time_module.time()-t, 'seconds')
- archive_assimilation(time)
+ archive_filteroutput(time)
+ archive_osq_final(time)
diff --git a/templates/input.nml b/templates/input.nml
index 0c0c1cee45f7d555bf6df602f35a43a122021e02..356d7258e25220b20d55fb62530cb22e2690bacc 100644
--- a/templates/input.nml
+++ b/templates/input.nml
@@ -51,7 +51,7 @@
output_timestamps = .false.,
trace_execution = .false.,
- stages_to_write = 'output'
+ stages_to_write = 'postassim', 'output'
output_members = .true.
output_mean = .true.
output_sd = .true.
@@ -114,6 +114,8 @@
convert_all_state_verticals_first = .true.
convert_all_obs_verticals_first = .true.
print_every_nth_obs = 0,
+ special_localization_obs_types = <list_obstypes>,
+ special_localization_cutoffs = <list_cutoffs>
/
&cov_cutoff_nml