synop obs, faster clouddep obserr, outlier thresh

config/cfg.py

@@ -9,7 +9,7 @@ class ExperimentConfiguration(object):
 
 
 exp = ExperimentConfiguration()
-exp.expname = "exp_v1.12_LMU_so_VIS2"
+exp.expname = "exp_v1.12_conv1"
 exp.model_dx = 2000
 exp.timestep = 10
 exp.n_ens = 40
@@ -17,17 +17,24 @@ exp.n_nodes = 10
 
 n_obs = 64  # radar: n_obs for each observation height level
 
-vis = dict(sat=True, channel=1, n_obs=n_obs, err_std=0.03,
+vis = dict(sat_channel=1, n_obs=n_obs, err_std=0.03,
            cov_loc_radius_km=10)
-wv = dict(sat=True, channel=6, n_obs=n_obs, err_std=5.,
+wv = dict(sat_channel=6, n_obs=n_obs, err_std=5.,
            cov_loc_radius_km=10)
-
-radar = dict(sat=False, kind='RADAR', n_obs=n_obs, err_std=5.,
-             cov_loc_radius_km=10, cov_loc_vert_km=1)
-psfc = dict(sat=False, kind='PSFC', n_obs=n_obs, err_std=50.,
+ir108 = dict(sat_channel=9, n_obs=n_obs, err_std=5.,
              cov_loc_radius_km=10)
 
-exp.observations = [vis, ]
+radar = dict(kind='RADAR', n_obs=n_obs, err_std=5.,
+             heights=np.arange(1000, 15001, 1000),
+             cov_loc_radius_km=10, cov_loc_vert_km=2)
+
+t2m = dict(kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=1.0, 
+           cov_loc_radius_km=32, cov_loc_vert_km=1)
+psfc = dict(kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
+           cov_loc_radius_km=32)
+           
+
+exp.observations = [t2m, psfc, ]
 
 # directory paths depend on the name of the experiment
 cluster.expname = exp.expname

scheduler.py

@@ -16,12 +16,18 @@ from scripts.utils import script_to_str, symlink
 # allow scripts to access the configuration
 symlink(cluster.scriptsdir+'/../config', cluster.scriptsdir+'/config')
 
+log_dir = cluster.archivedir()+'/logs/'
+slurm_scripts_dir = cluster.archivedir()+'/slurm-scripts/'
+print('logging to', log_dir)
+print('scripts, which are submitted to SLURM:', slurm_scripts_dir)
+
 def my_Slurm(*args, cfg_update=dict(), **kwargs):
-    """Shortcut to slurmpy's class; keep default kwargs
-    and only update according to `cfg_update`
+    """Shortcut to slurmpy's class; keep certain default kwargs
+    and only update some with kwarg `cfg_update`
     see https://github.com/brentp/slurmpy
     """
-    return Slurm(*args, slurm_kwargs=dict(cluster.slurm_cfg, **cfg_update), **kwargs)
+    return Slurm(*args, slurm_kwargs=dict(cluster.slurm_cfg, **cfg_update), 
+                 log_dir=log_dir, scripts_dir=slurm_scripts_dir, **kwargs)
 
 class Cmdline(object):
     def __init__(self, name, cfg_update):
@@ -31,18 +37,16 @@ class Cmdline(object):
         print('running', self.name, 'without SLURM')
         os.system(cmd)
 
-def clear_logs(backup_existing_to_archive=True):
-    dirs = ['/logs/', '/slurm-scripts/']
-    for d in dirs:
-        archdir = cluster.archivedir()+d
-        if backup_existing_to_archive:
-            os.makedirs(archdir, exist_ok=True)
-        dir = cluster.scriptsdir+'/../'+d
-        for f in os.listdir(dir):
-            if backup_existing_to_archive:
-                shutil.move(dir+f, archdir+f)
-            else:
-                os.remove(dir+f)
+def backup_scripts():
+    os.makedirs(cluster.archivedir(), exist_ok=True)
+    dir_backup_scripts = cluster.archivedir()+'/DART-WRF/'
+    os.makedirs(dir_backup_scripts, exist_ok=True)
+    os.makedirs(dir_backup_scripts+'/old/', exist_ok=True)
+    os.system('mv '+dir_backup_scripts+'/* '+dir_backup_scripts+'/old/')  # move old logs
+    shutil.copy(cluster.scriptsdir+'/../scheduler.py', dir_backup_scripts+'/scheduler.py')
+    shutil.copy(cluster.scriptsdir+'/../config/clusters.py', dir_backup_scripts+'/clusters.py')
+    shutil.copy(cluster.scriptsdir+'/../config/cfg.py', dir_backup_scripts+'/cfg.py')
+    os.system('cp -r '+cluster.scriptsdir+' '+dir_backup_scripts+'/')
 
 def prepare_wrfinput():
     """Create WRF/run directories and wrfinput files
@@ -71,7 +75,7 @@ done
     id = s.run(cmd, depends_on=[id])
     return id
 
-def update_wrfinput_from_archive(time, background_init_time, exppath, depends_on=None):
+def update_wrfinput_from_archive(valid_time, background_init_time, exppath, depends_on=None):
     """Given that directories with wrfinput files exist,
     update these wrfinput files according to wrfout files
     """
@@ -79,7 +83,7 @@ def update_wrfinput_from_archive(time, background_init_time, exppath, depends_on
 
     # path of initial conditions, <iens> is replaced by member index
     IC_path = exppath + background_init_time.strftime('/%Y-%m-%d_%H:%M/')  \
-              +'*iens*/'+time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
+              +'*iens*/'+valid_time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
     id = s.run(cluster.python+' '+cluster.scriptsdir+'/update_wrfinput_from_wrfout.py '
                 +IC_path, depends_on=[depends_on])
     return id
@@ -92,12 +96,14 @@ def run_ENS(begin, end, depends_on=None):
 
     # first minute forecast (needed for validating an assimilation)
     hist_interval = 1
+    radt = 1  # calc CFRAC also in first minute
     begin_plus1 = begin+dt.timedelta(minutes=1)
     s = my_Slurm("preWRF1", cfg_update=dict(time="2"))
-    id = s.run(cluster.python+' '+cluster.scriptsdir+'/prepare_namelist.py '
-            + begin.strftime('%Y-%m-%d_%H:%M')+' '
-            + begin_plus1.strftime('%Y-%m-%d_%H:%M')+' '
-            + str(hist_interval), 
+    id = s.run(' '.join([cluster.python,
+               cluster.scriptsdir+'/prepare_namelist.py',
+               begin.strftime('%Y-%m-%d_%H:%M'),
+               begin_plus1.strftime('%Y-%m-%d_%H:%M'),
+               str(hist_interval), str(radt),]) 
             depends_on=[id])
 
     s = my_Slurm("runWRF1", cfg_update={"nodes": "1", "array": "1-"+str(exp.n_nodes),
@@ -114,11 +120,13 @@ def run_ENS(begin, end, depends_on=None):
 
     # whole forecast timespan
     hist_interval = 5
+    radt = 5
     s = my_Slurm("preWRF2", cfg_update=dict(time="2"))
-    id = s.run(cluster.python+' '+cluster.scriptsdir+'/prepare_namelist.py '
-            + begin.strftime('%Y-%m-%d_%H:%M')+' '
-            + end.strftime('%Y-%m-%d_%H:%M')+' '
-            + str(hist_interval), 
+    id = s.run(' '.join([cluster.python,
+               cluster.scriptsdir+'/prepare_namelist.py',
+               begin.strftime('%Y-%m-%d_%H:%M'),
+               begin_plus1.strftime('%Y-%m-%d_%H:%M'),
+               str(hist_interval), str(radt),])
             depends_on=[id])
 
     time_in_simulation_hours = (end-begin).total_seconds()/3600
@@ -134,43 +142,40 @@ def run_ENS(begin, end, depends_on=None):
     #           + begin.strftime('%Y-%m-%d_%H:%M'), depends_on=[id2])
     return id
 
-def assimilate(assim_time, background_init_time,
-               prior_from_different_exp=False, depends_on=None):
+def assimilate(assim_time, prior_init_time,
+               prior_path_exp=False, depends_on=None):
     """Creates observations from a nature run and assimilates them.
 
     Args:
         assim_time (dt.datetime): timestamp of prior wrfout files
-        background_init_time (dt.datetime): 
+        prior_init_time (dt.datetime): 
             timestamp to find the directory where the prior wrfout files are
-        prior_from_different_exp (bool or str):
+        prior_path_exp (bool or str):
             put a `str` to take the prior from a different experiment
             if False: use `archivedir` (defined in config) to get prior state
             if str: use this directory to get prior state
     """
-    id = depends_on
-
-    if prior_from_different_exp:
-        prior_expdir = prior_from_different_exp
-    else:
+    if prior_path_exp == False:
         prior_expdir = cluster.archivedir()
+    elif not isinstance(prior_path_exp, str):
+        raise TypeError('prior_path_exp either str or False, is '+str(type(prior_path_exp)))
 
     # prepare state of nature run, from which observation is sampled
-    s = my_Slurm("prepNature", cfg_update=dict(time="2"))
-    id = s.run(cluster.python+' '+cluster.scriptsdir+'/prepare_nature.py '
-               +time.strftime('%Y-%m-%d_%H:%M'), depends_on=[depends_on])
+    #s = my_Slurm("prepNature", cfg_update=dict(time="2"))
+    #id = s.run(cluster.python+' '+cluster.scriptsdir+'/prepare_nature.py '
+    #           +time.strftime('%Y-%m-%d_%H:%M'), depends_on=[depends_on])
 
     # prepare prior model state
     s = my_Slurm("preAssim", cfg_update=dict(time="2"))
     id = s.run(cluster.python+' '+cluster.scriptsdir+'/pre_assim.py '
                +assim_time.strftime('%Y-%m-%d_%H:%M ')
-               +background_init_time.strftime('%Y-%m-%d_%H:%M ')
-               +prior_expdir, depends_on=[id])
+               +prior_init_time.strftime('%Y-%m-%d_%H:%M ')
+               +prior_path_exp, depends_on=[depends_on])
 
-    # generate observations
-    s = my_Slurm("genSynthObs", cfg_update=dict(ntasks="48", time="10"))
+    # prepare nature run, generate observations
+    s = my_Slurm("genSynthObs", cfg_update=dict(ntasks="48", time="12"))
     id = s.run(cluster.python+' '+cluster.scriptsdir+'/gen_synth_obs.py '
-               +time.strftime('%Y-%m-%d_%H:%M'),
-               depends_on=[id])
+               +time.strftime('%Y-%m-%d_%H:%M'), depends_on=[id])
  
     # actuall assimilation step
     s = my_Slurm("Assim", cfg_update=dict(ntasks="48", time="50", mem="200G"))
@@ -184,8 +189,8 @@ def assimilate(assim_time, background_init_time,
     s = my_Slurm("updateIC", cfg_update=dict(time="8"))
     id = s.run(cluster.python+' '+cluster.scriptsdir+'/update_wrfinput_from_filteroutput.py '
                 +assim_time.strftime('%Y-%m-%d_%H:%M ')
-                +background_init_time.strftime('%Y-%m-%d_%H:%M ')
-                +prior_expdir, depends_on=[id])
+                +prior_init_time.strftime('%Y-%m-%d_%H:%M ')
+                +prior_path_exp, depends_on=[id])
     return id
 
 
@@ -206,7 +211,7 @@ print('starting osse')
 timedelta_integrate = dt.timedelta(minutes=45)
 timedelta_btw_assim = dt.timedelta(minutes=30)
 
-clear_logs(backup_existing_to_archive=False)
+backup_scripts()
 id = None
 
 start_from_existing_state = True
@@ -225,32 +230,33 @@ if is_new_run:
     first_guess = False
     
 elif start_from_existing_state:
-    id = prepare_wrfinput()  # create initial conditions
+    #id = prepare_wrfinput()  # create initial conditions
     
     # get initial conditions from archive
     background_init_time = dt.datetime(2008, 7, 30, 11, 30)
     time = dt.datetime(2008, 7, 30, 12)
-    exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.12_LMU_so_VIS2'
+    exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.12_LMU_so_VIS'
     first_guess = exppath_arch
-    id = update_wrfinput_from_archive(time, background_init_time, exppath_arch, depends_on=id)
+    #id = update_wrfinput_from_archive(time, background_init_time, exppath_arch, depends_on=id)
 
-while time <= dt.datetime(2008, 7, 30, 15):
+while time <= dt.datetime(2008, 7, 30, 14):
     assim_time = time
+    prev_forecast_init = assim_time - timedelta_btw_assim
+    this_forecast_init = assim_time  # start integration from here
+    this_forecast_end = assim_time + timedelta_integrate
+
     id = assimilate(assim_time,
-                     background_init_time,
-                     prior_from_different_exp=first_guess,
+                    prev_forecast_init,
+                    prior_path_exp=first_guess,
                     depends_on=id)
-    first_guess = False
+    #first_guess = False
 
-    prev_forecast_init = assim_time - timedelta_btw_assim
-    this_forecast_init = assim_time  # start integration from here
-    this_forecast_end = assim_time + timedelta_btw_assim
+    #id = update_wrfinput_from_archive(this_forecast_init, this_forecast_init, cluster.archivedir(), depends_on=id)
 
     id = run_ENS(begin=this_forecast_init,
-                end=this_forecast_end,
+                 end=this_forecast_init+dt.timedelta(minutes=10),
                  depends_on=id)
     time += timedelta_btw_assim
-
-    create_satimages(depends_on=id)
+    #create_satimages(depends_on=id)
 
 mailme(id)

scripts/create_obsseq.py

@@ -297,23 +297,32 @@ def calc_obs_locations_3d(coords, heights):
     return coords2
 
 
-def generic_obs(obs_type, time_dt, coords, error_var, output_path='./'):
-
-    obs_codes = {'RASO_T': {'name': 'RADIOSONDE_TEMPERATURE', 'nr': '5'},
-                 'RADAR': {'name': 'RADAR_REFLECTIVITY', 'nr': '37'},
-                 'PSFC': {'name': 'LAND_SFC_PRESSURE', 'nr': '29'},}
+def generic_obs(obs_kind, time_dt, coords, error_var, heights=False, output_path='./'):
+
+    obs_kind_nrs = {'RADIOSONDE_TEMPERATURE': '5',
+                    'RADAR_REFLECTIVITY': '37',
+                    'SYNOP_SURFACE_PRESSURE': '94',
+                    'SYNOP_SPECIFIC_HUMIDITY': '95',
+                    'SYNOP_TEMPERATURE': '96',
+                    }
+
+    if 'SYNOP' in obs_kind:
+        is_sfc_obs = True
+        heights = [2,]
+    else:
+        is_sfc_obs = False
 
-    heights = [5000., ] #np.arange(5000, 15001, 1000)
+    if not heights:
+        heights = [5000., ]
     coords = calc_obs_locations_3d(coords, heights)
 
     dart_date_day, secs_thatday = get_dart_date(add_timezone_UTC(time_dt))
     print('secs, days:', secs_thatday, dart_date_day)
 
-    obs_name = obs_codes[obs_type]['name']
-    obs_kind_nr = obs_codes[obs_type]['nr']
-
-    write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords,
-                         dart_date_day, secs_thatday, output_path)
+    obs_kind_nr = obs_kind_nrs[obs_kind]
+    write_generic_obsseq(obs_kind, obs_kind_nr, error_var, coords,
+                         dart_date_day, secs_thatday, output_path,
+                         vert_coord_sfc=is_sfc_obs)
 
 
 if __name__ == '__main__':

scripts/gen_synth_obs.py

@@ -2,9 +2,11 @@ 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 as osq
+import wrfout_add_geo
 
 earth_radius_km = 6370
 
@@ -51,6 +53,19 @@ def read_prior_obs(f_obs_prior):
             OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens)))
     return OBSs
 
+def read_obsseqout(f):
+    obsseq = open(f, 'r').readlines()
+    true = []
+    obs = []
+    # read observations from obs_seq.out
+    for i, line in enumerate(obsseq):
+        if ' OBS ' in line:
+            observed = float(obsseq[i+1].strip())
+            truth = float(obsseq[i+2].strip())
+            true.append(truth)
+            obs.append(observed)
+    return true, obs
+
 def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
     """
     overwrite observation errors in a obs_seq.out file
@@ -107,6 +122,41 @@ def set_input_nml(sat_channel=False, just_prior_values=False,
         append_file(cluster.dartrundir+'/input.nml', rttov_nml)
 
 
+def obs_operator_ensemble():
+    os.chdir(cluster.dartrundir)
+
+    if sat_channel:
+        list_ensemble_truths = []
+
+        for iens in range(1, exp.n_ens+1):
+            print('observation operator for ens #'+str(iens))
+            # ens members are already linked to advance_temp<i>/wrfout_d01
+            copy(cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01',
+                 cluster.dartrundir+'/wrfout_d01')
+            
+            t = dt.datetime.now()
+            wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
+            print((dt.datetime.now()-t).total_seconds(), 'secs for adding geodata')
+
+            # DART may need a wrfinput file as well, which serves as a template for dimension sizes
+            symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
+
+            # run perfect_model obs (forward operator)
+            os.system('mpirun -np 12 ./perfect_model_obs > /dev/null')
+
+            # truth values in obs_seq.out are H(x) values
+            vals, _ = read_obsseqout(cluster.dartrundir+'/obs_seq.out')
+            list_ensemble_truths.append(vals)
+        
+        n_obs = len(list_ensemble_truths[0])
+        np_array = np.full((exp.n_ens, n_obs), np.nan)
+        for i in range(exp.n_ens):
+            np_array[i, :] = list_ensemble_truths[i]
+        return np_array
+    else:
+        raise NotImplementedError()
+
+
 if __name__ == "__main__":
 
     time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
@@ -115,80 +165,101 @@ if __name__ == "__main__":
     # remove any existing observation files
     os.chdir(cluster.dartrundir); os.system('rm -f obs_seq_*.out obs_seq.in obs_seq.final')
 
+    def prepare_nature_dart():
+        # get wrfout_d01 from nature run
+        shutil.copy(time.strftime(cluster.nature_wrfout),
+                    cluster.dartrundir+'/wrfout_d01')
+
+        wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
+
+        # DART may need a wrfinput file as well, which serves as a template for dimension sizes
+        symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
+
+    prepare_nature_dart()
+
     # loop over observation types
     for i_obs, obscfg in enumerate(exp.observations):
 
         n_obs = obscfg['n_obs']
         error_var = (obscfg['err_std'])**2
-        sat_channel = obscfg.get('channel', False)
+        sat_channel = obscfg.get('sat_channel', False)
         cov_loc = obscfg['cov_loc_radius_km']
         dist_obs = obscfg.get('distance_between_obs_km', False)
         cov_loc_vert_km = obscfg.get('cov_loc_vert_km', False)
+        heights = obscfg.get('heights', False)
 
         # generate obs_seq.in
         obs_coords = osq.calc_obs_locations(n_obs, coords_from_domaincenter=False, 
                                             distance_between_obs_km=dist_obs, 
                                             fpath_obs_locations=fpath_obs_coords)
 
-        if obscfg['sat']:
+        if sat_channel:
             osq.sat(time, sat_channel, obs_coords, error_var,
                     output_path=cluster.dartrundir)
         else:
             osq.generic_obs(obscfg['kind'], time, obs_coords, error_var,
+                            heights=heights, 
                             output_path=cluster.dartrundir)
 
         if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
             raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
 
-        # generate observations (obs_seq.out)
+        os.chdir(cluster.dartrundir)
+        if sat_channel == 6:
+            """cloud dependent observation error
+
+            # methodologically:
+            1) gen H(x_nature)
+            2) gen H(x_prior)
+            3) find the observation error for a pair of (H(x_nature), H(x_background))
+            4) generate actual observation 
+                with the observation error as function of H(x_nature) and H(x_background)
+
+            # technically:
+            4) the file 'obs_seq.in' needs to be edited to show corrected observation errors
+            """
+            # 1) gen H(x_nature)
             set_input_nml(sat_channel=sat_channel,
                             cov_loc_radius_km=cov_loc,
                             cov_loc_vert_km=cov_loc_vert_km)
-        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())
-
-        # cloud dependent observation error
-        if obscfg['sat']:
-            if sat_channel == 6:
-                # run ./filter to have prior observation estimates from model state
-                set_input_nml(sat_channel=sat_channel,
-                              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())
+            Hx_nature, _ = read_obsseqout(cluster.dartrundir+'/obs_seq.out')
 
-                # 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)
+            # 2) gen H(x_prior) for the whole ensemble 
+            Hx_prior = obs_operator_ensemble()  # files are already linked to DART directory
 
-                # compute the observation error necessary
-                # to achieve a certain operational FGD distribution
+            # 3) find the observation error for a pair of (H(x_nature), H(x_background))
+            # necessary to achieve a certain FGD distribution which is near to operational
+            n_obs = len(Hx_nature)
             OEs = []
-                for obs in OBSs:
-                    bt_y = obs['truth']
-                    bt_x_ens = obs['prior_ens']
+            for iobs in range(n_obs):
+
+                bt_y = Hx_nature[iobs]
+                bt_x_ens = Hx_prior[:,iobs]
                 CIs = [cloudimpact_73(bt_x, bt_y) for bt_x in bt_x_ens]
+                mean_CI = np.mean(CIs)
 
-                    oe_nature = oe_73(np.mean(CIs))
+                oe_nature = oe_73(mean_CI)
+                print('oe_nature=', oe_nature, ' K')
                 OEs.append(oe_nature)
 
-                # write obs_seq.out
+            # correct obs_err in obs_seq.in (to produce actual observations later on)
             fpath_obsseqout = cluster.dartrundir+'/obs_seq.in'
             edit_obserr_in_obsseq(fpath_obsseqout, OEs)
 
-                # 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())
+            # ensure correct nature file linked 
+            # nature should be the real nature again (was changed in the meantime)
+            prepare_nature_dart()   
                 
         # correct input.nml for actual assimilation later on
         set_input_nml(sat_channel=sat_channel,
                         cov_loc_radius_km=cov_loc,
                         cov_loc_vert_km=cov_loc_vert_km)
 
+        # 4) generate actual observations (with correct error)
+        os.chdir(cluster.dartrundir)
+        os.system('mpirun -np 12 ./perfect_model_obs')
+
         # rename according to i_obs
         os.rename(cluster.dartrundir+'/obs_seq.out', 
                   cluster.dartrundir+'/obs_seq_'+str(i_obs)+'.out')

scripts/prepare_namelist.py

@@ -5,7 +5,7 @@ sys.path.append(os.getcwd())
 from config.cfg import exp, cluster
 from utils import sed_inplace, copy, symlink, mkdir
 
-def run(cluster, iens, begin, end, hist_interval=5):
+def run(cluster, iens, begin, end, hist_interval=5, radt=5):
     rundir = cluster.wrf_rundir(iens)
     print(rundir)
     copy(cluster.namelist, rundir+'/namelist.input')
@@ -13,6 +13,7 @@ def run(cluster, iens, begin, end, hist_interval=5):
     sed_inplace(rundir+'/namelist.input', '<dx>', str(int(exp.model_dx)))
     sed_inplace(rundir+'/namelist.input', '<timestep>', str(int(exp.timestep)))
     sed_inplace(rundir+'/namelist.input', '<hist_interval>', str(int(hist_interval)))
+    sed_inplace(rundir+'/namelist.input', '<radt>', str(int(radt)))
 
     archdir = cluster.archivedir()+begin.strftime('/%Y-%m-%d_%H:%M/'+str(iens)+'/')
     print('namelist for run from', begin, end, 'output to', archdir)
@@ -42,8 +43,9 @@ def run(cluster, iens, begin, end, hist_interval=5):
 if __name__ == '__main__':
     begin = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
     end = dt.datetime.strptime(sys.argv[2], '%Y-%m-%d_%H:%M')
-    hist_interval = int(sys.argv[3])
+    intv = int(sys.argv[3])
+    radt = int(sys.argv[4])
 
     print('prepare namelists for all ens members')
     for iens in range(1, exp.n_ens+1):
-        run(cluster, iens, begin, end, hist_interval)
+        run(cluster, iens, begin, end, hist_interval=intv, radt=radt)

scripts/prepare_wrfinput.py

@@ -5,11 +5,6 @@ from utils import symlink, copy, link_contents
 
 import prepare_namelist
 
-# archive configuration
-os.makedirs(cluster.archivedir(), exist_ok=True)
-shutil.copy(cluster.scriptsdir+'/../config/clusters.py', cluster.archivedir()+'/clusters.py')
-shutil.copy(cluster.scriptsdir+'/../config/cfg.py', cluster.archivedir()+'/cfg.py')
-
 for iens in range(1, exp.n_ens+1):
     print('preparing ens', iens)
     input_prof = (cluster.input_profile).replace('<iens>', str(iens).zfill(3))

scripts/update_wrfinput_from_wrfout.py

@@ -18,7 +18,7 @@ for iens in range(1, exp.n_ens+1):
     wrfin = cluster.wrf_rundir(iens)+'/wrfinput_d01'
 
     print('updating', wrfin, 'to state in', wrfout)
-    assert os.path.isfile(wrfout)
+    assert os.path.isfile(wrfout), wrfout
 
     # overwrite variables in wrfinput file
     os.system(cluster.ncks+' -A -v '+vars+' '+wrfout+' '+wrfin)

templates/input.nml

@@ -71,7 +71,7 @@
 
 &quality_control_nml
    input_qc_threshold          = 3.0,
-   outlier_threshold           = 3.0,
+   outlier_threshold           = -1.0,
    enable_special_outlier_code = .false.
    /
 
@@ -281,7 +281,7 @@
    layout = 1,
    tasks_per_node = 48
    communication_configuration = 1
-   debug                       = .true.
+   debug                       = .false.
    /
 
 &obs_def_gps_nml

templates/namelist.input

@@ -52,7 +52,7 @@
  mp_physics                          = 8,     1,     1,
  ra_lw_physics                       = 4,     0,     0,
  ra_sw_physics                       = 4,     0,     0,
- radt                                = 5,    30,    30,
+ radt                                = <radt>,    30,    30,
  sf_sfclay_physics                   = 5,     0,     0,
  sf_surface_physics                  = 2,     0,     0,
  bl_pbl_physics                      = 5,     0,     0,