From 7bfd329fd25f0e42adea7430a78a8ce5703d969e Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Wed, 3 May 2023 15:40:48 +0200
Subject: [PATCH] consistent configurations
---
config/{cfg.py => exp_example.py} | 74 +++++-----
cycled_exp.py | 12 +-
dartwrf/assim_synth_obs.py | 70 +---------
dartwrf/create_obsseq.py | 2 -
dartwrf/dart_nml.py | 216 ++++++++++++++++++++++++++++++
5 files changed, 268 insertions(+), 106 deletions(-)
rename config/{cfg.py => exp_example.py} (53%)
create mode 100644 dartwrf/dart_nml.py
diff --git a/config/cfg.py b/config/exp_example.py
similarity index 53%
rename from config/cfg.py
rename to config/exp_example.py
index 8a833a6..2f158be 100755
--- a/config/cfg.py
+++ b/config/exp_example.py
@@ -1,86 +1,98 @@
from dartwrf import utils
exp = utils.Experiment()
-exp.expname = "obs-cold_localization-narrow"
+exp.expname = "test_newcode"
exp.model_dx = 2000
-exp.n_ens = 10
-
-exp.filter_kind = 1
-exp.prior_inflation = 0
-exp.post_inflation = 4
-exp.sec = True
-exp.cov_loc_vert_km_horiz_km = (4, 40)
+exp.n_ens = 40
exp.superob_km = False # False or int (spatial averaging of observations)
-exp.adjust_obs_impact = False
exp.use_existing_obsseq = False # False or pathname (use precomputed obs_seq.out files)
-exp.use_existing_obsseq = '/users/students/lehre/advDA_s2023/dartwrf_tutorial/very_cold_observation.out'
+#exp.use_existing_obsseq = '/users/students/lehre/advDA_s2023/dartwrf_tutorial/very_cold_observation.out'
# path to the nature run, where we take observations from
-exp.nature = '/users/students/lehre/advDA_s2023/data/sample_nature/'
+exp.nature_wrfout = '/jetfs/home/lkugler/data/sim_archive/exp_v1.18_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
exp.input_profile = '/mnt/jetfs/home/lkugler/data/initial_profiles/wrf/ens/2022-03-31/raso.fc.<iens>.wrfprof'
+
+exp.dart_nml = {'&assim_tools_nml':
+ dict(assim_tools_nml='.false.',
+ filter_kind='1',
+ sampling_error_correction='.true.',
+ # obs_impact_filename='/jetfs/home/lkugler/DART-WRF/templates/impactfactor_T.txt',
+ ),
+ '&filter_nml':
+ dict(ens_size=str(exp.n_ens),
+ num_output_state_members=str(exp.n_ens),
+ num_output_obs_members=str(exp.n_ens),
+ inf_flavor=['2', '0'],
+ ),
+ '&location_nml':
+ dict(horiz_dist_only='.true.',
+ ),
+ }
+
+
# n_obs can be 22500: 2km, 5776: 4km, 121: 30km, 256:16x16 (20km); 961: 10km resoltn
# if radar: then n_obs is for each observation height level
-vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
+vis = dict(var_name='VIS 0.6µm', unit='[1]',
kind='MSG_4_SEVIRI_BDRF', sat_channel=1,
n_obs=961, obs_locations='square_array_evenly_on_grid',
# n_obs=1, obs_locations=[(44.141, -0.99)],
error_generate=0.03, error_assimilate=0.03,
- cov_loc_radius_km=20)
+ loc_horiz_km=20)
-wv62 = dict(plotname='Brightness temperature WV 6.2µm', plotunits='[K]',
+wv62 = dict(var_name='Brightness temperature WV 6.2µm', unit='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=5,
n_obs=961, obs_locations='square_array_evenly_on_grid',
error_generate=1., error_assimilate=2.,
- cov_loc_radius_km=20)
+ loc_horiz_km=20)
-wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
+wv73 = dict(var_name='Brightness temperature WV 7.3µm', unit='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=6,
n_obs=961, obs_locations='square_array_evenly_on_grid',
error_generate=1., error_assimilate=3.,
- cov_loc_radius_km=20)
+ loc_horiz_km=20)
-ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
+ir108 = dict(var_name='Brightness temperature IR 10.8µm', unit='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=9,
n_obs=1, obs_locations='square_array_evenly_on_grid',
error_generate=5., error_assimilate=10.,
- cov_loc_radius_km=32)
+ loc_horiz_km=32)
-radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
+radar = dict(var_name='Radar reflectivity', unit='[dBz]',
kind='RADAR_REFLECTIVITY',
n_obs=5776, obs_locations='square_array_evenly_on_grid',
error_generate=2.5, error_assimilate=2.5,
heights=range(2000, 14001, 2000),
- cov_loc_radius_km=1)
+ loc_horiz_km=20, loc_vert_km=2.5)
-t = dict(plotname='Temperature', plotunits='[K]',
+t = dict(var_name='Temperature', unit='[K]',
kind='RADIOSONDE_TEMPERATURE',
#n_obs=22500, obs_locations='square_array_evenly_on_grid',
n_obs=1, obs_locations=[(45., 0.)],
error_generate=0.2, error_assimilate=0.2,
heights=[1000,], #range(1000, 17001, 2000),
- cov_loc_radius_km=50)
+ loc_horiz_km=50, loc_vert_km=2.5)
-q = dict(plotname='Specific humidity', plotunits='[kg/kg]',
+q = dict(var_name='Specific humidity', unit='[kg/kg]',
kind='RADIOSONDE_SPECIFIC_HUMIDITY', n_obs=1,
error_generate=0., error_assimilate=5*1e-5,
heights=[1000], #range(1000, 17001, 2000),
- cov_loc_radius_km=0.1)
+ loc_horiz_km=0.1, loc_vert_km=2.5)
-t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
+t2m = dict(var_name='SYNOP Temperature', unit='[K]',
kind='SYNOP_TEMPERATURE', n_obs=1,
error_generate=0.1, error_assimilate=0.1,
- cov_loc_radius_km=40)
+ loc_horiz_km=40, loc_vert_km=2.5)
-psfc = dict(plotname='SYNOP Pressure', plotunits='[Pa]',
+psfc = dict(var_name='SYNOP Pressure', unit='[Pa]',
kind='SYNOP_SURFACE_PRESSURE', n_obs=1,
error_generate=50., error_assimilate=100.,
- cov_loc_radius_km=32)
+ loc_horiz_km=32, loc_vert_km=5)
-exp.observations = [t2m]
+exp.observations = [vis]
exp.update_vars = ['U', 'V', 'W', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC']
-#exp.update_vars = ['U', 'V', 'W', 'T', 'PH', 'MU', 'QVAPOR', 'PSFC']
+
diff --git a/cycled_exp.py b/cycled_exp.py
index 243ff26..8e28725 100755
--- a/cycled_exp.py
+++ b/cycled_exp.py
@@ -1,22 +1,20 @@
#!/usr/bin/python3
import os, sys, shutil, glob, warnings
import datetime as dt
-
-from dartwrf.utils import script_to_str
from dartwrf.workflows import WorkFlows
if __name__ == "__main__":
"""
Run a cycled OSSE with WRF and DART.
"""
- w = WorkFlows(exp_config='cfg.py', server_config='jet.py')
+ w = WorkFlows(exp_config='exp_example.py', server_config='jet.py')
timedelta_integrate = dt.timedelta(minutes=15)
timedelta_btw_assim = dt.timedelta(minutes=15)
id = None
- if True: # warm bubble
+ if False: # warm bubble
prior_path_exp = '/jetfs/home/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_noDA'
init_time = dt.datetime(2008, 7, 30, 12)
@@ -28,13 +26,13 @@ if __name__ == "__main__":
# id = w.run_ideal(depends_on=id)
# id = w.wrfinput_insert_wbubble(depends_on=id)
- if False: # random
+ if True: # random
prior_path_exp = '/jetfs/home/lkugler/data/sim_archive/exp_v1.19_P2_noDA'
init_time = dt.datetime(2008, 7, 30, 12)
time = dt.datetime(2008, 7, 30, 13)
- last_assim_time = dt.datetime(2008, 7, 30, 14)
- forecast_until = dt.datetime(2008, 7, 30, 18)
+ last_assim_time = dt.datetime(2008, 7, 30, 13)
+ forecast_until = dt.datetime(2008, 7, 30, 13, 15)
w.prepare_WRFrundir(init_time)
# id = w.run_ideal(depends_on=id)
diff --git a/dartwrf/assim_synth_obs.py b/dartwrf/assim_synth_obs.py
index 7fcc052..ddc39b3 100755
--- a/dartwrf/assim_synth_obs.py
+++ b/dartwrf/assim_synth_obs.py
@@ -8,75 +8,13 @@ from dartwrf.obs import error_models as err
import dartwrf.create_obsseq as osq
from dartwrf import wrfout_add_geo
from dartwrf import obsseq
+from dartwrf import dart_nml
from config.cfg import exp
from config.cluster import cluster
-earth_radius_km = 6370
wrfout_format = 'wrfout_d01_%Y-%m-%d_%H:%M:%S'
-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: # if not compute posterior
- template = cluster.scriptsdir + "/../templates/input.eval.nml"
- else:
- template = cluster.scriptsdir + "/../templates/input.nml"
- copy(template, cluster.dartrundir + "/input.nml")
-
- # impact_factor
- if exp.adjust_obs_impact:
- copy(cluster.obs_impact_filename, cluster.dartrundir + "/control_impact_runtime.txt")
-
- # The keys in `options` are placeholders in input.nml which will be replaced.
- # How? This is defined here
- options = {
- "<adjust_obs_impact>": '.true.' if exp.adjust_obs_impact else '.false.',
- "<filter_kind>": str(int(exp.filter_kind)),
- "<sampling_error_correction>": '.true.' if exp.sec else '.false.',
- "<prior_inflation>": str(exp.prior_inflation),
- "<post_inflation>": str(exp.post_inflation),
- "<n_ens>": str(int(exp.n_ens)),
- "<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),
- }
-
- # fail if horiz_dist_only == false but observations contain a satellite channel
- if exp.cov_loc_vert_km_horiz_km != False:
- for obscfg in exp.observations:
- if hasattr(obscfg, "sat_channel"):
- raise ValueError("Selected vertical localization, but observations contain satellite obs -> Not possible.")
-
- # Note: only one value of vertical localization possible
- try:
- 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)
- options["<horiz_dist_only>"] = ".false."
- except Exception as e:
- warnings.warn(str(e)+' - not using vertical localization.')
- options["<horiz_dist_only>"] = ".true."
- options["<vert_norm_hgt>"] = "99999.0" # dummy value
-
- for key, value in options.items():
- sed_inplace(cluster.dartrundir + "/input.nml", key, value)
-
- # input.nml for RTTOV
- rttov_nml = cluster.scriptsdir + "/../templates/obs_def_rttov.VIS.nml"
- append_file(cluster.dartrundir + "/input.nml", rttov_nml)
-
-
def link_nature_to_dart_truth(time):
"""Set a symlink from the WRFout file to be used as nature to the run_DART folder
@@ -363,7 +301,7 @@ def evaluate(assim_time,
if not os.path.isfile(cluster.dartrundir+'/obs_seq.out'):
raise RuntimeError(cluster.dartrundir+'/obs_seq.out does not exist')
- set_DART_nml(just_prior_values=True)
+ dart_nml.write_namelist(just_prior_values=True)
filter(nproc=6)
# archiving
@@ -531,7 +469,7 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
# remove any existing observation files
os.system("rm -f input.nml obs_seq.in obs_seq.out obs_seq.out-orig obs_seq.final")
- set_DART_nml()
+ dart_nml.write_namelist()
print("prepare nature")
prepare_nature_dart(time) # link WRF files to DART directory
@@ -556,7 +494,7 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
prepare_inflation_2(time, prior_init_time)
print(" 3) run filter ")
- set_DART_nml()
+ dart_nml.write_namelist()
filter(nproc=nproc)
archive_filteroutput(time)
diff --git a/dartwrf/create_obsseq.py b/dartwrf/create_obsseq.py
index c0a8a13..12aeab6 100755
--- a/dartwrf/create_obsseq.py
+++ b/dartwrf/create_obsseq.py
@@ -5,10 +5,8 @@ according to which observations are generated and subsequently assimilated.
import os, sys, warnings
import numpy as np
import datetime as dt
-import csv
from pysolar.solar import get_altitude, get_azimuth
-from config.cfg import exp
from config.cluster import cluster
from dartwrf.obs import calculate_obs_locations as col
from dartwrf import utils
diff --git a/dartwrf/dart_nml.py b/dartwrf/dart_nml.py
new file mode 100644
index 0000000..2271506
--- /dev/null
+++ b/dartwrf/dart_nml.py
@@ -0,0 +1,216 @@
+from dartwrf.utils import append_file
+
+from config.cfg import exp
+from config.cluster import cluster
+
+earth_radius_km = 6370
+
+def read_namelist(filepath):
+ """Read the DART namelist file into a dictionary.
+
+ Args:
+ filepath (str): Path to namelist file
+
+ Returns:
+ dict: keys are namelist sections, values are dictionaries of namelist variables
+ """
+
+ d = dict()
+ # read file into a list of strings
+ with open(filepath, 'r') as f:
+ lines = f.readlines()
+
+ for line in lines:
+ # skip whitespace
+ line = line.strip()
+
+ # skip comments
+ if not line.startswith('#'):
+
+ # skip empty lines
+ if len(line) > 0:
+
+ # namelist section
+ if line.startswith('&'):
+ section = line.lower()
+ d[section] = dict()
+
+ # namelist variable
+ else:
+ # split line into variable name and value
+ var, val = line.split('=')
+ var = var.strip().lower()
+ val = val.strip()
+
+ # split value into list if possible
+ if ',' in val:
+ val = val.split(',')
+ val = [v.strip() for v in val]
+
+ # add variable to dictionary
+ d[section][var] = val
+ return d
+
+
+def write_namelist_from_dict(d, filepath):
+ """Write a DART namelist dictionary to a file.
+
+ Args:
+ d (dict): keys are namelist sections, values are dictionaries of namelist variables
+ filepath (str): Path to namelist file
+ """
+ with open(filepath, 'w') as f:
+ for section in d:
+ f.write(section+'\n')
+ for var in d[section]:
+ val = d[section][var]
+ if isinstance(val, list):
+ val = ', '.join(val)
+ f.write('\t '+var+' = '+str(val)+'\n')
+ f.write('\t /\n\n')
+
+
+def _get_list_of_localizations():
+ """Compile the list of localizations for the DART namelist variables
+
+ Vertical localization can be defined in section &location_nml of 'input.nml'
+ using following namelist variables:
+ special_vert_normalization_obs_types (list of str)
+ special_vert_normalization_pressures (list of float)
+ special_vert_normalization_heights (list of float)
+ special_vert_normalization_levels (list of float)
+ special_vert_normalization_scale_heights (list of float)
+
+ To use scale height normalization, set obsdict['loc_vert_scaleheight'] = 0.5
+ To use height normalization, set obsdict['loc_vert_km'] = 3.0
+
+ Args:
+ exp (Experiment): Experiment object
+
+ Returns:
+ l_obstypes (list of str): entries for `special_vert_normalization_obs_types`
+ l_loc_horiz_rad (list of str): entries for `special_localization_cutoffs`
+ l_loc_vert_km (list of str): entries for `special_vert_normalization_heights`
+ l_loc_vert_scaleheight (list of str): entries for `special_vert_normalization_scale_heights`
+ """
+ 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
+
+ l_obstypes = []
+ l_loc_horiz_rad = []
+ l_loc_vert_km = []
+ l_loc_vert_scaleheight = []
+
+ for obscfg in exp.observations:
+
+ l_obstypes.append(obscfg["kind"])
+ loc_horiz_km = obscfg["loc_horiz_km"]
+
+ # compute horizontal localization
+ loc_horiz_rad = str(to_radian_horizontal(loc_horiz_km))
+ l_loc_horiz_rad.append(loc_horiz_rad)
+
+ # compute vertical localization
+
+ # choose either localization by height or by scale height
+ if hasattr(obscfg, "loc_vert_km") and hasattr(obscfg, "loc_vert_scaleheight"):
+ raise ValueError("Observation config contains both loc_vert_km and loc_vert_scaleheight. Please choose one.")
+
+ elif hasattr(obscfg, "loc_vert_km"): # localization by height
+ loc_vert_km = str(obscfg["loc_vert_km"])
+
+ vert_norm_hgt = to_vertical_normalization(loc_vert_km, loc_horiz_km)
+ l_loc_vert_km.append(vert_norm_hgt)
+
+ elif hasattr(obscfg, "loc_vert_scaleheight"): # localization by scale height
+ loc_vert_scaleheight = str(obscfg["loc_vert_scaleheight"])
+
+ # no conversion necessary, take the values as defined in obscfg
+ l_loc_vert_scaleheight.append(loc_vert_scaleheight)
+
+
+ # fail when both localization by height and scale height are requested
+ if len(l_loc_vert_km) > 0 and len(l_loc_vert_scaleheight) > 0:
+ raise ValueError("List of observation configurations contain both height and scale-height localization. Please choose one.")
+
+ # set the other (unused) list to a dummy value
+ if len(l_loc_vert_km) > 0:
+ l_loc_vert_scaleheight = [-1,]
+ else:
+ l_loc_vert_km = [-1,]
+
+ return l_obstypes, l_loc_horiz_rad, l_loc_vert_km, l_loc_vert_scaleheight
+
+
+def write_namelist(just_prior_values=False):
+ """Set DART namelist variables in 'input.nml' file.
+
+ 1. Takes the default namelist is the one already defined in the DART source code
+ 2. Calculates localization parameters from the observation configurations
+ 3. Overwrites other parameters as defined in the experiment configuration
+ 4. Writes the namelist to the DART run directory
+
+ Args:
+ just_prior_values (bool, optional): If True, only compute prior values, not posterior. Defaults to False.
+
+ Raises:
+ ValueError: If both height and scale-height localization are requested
+
+ Returns:
+ None
+ """
+ list_obstypes, list_loc_horiz_rad, list_loc_vert_km, list_loc_vert_scaleheight = _get_list_of_localizations()
+
+ nml = read_namelist(cluster.dart_srcdir + "/input.nml")
+
+ # make sure that observations defined in `exp.observations` are assimilated
+ nml['&obs_kind_nml']['assimilate_these_obs_types'] = list_obstypes
+
+ # dont compute posterior, just evaluate prior
+ if just_prior_values:
+ nml['&filter_nml']['compute_posterior'] = '.false.'
+ nml['&filter_nml']['output_members'] = '.false.'
+ nml['&filter_nml']['output_mean'] = '.false.'
+ nml['&filter_nml']['output_sd'] = '.false.'
+ nml['&obs_kind_nml']['assimilate_these_obs_types'] = []
+ nml['&obs_kind_nml']['evaluate_these_obs_types'] = list_obstypes
+
+
+ # write localization variables
+ nml['&assim_tools_nml']['special_localization_obs_types'] = list_obstypes
+ nml['&assim_tools_nml']['special_localization_cutoffs'] = list_loc_horiz_rad
+
+ nml['&location_nml']['special_vert_normalization_obs_types'] = list_obstypes
+ nml['&location_nml']['special_vert_normalization_heights'] = list_loc_vert_km
+ nml['&location_nml']['special_vert_normalization_scale_heights'] = list_loc_vert_scaleheight
+
+
+ # overwrite namelist with DART-WRF/config/ configuration
+ for key, value in exp.dart_nml.items():
+
+ # if key is not in namelist, add it
+ if key not in nml:
+ nml[key] = {}
+
+ # overwrite entry in each dictionary
+ nml[key] = value
+
+
+ # final checks
+ # fail if horiz_dist_only == false but observations contain a satellite channel
+ if nml['&location_nml']['horiz_dist_only'] == '.false.':
+ for obscfg in exp.observations:
+ if hasattr(obscfg, "sat_channel"):
+ raise ValueError("Selected vertical localization, but observations contain satellite obs -> Not possible.")
+
+ # write to file
+ write_namelist_from_dict(nml, cluster.dartrundir + "/input.nml")
+
+ # append section for RTTOV
+ rttov_nml = cluster.scriptsdir + "/../templates/obs_def_rttov.VIS.nml"
+ append_file(cluster.dartrundir + "/input.nml", rttov_nml)
\ No newline at end of file
--
GitLab