diff --git a/config/cfg.py b/config/cfg.py
index e0ffa92219594d675e78ae4e65e319790960eb2e..5f3a21bb427c55104b492e7451cf8b0b8878649f 100755
--- a/config/cfg.py
+++ b/config/cfg.py
@@ -7,16 +7,17 @@ class ExperimentConfiguration(object):
pass
exp = ExperimentConfiguration()
-exp.expname = "exp_v1.11_LMU_filter_domainobs"
+exp.expname = "exp_v1.12_LMU_WV73_cde"
exp.model_dx = 2000
exp.timestep = 10
-exp.n_ens = 40
-exp.n_nodes = 10
-exp.n_obs = 100
-
-exp.sat_channels = [1,]
-exp.error_variance = 0.0009
+exp.n_ens = 20
+exp.n_nodes = 5
+exp.n_obs = 100 # radar: n_obs for each observation height level
+exp.sat_channels = [6,]
+exp.sat_err = 0.03
+exp.radar_err = 5.
+exp.distance_between_obs_meters = 10000
# directory paths depend on the name of the experiment
cluster.expname = exp.expname
diff --git a/config/clusters.py b/config/clusters.py
index 82af2b5bdf8ae5a623025d09e7c36caea9096973..56d0c9ac854ed8fa85dcafbcba708be75fbca7df 100755
--- a/config/clusters.py
+++ b/config/clusters.py
@@ -17,7 +17,7 @@ class ClusterConfig(object):
vsc = ClusterConfig()
vsc.name = 'vsc'
-vsc.python = '/home/fs71386/lkugler/miniconda3/bin/python'
+vsc.python = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/python'
vsc.ncks = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/ncks'
vsc.userdir = '/home/fs71386/lkugler'
vsc.srcdir = '/home/fs71386/lkugler/compile/WRF/WRF-4.2.1/run'
diff --git a/scripts/archive_assim.py b/scripts/archive_assim.py
index dbb88d90af66cb3d964a142f018091e01c9f0993..ad2c81bb2f2ec5f24dcc370a740d9df1c412b7ab 100755
--- a/scripts/archive_assim.py
+++ b/scripts/archive_assim.py
@@ -34,15 +34,17 @@ try:
for iens in range(1, exp.n_ens+1):
savedir = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/')+str(iens)
mkdir(savedir)
-
+
copy(cluster.dartrundir+'/input.nml',
cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/input.nml'))
filter_in = cluster.dartrundir+'/preassim_member_'+str(iens).zfill(4)+'.nc'
filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
- copy(filter_in, savedir+time.strftime('/%Y-%m-%d_%H:%M_prior'))
- copy(filter_out, savedir+time.strftime('/%Y-%m-%d_%H:%M_posterior'))
+ # files not necessary as wrfouts exist
+ #copy(filter_in, savedir+time.strftime('/%Y-%m-%d_%H:%M_prior'))
+ #copy(filter_out, savedir+time.strftime('/%Y-%m-%d_%H:%M_posterior'))
+
# copy mean and sd to archive
print('copy preassim, postassim mean and sd')
diff --git a/scripts/gen_synth_obs.py b/scripts/gen_synth_obs.py
new file mode 100755
index 0000000000000000000000000000000000000000..7a7b9ee3428aab9254a77c82fd5d4e4cdf053a7e
--- /dev/null
+++ b/scripts/gen_synth_obs.py
@@ -0,0 +1,171 @@
+import os, sys, shutil
+import datetime as dt
+import numpy as np
+from scipy.interpolate import interp1d
+from config.cfg import exp, cluster
+from utils import symlink, copy, sed_inplace, append_file
+import create_obsseq
+
+# fit of Fig 7, Harnisch 2016
+x_ci = [0, 5, 10.5, 13, 16]
+y_oe = [1, 4.5, 10, 12, 13] # Kelvin
+oe_73_linear = interp1d(x_ci, y_oe, assume_sorted=True)
+
+def oe_73(ci):
+ if ci < 13:
+ return oe_73_linear(ci)
+ else:
+ return 16.
+
+def cloudimpact_73(bt_mod, bt_obs):
+ """
+ follows Harnisch 2016
+ """
+ biascor_obs = 0
+ bt_lim = 255 # Kelvin for 7.3 micron WV channel
+
+ ci_obs = max(0, bt_lim-(bt_obs - biascor_obs))
+ ci_mod = max(0, bt_lim-bt_mod)
+ ci = (ci_obs+ci_mod)/2
+ return ci
+
+def read_prior_obs(f_obs_prior):
+ """
+ docstring
+ """
+ obsseq = open(f_obs_prior, 'r').readlines()
+ OBSs = []
+ # read observations from obs_seq.final
+ for i, line in enumerate(obsseq):
+ if ' OBS ' in line:
+ observed = float(obsseq[i+1].strip())
+ truth = float(obsseq[i+2].strip())
+ prior_ensmean = float(obsseq[i+3].strip())
+ prior_enssd = float(obsseq[i+4].strip())
+ prior_ens = []
+ for j in range(5, 5+exp.n_ens):
+ prior_ens.append(float(obsseq[i+j].strip()))
+
+ OBSs.append(dict(observed=observed, truth=truth,
+ prior_ens=np.array(prior_ens)))
+ return OBSs
+
+def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
+ """
+ overwrite observation errors in a obs_seq.out file
+ according to the values in OEs
+ """
+ # write to txt (write whole obs_seq.out again)
+ obsseq = open(fpath_obsseqin, 'r').readlines()
+ obsseq_new = obsseq.copy()
+ i_obs = 0
+ for i, line in enumerate(obsseq):
+ if 'kind\n' in line:
+ i_line_oe = i+9 # 9 for satellite obs
+ obsseq_new[i_line_oe] = ' '+str(OEs[i_obs])+' \n'
+ i_obs += 1
+
+ os.rename(fpath_obsseqin, fpath_obsseqin+'-bak') # backup
+ # write cloud dependent errors (actually whole file)
+ with open(fpath_obsseqin, 'w') as f:
+ for line in obsseq_new:
+ f.write(line)
+
+def set_input_nml(sat_channel=False, just_prior_values=False):
+ """descr"""
+ if just_prior_values:
+ template = cluster.scriptsdir+'/../templates/input.prioronly.nml'
+ else:
+ template = cluster.scriptsdir+'/../templates/input.nml'
+ copy(template, cluster.dartrundir+'/input.nml')
+ sed_inplace(cluster.dartrundir+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
+
+ # input.nml for RTTOV
+ if sat_channel > 0:
+ if sat_channel in [1, 2, 3, 12]:
+ rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.VIS.nml'
+ else:
+ rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
+ append_file(cluster.dartrundir+'/input.nml', rttov_nml)
+
+if __name__ == "__main__":
+
+ time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
+ print(dt.datetime.now())
+
+ # error_var = (5.)**2
+ # kind = 'RADAR'
+
+ # error_var = (0.5)**2
+ # kind = 'RASO_T'
+
+ #obs.generic_obs(kind, time, exp.n_obs, exp.radar_err**2, exp.distance_between_obs_meters,
+ # output_path=cluster.dartrundir,
+ # fpath_obs_locations=cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M')
+ # +'/obs_coords.pkl')
+
+
+ for channel_id in exp.sat_channels:
+
+ n_obs = 100
+ channel_id = 6
+ distance_between_obs_meters = 10000
+
+ error_var = (1.)**2
+ create_obsseq.sat(time, channel_id, n_obs, error_var,
+ distance_between_obs_meters,
+ output_path=cluster.dartrundir,
+ fpath_obs_locations='./domain.pkl')
+ print(dt.datetime.now())
+
+ if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
+ raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
+
+ # generate observations
+ set_input_nml(sat_channel=channel_id, just_prior_values=True)
+ os.chdir(cluster.dartrundir)
+ t = dt.datetime.now()
+ os.system('mpirun -np 12 ./perfect_model_obs')
+ print('1st perfect_model_obs', (dt.datetime.now()-t).total_seconds())
+
+ if channel_id == 6:
+ # run ./filter to have prior observation estimates from model state
+ set_input_nml(sat_channel=channel_id, just_prior_values=True)
+ t = dt.datetime.now()
+ os.system('mv obs_seq.out obs_seq_all.out; mpirun -np 20 ./filter')
+ print('1st filter', (dt.datetime.now()-t).total_seconds())
+
+ # find the observation error for a pair of (H(x_nature), H(x_background))
+ f_obs_prior = cluster.dartrundir+'/obs_seq.final'
+ OBSs = read_prior_obs(f_obs_prior)
+
+ # compute the observation error necessary
+ # to achieve a certain operational FGD distribution
+ OEs = []
+ for obs in OBSs:
+ bt_y = obs['truth']
+ bt_x_ens = obs['prior_ens']
+ CIs = [cloudimpact_73(bt_x, bt_y) for bt_x in bt_x_ens]
+
+ oe_nature = oe_73(np.mean(CIs))
+ OEs.append(oe_nature)
+
+ # write obs_seq.out
+ fpath_obsseqout = cluster.dartrundir+'/obs_seq.in'
+ edit_obserr_in_obsseq(fpath_obsseqout, OEs)
+ print('after editing oe', dt.datetime.now())
+
+ # generate actual observations (with correct error)
+ os.chdir(cluster.dartrundir)
+ t = dt.datetime.now()
+ os.system('mpirun -np 12 ./perfect_model_obs')
+ print('real obs gen', (dt.datetime.now()-t).total_seconds())
+
+ # correct input.nml for actual assimilation later on
+ set_input_nml(sat_channel=channel_id, just_prior_values=False)
+
+ # FIXME: missing is the concatenation of
+ os.chdir(cluster.dartrundir)
+ os.system('cat obs_seq.out >> obs_seq_all.out')
+ print(dt.datetime.now())
+ # FIXME: what if different channels in one obs_seq.out -> need different input.nml for different channels
diff --git a/scripts/pre_assim.py b/scripts/pre_assim.py
index 7f2667c9593360c0c6a087cba65cfae79a53d023..9438233acddaa48edc26829af0ede5770461607b 100755
--- a/scripts/pre_assim.py
+++ b/scripts/pre_assim.py
@@ -16,7 +16,7 @@ for iens in range(1, exp.n_ens+1):
print('link wrfout file to DART background file')
wrfout_run = exppath_firstguess+background_init_time.strftime('/%Y-%m-%d_%H:%M/') \
- +str(iens)+'/'+assim_time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
+ +str(iens)+assim_time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
dart_ensdir = cluster.dartrundir+'/advance_temp'+str(iens)
wrfout_dart = dart_ensdir+'/wrfout_d01'
@@ -49,5 +49,5 @@ os.system('rm -rf '+cluster.dartrundir+'/output_mean*')
os.system('rm -rf '+cluster.dartrundir+'/output_sd*')
os.system('rm -rf '+cluster.dartrundir+'/perfect_output_*')
-print('replace measurement error with obs error for assimilation') # FIXME !!! temporary only
-sed_inplace(cluster.dartrundir+'/obs_seq.out', '9.000000000000000E-004', '0.04')
+#print('replace measurement error with obs error for assimilation') # not needed anymore
+#sed_inplace(cluster.dartrundir+'/obs_seq.out', '9.000000000000000E-004', '0.04')
diff --git a/scripts/prepare_namelist.py b/scripts/prepare_namelist.py
index ea6fee41108862ef37dc6ed9b002be6f331a34f2..0fd483a67722e18b459b9734eea3b0ca3672b3ff 100755
--- a/scripts/prepare_namelist.py
+++ b/scripts/prepare_namelist.py
@@ -1,9 +1,9 @@
-import os, sys, shutil
+import os, sys, shutil, warnings
import datetime as dt
sys.path.append(os.getcwd())
from config.cfg import exp, cluster
-from utils import sed_inplace, copy, symlink
+from utils import sed_inplace, copy, symlink, mkdir
def run(cluster, iens, begin, end):
rundir = cluster.wrf_rundir(iens)
@@ -26,6 +26,17 @@ def run(cluster, iens, begin, end):
'<HH2>': '%H', '<MM2>': '%M'}.items():
sed_inplace(rundir+'/namelist.input', k, end.strftime(v))
+ #########################
+ try:
+ print('copy wrfinput of this run to archive')
+ wrfin_old = rundir+'/wrfinput_d01'
+ init_dir = cluster.archivedir()+begin.strftime('/%Y-%m-%d_%H:%M/')+str(iens)
+ os.makedirs(init_dir, exist_ok=True)
+ wrfin_arch = init_dir+'/wrfinput_d01'
+ copy(wrfin_old, wrfin_arch)
+ except Exception as e:
+ warnings.warn(str(e))
+
if __name__ == '__main__':
begin = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
diff --git a/scripts/run_obs_diag.py b/scripts/run_obs_diag.py
index cf8a78dbbb463ed713480739e97b157465573ba0..5fd07f90752f0a62f16e59eb23bbb25832663b43 100644
--- a/scripts/run_obs_diag.py
+++ b/scripts/run_obs_diag.py
@@ -28,19 +28,25 @@ def run(folder_obs_seq_final):
print('running obs_diag program')
os.chdir(rundir_program)
symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
- os.system('./obs_diag >& obs_diag.log')
+ try:
+ os.remove(rundir_program+'/obs_seq_to_netcdf')
+ except:
+ pass
+ os.system('./obs_diag >& obs_diag.log') # caution, this overwrites obs_seq_to_netcdf
outdir = '/'.join(folder_obs_seq_final.split('/')[:-1])
print('moving output to', outdir+'/obs_diag_output.nc')
copy(rundir_program+'/obs_diag_output.nc', outdir+'/obs_diag_output.nc')
print('running obs_seq_to_netcdf program')
+ shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', cluster.dart_srcdir+'/obs_seq_to_netcdf')
symlink(cluster.dart_srcdir+'/obs_seq_to_netcdf', rundir_program+'/obs_seq_to_netcdf')
- os.system('./obs_seq_to_netcdf >& obs_seq_to_netcdf.log')
- print('moving output to', outdir+'/obs_seq_output.nc')
- copy(rundir_program+'/obs_diag_output.nc', outdir+'/obs_seq_output.nc')
+ os.system('./obs_seq_to_netcdf >& obs_seq_to_netcdf.log') # caution, overwrites its own binary?!
+ print('moving output to', outdir+'/obs_seq...')
+ os.system('mv '+rundir_program+'/obs_epoch_*.nc '+outdir+'/')
if __name__ == '__main__':
- folder_obs_seq_final = '/home/fs71386/lkugler/data/sim_archive/exp_v1.11_LMU_filter_domainobs/obs_seq_final/'
+ #folder_obs_seq_final = '/home/fs71386/lkugler/data/sim_archive/exp_v1.11_LMU_filter2/obs_seq_final/'
+ folder_obs_seq_final = str(sys.argv[1])
run(folder_obs_seq_final)
diff --git a/scripts/update_wrfinput_from_filteroutput.py b/scripts/update_wrfinput_from_filteroutput.py
index fe9869b42cf21d21cc88907e2308eff9f5575924..9ccce8a32095ac2ab319396859a9d8ab5d8f17cf 100755
--- a/scripts/update_wrfinput_from_filteroutput.py
+++ b/scripts/update_wrfinput_from_filteroutput.py
@@ -1,27 +1,51 @@
import os, sys, warnings
import datetime as dt
+import netCDF4 as nc
from config.cfg import exp, cluster
from utils import symlink, copy_scp_srvx8, copy, mkdir, mkdir_srvx8, clean_wrfdir
time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
+background_init_time = dt.datetime.strptime(sys.argv[2], '%Y-%m-%d_%H:%M')
+exppath_firstguess = str(sys.argv[3])
#if cluster.name != 'srvx8':
# copy = copy_scp_srvx8
# mkdir = mkdir_srvx8
+cycle_vars = ['U', 'V', 'P', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW',
+ 'QGRAUP', 'QNICE', 'QNRAIN', 'U10', 'V10', 'T2', 'Q2', 'PSFC', 'TSLB',
+ 'SMOIS', 'TSK']
-update_vars = ['Times', 'U', 'V', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
+update_vars = ['Times', 'U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
+
+# variables which are updated need not to be cycled
+for var in update_vars:
+ if var in cycle_vars:
+ cycle_vars.remove(var)
+
+cycles = ','.join(cycle_vars)
+updates = ','.join(update_vars)
print('move output to WRF dir as new initial conditions')
for iens in range(1, exp.n_ens+1):
clean_wrfdir(cluster.wrf_rundir(iens))
+ prior_wrf = exppath_firstguess + background_init_time.strftime('/%Y-%m-%d_%H:%M/') \
+ +str(iens)+time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
wrf_ic = cluster.wrf_rundir(iens) + '/wrfinput_d01'
- # overwrite variables in wrfinput file
- vars = ','.join(update_vars)
- print('updating', vars, 'in', wrf_ic)
- os.system(cluster.ncks+' -A -v '+vars+' '+filter_out+' '+wrf_ic)
+ # cycles variables from wrfout (prior state)
+ print('cycling', cycles, 'into', wrf_ic, 'from', prior_wrf)
+ os.system(cluster.ncks+' -A -v '+cycles+' '+prior_wrf+' '+wrf_ic)
+
+ print('updating', updates, 'in', wrf_ic, 'from', filter_out)
+ os.system(cluster.ncks+' -A -v '+updates+' '+filter_out+' '+wrf_ic)
+
+ print('writing T into THM of wrfinput')
+ thm_in = nc.Dataset(filter_out, 'r').variables['T'][:]
+ dsout = nc.Dataset(wrf_ic, 'r+')
+ dsout.variables['THM'][:] = thm_in
+ dsout.close()
# clean up
#try:
diff --git a/templates/input.nml b/templates/input.nml
index 2910776ec314ba1100d9bf975015a60ceb18d855..101c3fdc534fd8ef5e8e2a9a8762cfe5e8193174 100644
--- a/templates/input.nml
+++ b/templates/input.nml
@@ -155,7 +155,8 @@
'ACARS_V_WIND_COMPONENT',
'ACARS_TEMPERATURE',
'GPSRO_REFRACTIVITY',
- 'LAND_SFC_ALTIMETER',
+ 'DOPPLER_RADIAL_VELOCITY',
+ 'RADAR_REFLECTIVITY',
'MSG_4_SEVIRI_RADIANCE',
'MSG_4_SEVIRI_TB',
'MSG_4_SEVIRI_BDRF'
@@ -177,7 +178,7 @@
lowest_reflectivity_fwd_op = -10.0,
max_radial_vel_obs = 1000000,
allow_wet_graupel = .false.,
- microphysics_type = 2 ,
+ microphysics_type = 5 ,
allow_dbztowt_conv = .false.,
dielectric_factor = 0.224,
n0_rain = 8.0e6,
@@ -218,6 +219,7 @@
'PSFC', 'QTY_SURFACE_PRESSURE', 'TYPE_PSFC', 'UPDATE','999',
'T2', 'QTY_2M_TEMPERATURE', 'TYPE_T', 'UPDATE','999',
'TSK', 'QTY_SKIN_TEMPERATURE', 'TYPE_T', 'UPDATE','999',
+ 'REFL_10CM','QTY_RADAR_REFLECTIVITY','TYPE_REFL', 'UPDATE','999',
wrf_state_bounds = 'QVAPOR','0.0','NULL','CLAMP',
'QRAIN', '0.0','NULL','CLAMP',
'QCLOUD','0.0','NULL','CLAMP',
@@ -298,10 +300,10 @@
&obs_diag_nml
obs_sequence_name = '',
obs_sequence_list = 'obsdiag_inputlist.txt',
- first_bin_center = 2008, 7,30,11, 0, 0 ,
- last_bin_center = 2008, 7,30,16, 0, 0 ,
- bin_separation = 0, 0, 0, 0,15, 0 ,
- bin_width = 0, 0, 0, 0,15, 0 ,
+ first_bin_center = 2008, 7,30,10, 0, 0 ,
+ last_bin_center = 2008, 7,30,18, 0, 0 ,
+ bin_separation = 0, 0, 0, 0,30, 0 ,
+ bin_width = 0, 0, 0, 0,30, 0 ,
time_to_skip = 0, 0, 0, 0, 0, 0 ,
max_num_bins = 1000,
trusted_obs = 'null',
@@ -314,7 +316,7 @@
print_mismatched_locs = .false.,
create_rank_histogram = .true.,
outliers_in_histogram = .true.,
- use_zero_error_obs = .false.,
+ use_zero_error_obs = .true.,
verbose = .false.
/
@@ -332,7 +334,7 @@
&obs_seq_to_netcdf_nml
obs_sequence_name = '',
obs_sequence_list = 'obsdiag_inputlist.txt',
- append_to_netcdf = .false.,
+ append_to_netcdf = .true.,
lonlim1 = 0.0,
lonlim2 = 360.0,
latlim1 = -90.0,
diff --git a/templates/input.prioronly.nml b/templates/input.prioronly.nml
new file mode 100644
index 0000000000000000000000000000000000000000..53ab2df3c9ec4df3e704497c3fba4a7e363c185a
--- /dev/null
+++ b/templates/input.prioronly.nml
@@ -0,0 +1,572 @@
+&perfect_model_obs_nml
+ read_input_state_from_file = .true.
+ single_file_in = .false.
+ input_state_files = 'wrfout_d01'
+ init_time_days = -1
+ init_time_seconds = -1
+
+ write_output_state_to_file = .false.
+ single_file_out = .false.
+ output_state_files = 'perfect_output_d01.nc'
+ output_interval = 1
+
+ obs_seq_in_file_name = "obs_seq.in"
+ obs_seq_out_file_name = "obs_seq.out"
+ first_obs_days = -1
+ first_obs_seconds = -1
+ last_obs_days = -1
+ last_obs_seconds = -1
+
+ async = 0
+ adv_ens_command = "../shell_scripts/advance_model.csh"
+
+ trace_execution = .true.
+ output_timestamps = .false.
+ print_every_nth_obs = -1
+ output_forward_op_errors = .false.
+ silence = .false.
+ /
+
+&filter_nml
+ async = 0,
+ adv_ens_command = "../shell_scripts/advance_model.csh",
+ ens_size = <n_ens>,
+ obs_sequence_in_name = "obs_seq_all.out",
+ obs_sequence_out_name = "obs_seq.final",
+ input_state_file_list = "input_list.txt"
+ output_state_file_list = "output_list.txt"
+ init_time_days = -1,
+ init_time_seconds = -1,
+ first_obs_days = -1,
+ first_obs_seconds = -1,
+ last_obs_days = -1,
+ last_obs_seconds = -1,
+ num_output_state_members = <n_ens>,
+ num_output_obs_members = <n_ens>,
+ output_interval = 1,
+ num_groups = 1,
+ distributed_state = .true.
+ compute_posterior = .false.
+ output_forward_op_errors = .false.,
+ output_timestamps = .false.,
+ trace_execution = .false.,
+
+ stages_to_write = 'preassim'
+ output_members = .false.
+ output_mean = .false.
+ output_sd = .false.
+ write_all_stages_at_end = .false.
+
+ inf_flavor = 0, 0,
+ inf_initial_from_restart = .true., .false.,
+ inf_sd_initial_from_restart = .true., .false.,
+ inf_initial = 1.0, 1.00,
+ inf_lower_bound = 1.0, 1.0,
+ inf_upper_bound = 1000000.0, 1000000.0,
+ inf_damping = 0.9, 1.0,
+ inf_sd_initial = 0.6, 0.0,
+ inf_sd_lower_bound = 0.6, 0.0,
+ inf_sd_max_change = 1.05, 1.05,
+ /
+
+&quality_control_nml
+ input_qc_threshold = 3.0,
+ outlier_threshold = 3.0,
+ enable_special_outlier_code = .false.
+ /
+
+&fill_inflation_restart_nml
+ write_prior_inf = .false.
+ prior_inf_mean = 1.00
+ prior_inf_sd = 0.6
+
+ write_post_inf = .false.
+ post_inf_mean = 1.00
+ post_inf_sd = 0.6
+
+ input_state_files = 'wrfinput_d01', 'wrfinput_d02'
+ single_file = .false.
+ verbose = .false.
+ /
+
+&smoother_nml
+ num_lags = 0,
+ start_from_restart = .false.,
+ output_restart = .false.,
+ restart_in_file_name = 'smoother_ics',
+ restart_out_file_name = 'smoother_restart',
+ /
+
+# cutoff is in radians; for the earth, 0.05 is about 300 km.
+# cutoff is defined to be the half-width of the localization radius,
+# so 0.05 radians for cutoff is about an 600 km effective
+# localization radius, where the influence of an obs decreases
+# to ~half at 300 km, and ~0 at the edges of the area.
+&assim_tools_nml
+ filter_kind = 1,
+ cutoff = 0.0001,
+ sort_obs_inc = .false.,
+ spread_restoration = .false.,
+ sampling_error_correction = .false.,
+ adaptive_localization_threshold = -1,
+ output_localization_diagnostics = .false.,
+ localization_diagnostics_file = 'localization_diagnostics',
+ convert_all_state_verticals_first = .true.
+ convert_all_obs_verticals_first = .true.
+ print_every_nth_obs = 0,
+ /
+
+&cov_cutoff_nml
+ select_localization = 1,
+ /
+
+&obs_sequence_nml
+ write_binary_obs_sequence = .false.,
+ /
+
+&preprocess_nml
+ overwrite_output = .true.,
+ input_obs_kind_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90',
+ output_obs_kind_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90',
+ input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90',
+ output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90',
+ input_files = '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90',
+ '../../../observations/forward_operators/obs_def_radar_mod.f90',
+ '../../../observations/forward_operators/obs_def_metar_mod.f90',
+ '../../../observations/forward_operators/obs_def_dew_point_mod.f90',
+ '../../../observations/forward_operators/obs_def_rel_humidity_mod.f90',
+ '../../../observations/forward_operators/obs_def_altimeter_mod.f90',
+ '../../../observations/forward_operators/obs_def_gps_mod.f90',
+ '../../../observations/forward_operators/obs_def_vortex_mod.f90',
+ '../../../observations/forward_operators/obs_def_gts_mod.f90',
+ '../../../observations/forward_operators/obs_def_rttov_mod.f90',
+ /
+
+&obs_kind_nml
+ assimilate_these_obs_types = 'RADIOSONDE_TEMPERATURE',
+ 'RADIOSONDE_U_WIND_COMPONENT',
+ 'RADIOSONDE_V_WIND_COMPONENT',
+ 'SAT_U_WIND_COMPONENT',
+ 'SAT_V_WIND_COMPONENT',
+ 'AIRCRAFT_U_WIND_COMPONENT',
+ 'AIRCRAFT_V_WIND_COMPONENT',
+ 'AIRCRAFT_TEMPERATURE',
+ 'ACARS_U_WIND_COMPONENT',
+ 'ACARS_V_WIND_COMPONENT',
+ 'ACARS_TEMPERATURE',
+ 'GPSRO_REFRACTIVITY',
+ 'DOPPLER_RADIAL_VELOCITY',
+ 'RADAR_REFLECTIVITY',
+ 'MSG_4_SEVIRI_RADIANCE',
+ 'MSG_4_SEVIRI_TB',
+ 'MSG_4_SEVIRI_BDRF'
+ evaluate_these_obs_types = 'RADIOSONDE_SPECIFIC_HUMIDITY',
+ /
+
+# Notes for obs_def_radar_mod_nml:
+# (1) Reflectivity limit can be applied to observations and/or forward operator.
+# (2) The default constants below match the WRF defaults. They will need to
+# be changed for other cases, depending on which microphysics scheme is used.
+#
+
+&obs_def_radar_mod_nml
+ apply_ref_limit_to_obs = .false.,
+ reflectivity_limit_obs = -10.0,
+ lowest_reflectivity_obs = -10.0,
+ apply_ref_limit_to_fwd_op = .false.,
+ reflectivity_limit_fwd_op = -10.0,
+ lowest_reflectivity_fwd_op = -10.0,
+ max_radial_vel_obs = 1000000,
+ allow_wet_graupel = .false.,
+ microphysics_type = 5 ,
+ allow_dbztowt_conv = .false.,
+ dielectric_factor = 0.224,
+ n0_rain = 8.0e6,
+ n0_graupel = 4.0e6,
+ n0_snow = 3.0e6,
+ rho_rain = 1000.0,
+ rho_graupel = 400.0,
+ rho_snow = 100.0,
+ /
+
+
+# Notes for model_nml:
+# (1) vert_localization_coord must be one of:
+# 1 = model level
+# 2 = pressure
+# 3 = height
+# 4 = scale height
+
+# set default_state_variables to .false. to use the explicit list.
+# otherwise it uses a hardcoded default list: U, V, W, PH, T, MU, QV only.
+# see ../wrf_state_variables_table for a full list of what wrf fields are
+# supported in the DART state vector, and what settings should be used here.
+# 'UPDATE' and 'NO_COPY_BACK' are supported in the 4th column; 'NO_UPDATE' is
+# not yet supported.
+
+&model_nml
+ default_state_variables = .false.,
+ wrf_state_variables = 'U', 'QTY_U_WIND_COMPONENT', 'TYPE_U', 'UPDATE','999',
+ 'V', 'QTY_V_WIND_COMPONENT', 'TYPE_V', 'UPDATE','999',
+ 'W', 'QTY_VERTICAL_VELOCITY', 'TYPE_W', 'UPDATE','999',
+ 'PH', 'QTY_GEOPOTENTIAL_HEIGHT', 'TYPE_GZ', 'UPDATE','999',
+ 'T', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
+ 'MU', 'QTY_PRESSURE', 'TYPE_MU', 'UPDATE','999',
+ 'QVAPOR','QTY_VAPOR_MIXING_RATIO', 'TYPE_QV', 'UPDATE','999',
+ 'QICE', 'QTY_ICE_MIXING_RATIO', 'TYPE_QI', 'UPDATE','999',
+ 'QCLOUD','QTY_CLOUDWATER_MIXING_RATIO','TYPE_QC', 'UPDATE','999',
+ 'CLDFRA','QTY_CLOUD_FRACTION', 'TYPE_CFRAC','UPDATE','999',
+ 'PSFC', 'QTY_SURFACE_PRESSURE', 'TYPE_PSFC', 'UPDATE','999',
+ 'T2', 'QTY_2M_TEMPERATURE', 'TYPE_T', 'UPDATE','999',
+ 'TSK', 'QTY_SKIN_TEMPERATURE', 'TYPE_T', 'UPDATE','999',
+ 'REFL_10CM','QTY_RADAR_REFLECTIVITY','TYPE_REFL', 'UPDATE','999',
+ wrf_state_bounds = 'QVAPOR','0.0','NULL','CLAMP',
+ 'QRAIN', '0.0','NULL','CLAMP',
+ 'QCLOUD','0.0','NULL','CLAMP',
+ 'QICE','0.0','NULL','CLAMP',
+ 'CLDFRA','0.0','1.0','CLAMP',
+ num_domains = 1,
+ calendar_type = 3,
+ assimilation_period_seconds = 21600,
+ vert_localization_coord = 3,
+ center_search_half_length = 500000.,
+ center_spline_grid_scale = 10,
+ sfc_elev_max_diff = -1.0,
+ circulation_pres_level = 80000.0,
+ circulation_radius = 108000.0,
+ allow_obs_below_vol = .false.
+ /
+
+# vert_normalization_X is amount of X equiv to 1 radian in horiz.
+# vert localization is 'cutoff' times the pressure/height/levels,
+# only if horiz_dist_only is set to .false. in the namelist below.
+# the default nlon/nlat should be good for most experiments. it sets
+# an internal structure that speeds up searches. don't change it
+# based on your grid size. nlon must be an odd number.
+&location_nml
+ horiz_dist_only = .false.,
+ vert_normalization_pressure = 6666666.7,
+ vert_normalization_height = 5000000.0,
+ vert_normalization_level = 2666.7,
+ vert_normalization_scale_height = 10.0,
+ approximate_distance = .false.,
+ nlon = 71,
+ nlat = 36,
+ output_box_info = .false.,
+ /
+
+&utilities_nml
+ TERMLEVEL = 1,
+ module_details = .false.,
+ logfilename = 'dart_log.out',
+ nmlfilename = 'dart_log.nml',
+ write_nml = 'file',
+ /
+
+&mpi_utilities_nml
+ /
+
+®_factor_nml
+ select_regression = 1,
+ input_reg_file = "time_mean_reg",
+ save_reg_diagnostics = .true.,
+ reg_diagnostics_file = "reg_diagnostics",
+ /
+
+# layout = 2 spreads the IO tasks across the nodes.
+# This can greatly improve the performance in IO if
+# tasks_per_node is set to match your hardware
+&ensemble_manager_nml
+ layout = 2,
+ tasks_per_node = 48
+ /
+
+&obs_def_gps_nml
+ max_gpsro_obs = 100000,
+ /
+
+&obs_def_tpw_nml
+ /
+
+# The times in the namelist for the obs_diag program are vectors
+# that follow the following sequence:
+# year month day hour minute second
+# max_num_bins can be used to specify a fixed number of bins,
+# in which case last_bin_center should be safely in the future.
+#
+# Acceptable latitudes range from [-90, 90]
+# Acceptable longitudes range from [ 0, Inf]
+
+&obs_diag_nml
+ obs_sequence_name = '',
+ obs_sequence_list = 'obsdiag_inputlist.txt',
+ first_bin_center = 2008, 7,30,10, 0, 0 ,
+ last_bin_center = 2008, 7,30,18, 0, 0 ,
+ bin_separation = 0, 0, 0, 0,30, 0 ,
+ bin_width = 0, 0, 0, 0,30, 0 ,
+ time_to_skip = 0, 0, 0, 0, 0, 0 ,
+ max_num_bins = 1000,
+ trusted_obs = 'null',
+ Nregions = 1,
+ lonlim1 = 0.0, 0.0, 0.0, 235.0,
+ lonlim2 = 360.0, 360.0, 360.0, 295.0,
+ latlim1 = 20.0, -80.0, -20.0, 25.0,
+ latlim2 = 80.0, -20.0, 20.0, 55.0,
+ reg_names = 'all',
+ print_mismatched_locs = .false.,
+ create_rank_histogram = .true.,
+ outliers_in_histogram = .true.,
+ use_zero_error_obs = .true.,
+ verbose = .false.
+ /
+
+&schedule_nml
+ calendar = 'Gregorian',
+ first_bin_start = 1601, 1, 1, 0, 0, 0,
+ first_bin_end = 2999, 1, 1, 0, 0, 0,
+ last_bin_end = 2999, 1, 1, 0, 0, 0,
+ bin_interval_days = 1000000,
+ bin_interval_seconds = 0,
+ max_num_bins = 1000,
+ print_table = .true.,
+ /
+
+&obs_seq_to_netcdf_nml
+ obs_sequence_name = '',
+ obs_sequence_list = 'obsdiag_inputlist.txt',
+ append_to_netcdf = .true.,
+ lonlim1 = 0.0,
+ lonlim2 = 360.0,
+ latlim1 = -90.0,
+ latlim2 = 90.0,
+ verbose = .false.,
+ /
+
+# There is one GIGANTIC difference between the obsdef_mask.txt and .nc
+# The netCDF file intentionally ignores the effect of nTmin/nTmax.
+# The netCDF file has ALL matching stations, regardless of temporal coverage.
+
+&obs_seq_coverage_nml
+ obs_sequences = ''
+ obs_sequence_list = 'obs_coverage_list.txt'
+ obs_of_interest = 'METAR_U_10_METER_WIND'
+ textfile_out = 'METAR_U_10_METER_WIND_obsdef_mask.txt'
+ netcdf_out = 'METAR_U_10_METER_WIND_obsdef_mask.nc'
+ first_analysis = 2003, 1, 1, 0, 0, 0
+ last_analysis = 2003, 1, 2, 0, 0, 0
+ forecast_length_days = 1
+ forecast_length_seconds = 0
+ verification_interval_seconds = 21600
+ temporal_coverage_percent = 100.0
+ lonlim1 = 0.0
+ lonlim2 = 360.0
+ latlim1 = -90.0
+ latlim2 = 90.0
+ verbose = .true.
+ /
+
+# selections_file is a list of obs_defs output
+# from the obs_seq_coverage utility.
+
+&obs_selection_nml
+ filename_seq = 'obs_seq.out',
+ filename_seq_list = '',
+ filename_out = 'obs_seq.processed',
+ selections_file = 'obsdef_mask.txt',
+ selections_is_obs_seq = .false.,
+ print_only = .false.,
+ calendar = 'gregorian',
+ /
+
+&obs_seq_verify_nml
+ obs_sequences = ''
+ obs_sequence_list = 'obs_verify_list.txt'
+ input_template = 'obsdef_mask.nc'
+ netcdf_out = 'forecast.nc'
+ obtype_string = 'METAR_U_10_METER_WIND'
+ print_every = 10000
+ verbose = .true.
+ debug = .false.
+ /
+
+&obs_sequence_tool_nml
+ num_input_files = 1,
+ filename_seq = 'obs_seq.out',
+ filename_out = 'obs_seq.processed',
+ first_obs_days = -1,
+ first_obs_seconds = -1,
+ last_obs_days = -1,
+ last_obs_seconds = -1,
+ obs_types = '',
+ keep_types = .false.,
+ print_only = .false.,
+ min_lat = -90.0,
+ max_lat = 90.0,
+ min_lon = 0.0,
+ max_lon = 360.0,
+ /
+
+&replace_wrf_fields_nml
+ debug = .false.,
+ fail_on_missing_field = .false.,
+ fieldnames = "SNOWC",
+ "ALBBCK",
+ "TMN",
+ "TSK",
+ "SH2O",
+ "SMOIS",
+ "SEAICE",
+ "HGT_d01",
+ "TSLB",
+ "SST",
+ "SNOWH",
+ "SNOW",
+ fieldlist_file = '',
+ /
+
+&obs_common_subset_nml
+ num_to_compare_at_once = 2,
+ filename_seq = 'obs_seq1.final', 'obs_seq2.final',
+ filename_seq_list = '',
+ filename_out_suffix = '.common' ,
+ calendar = 'Gregorian',
+ print_every = 1000,
+ dart_qc_threshold = 3,
+ print_only = .false.,
+ /
+
+&wrf_dart_to_fields_nml
+ include_slp = .true.,
+ include_wind_components = .true.,
+ include_height_on_pres = .true.,
+ include_temperature = .true.,
+ include_rel_humidity = .true.,
+ include_surface_fields = .false.,
+ include_sat_ir_temp = .false.,
+ pres_levels = 70000.,
+ /
+
+&ncepobs_nml
+ year = 2010,
+ month = 06,
+ day = 00,
+ tot_days = 1,
+ max_num = 1000000,
+ ObsBase = 'temp_obs.',
+ select_obs = 0,
+ ADPUPA = .false.,
+ AIRCAR = .false.,
+ AIRCFT = .false.,
+ SATEMP = .false.,
+ SFCSHP = .false.,
+ ADPSFC = .false.,
+ SATWND = .true.,
+ obs_U = .false.,
+ obs_V = .false.,
+ obs_T = .false.,
+ obs_PS = .false.,
+ obs_QV = .false.,
+ daily_file = .true.,
+ obs_time = .false.,
+ lat1 = 10.00,
+ lat2 = 60.00,
+ lon1 = 210.0,
+ lon2 = 300.0
+ /
+
+&prep_bufr_nml
+ obs_window_upa = 1.0,
+ obs_window_air = 1.0,
+ obs_window_cw = 1.0,
+ otype_use = 242.0, 243.0, 245.0, 246.0, 251.0, 252.0, 253.0, 257.0, 259.0
+ qctype_use = 0, 1, 2, 3, 4, 9, 15
+ /
+
+&convert_cosmic_gps_nml
+ gpsro_netcdf_file = '',
+ gpsro_netcdf_filelist = 'flist',
+ gpsro_out_file = 'obs_seq.gpsro',
+ local_operator = .true.,
+ obs_levels = 0.22, 0.55, 1.1, 1.8, 2.7, 3.7, 4.9,
+ 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0,
+ ray_ds = 5000.0,
+ ray_htop = 13000.1,
+ /
+
+&wrf_obs_preproc_nml
+
+ file_name_input = 'obs_seq20110901'
+ file_name_output = 'obs_seq.europe.prev'
+
+ overwrite_obs_time = .false.
+
+ obs_boundary = 0.0
+ increase_bdy_error = .false.
+ maxobsfac = 2.5
+ obsdistbdy = 1.0
+
+ sfc_elevation_check = .false.
+ sfc_elevation_tol = 3000.0
+ obs_pressure_top = 0.0
+ obs_height_top = 2.0e10
+
+ include_sig_data = .true.
+ tc_sonde_radii = -1.0
+
+ superob_aircraft = .true.
+ aircraft_horiz_int = 800.0
+ aircraft_pres_int = 25000.0
+
+ superob_sat_winds = .true.
+ sat_wind_horiz_int = 800.0
+ sat_wind_pres_int = 25000.0
+
+ overwrite_ncep_satwnd_qc = .false.
+ overwrite_ncep_sfc_qc = .false.
+ /
+
+! sonde_extra = 'obs_seq.rawin'
+! land_sfc_extra = 'obs_seq.land_sfc'
+! metar_extra = 'obs_seq.metar'
+! marine_sfc_extra = 'obs_seq.marine'
+! sat_wind_extra = 'obs_seq.satwnd'
+! profiler_extra = 'obs_seq.profiler'
+! gpsro_extra = 'obs_seq.gpsro'
+! acars_extra = 'obs_seq.acars'
+! trop_cyclone_extra = 'obs_seq.tc'
+
+&state_vector_io_nml
+ single_precision_output = .true.,
+ /
+
+&compare_states_nml
+ /
+
+&closest_member_tool_nml
+ input_restart_file_list = 'input_file_list_d01.txt',
+ output_file_name = 'closest_results.txt'
+ ens_size = 3,
+ single_restart_file_in = .false.,
+ difference_method = 4,
+ use_only_qtys = 'QTY_U_WIND_COMPONENT'
+ /
+
+# To test both domains, you must change 'model_nml:num_domains = 2'
+
+&model_mod_check_nml
+ input_state_files = 'wrfinput_d01', 'wrfinput_d02'
+ output_state_files = 'mmc_output1.nc', 'mmc_output2.nc'
+ test1thru = 0
+ run_tests = 1,2,3,4,5
+ x_ind = 87370
+ loc_of_interest = 231.0, 40.0, 10.0
+ quantity_of_interest = 'QTY_U_WIND_COMPONENT'
+ interp_test_dlon = 0.1
+ interp_test_dlat = 0.1
+ interp_test_dvert = 1000.0
+ interp_test_lonrange = 250.0, 260.0
+ interp_test_latrange = 30.0, 45.0
+ interp_test_vertrange = 2000.0, 4000.0
+ interp_test_vertcoord = 'VERTISHEIGHT'
+ verbose = .false.
+ /
diff --git a/templates/namelist.input b/templates/namelist.input
index bfd3b8a6e887523b2cd74973652835e8a8738a27..74bc39c72d3f07d6281e73c44bfcfeb60ecd09eb 100644
--- a/templates/namelist.input
+++ b/templates/namelist.input
@@ -90,6 +90,7 @@
tke_adv_opt = 1, 1, 1,
non_hydrostatic = .true., .true., .true.,
mix_full_fields = .true., .true., .true.,
+ use_theta_m = 0,
/
&bdy_control