consistent configurations

config/cfg.py → config/exp_example.py

 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']
+
 

cycled_exp.py

 #!/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)

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)
 

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

dartwrf/dart_nml.py

+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