unfinished

config/jet.py

@@ -2,7 +2,7 @@
 
 cluster_defaults = dict(
     max_nproc = 20,
-    max_nproc_for_each_ensemble_member = 16,
+    max_nproc_for_each_ensemble_member = 9,
     use_slurm = True,
 
     # binaries
@@ -20,7 +20,7 @@ cluster_defaults = dict(
 
     # paths used as input
     dir_wrf_src = '/jetfs/home/lkugler/data/compile/WRF-4.3/run',
-    dir_dart_src = '/jetfs/home/lkugler/data/compile/DART/DART-10.8.3_10pct/models/wrf/work/',
+    dir_dart_src = '/jetfs/home/lkugler/data/compile/DART/DART-10.8.3/models/wrf/work/',
     dir_rttov_src = '/jetfs/home/lkugler/data/compile/RTTOV13/rtcoef_rttov13/',
     dir_dartwrf_dev = '/jetfs/home/lkugler/DART-WRF/',
     

dartwrf/evaluate_obs_space.py

@@ -21,9 +21,9 @@ def evaluate_one_time(assim_time, valid_time, use_other_obsseq=False):
      aso.prepare_prior_ensemble(valid_time, prior_init_time=assim_time, prior_valid_time=valid_time, prior_path_exp=cluster.archivedir)
      dart_nml.write_namelist()
 
-     if os.path.isfile(exp.use_existing_obsseq):
+     if os.path.isfile(exp.assimilate_existing_obsseq):
           # use the existing obs_seq.out file
-          shutil.copy(exp.use_existing_obsseq, cluster.dart_rundir+'/obs_seq.out')
+          shutil.copy(exp.assimilate_existing_obsseq, cluster.dart_rundir+'/obs_seq.out')
      else:
           # is there a pre-existing obs file from assimilation before?
           f_oso = valid_time.strftime(cluster.pattern_obs_seq_out)

dartwrf/evaluate_prior.py

@@ -21,8 +21,8 @@ if __name__ == "__main__":
      aso.prepare_prior_ensemble(assim_time, prior_init, prior_valid_time, cluster.archive_base+'/'+prior_exp)
 
     # prepare an obsseq without rejected observations
-     if exp.use_existing_obsseq:  # from another exp
-          oso_input = assim_time.strftime(exp.use_existing_obsseq)
+     if exp.assimilate_existing_obsseq:  # from another exp
+          oso_input = assim_time.strftime(exp.assimilate_existing_obsseq)
      else:  # from same exp
           oso_input = cluster.archivedir+'/obs_seq_out' + assim_time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out-beforeQC")
           if not os.path.isfile(oso_input):

docs/source/tutorial1.rst

@@ -61,7 +61,7 @@ In case you want to generate new observations, like for an observing system simu
 
 .. code-block:: python
 
-    exp.use_existing_obsseq = False
+    exp.assimilate_existing_obsseq = False
     
 in this case, you need to set the path to WRF nature run files from where DART can generate observations:
 
@@ -77,7 +77,7 @@ You can use pre-existing observation files with
 
 .. code-block:: python
 
-    exp.use_existing_obsseq = '/usr/data/sim_archive/exp_ABC/obs_seq_out/%Y-%m-%d_%H:%M_obs_seq.out'
+    exp.assimilate_existing_obsseq = '/usr/data/sim_archive/exp_ABC/obs_seq_out/%Y-%m-%d_%H:%M_obs_seq.out'
     
 where time-placeholders (``%Y-%m-%d_%H:%M``) are filled in later, depending on the assimilation time.
 

free_forecast2.py

+import sys
+import datetime as dt
+import pandas as pd
+from dartwrf.workflows import WorkFlows
+from dartwrf.utils import Config
+
+
+# import default config for jet
+from config.jet import cluster_defaults
+from config.defaults import dart_nml, CF_config
+
+
+ensemble_size = 1
+t0 = dt.datetime(2008, 7, 30, 8)
+
+for c_s in [0.15, 0.2]:
+
+    id = None
+    cfg = Config(name='exp_nat250m_Cs'+str(c_s), 
+        model_dx = 2000,
+        ensemble_size = ensemble_size,
+        geo_em_forecast = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
+        time = t0,
+        c_s = str(c_s),
+
+        input_profile = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/2008-07-30_08:00/1/input_sounding',
+        verify_against = 'nat_250m_blockavg2km',
+        **cluster_defaults, # type: ignore
+    )
+
+
+    w = WorkFlows(cfg)
+    w.prepare_WRFrundir(cfg)
+    id = w.run_ideal(cfg, depends_on=id)
+
+    timedelta_integrate = dt.timedelta(hours=10)
+    restart_interval = 9999
+    
+    cfg.update( WRF_start=t0, 
+                WRF_end=t0+timedelta_integrate, 
+                restart=False, 
+                restart_interval=restart_interval,
+                hist_interval_s=300,
+    )
+    # only verification depends on the output
+    id_long = w.run_WRF(cfg, depends_on=id)
+    id_long = w.run_RTTOV(cfg, depends_on=id_long)
+        
+    w.verify(cfg, depends_on=id_long)
\ No newline at end of file

multiple_exps.py

+import sys
 import datetime as dt
 import pandas as pd
 from dartwrf.workflows import WorkFlows
@@ -6,89 +7,103 @@ from dartwrf.utils import Config
 
 # import default config for jet
 from config.jet import cluster_defaults
-from config.defaults import dart_nml
+from config.defaults import dart_nml, CF_config, vis
 
 # test multiple assimilation windows (11-12, 12-13, 13-14, )
 timedelta_btw_assim = dt.timedelta(minutes=15)
-assim_times_start = pd.date_range('2008-07-30 11:00', '2008-07-30 13:00', freq='1h')
 ensemble_size = 40
 
 dart_nml['&filter_nml'].update(num_output_state_members=ensemble_size,
                                 ens_size=ensemble_size)
 
-vis = dict(
-        kind='MSG_4_SEVIRI_BDRF', sat_channel=1,
-        km_between_obs=12, skip_border_km=8.0,
-        error_generate=0.03, error_assimilate=0.06,
-        loc_horiz_km=20,
-        height=6000, loc_vert_km=6,
-        )
+# which DART version to use?
+assimilate_cloudfractions = True
+cf1 = dict(kind='CF192km', loc_horiz_km=9999)
+cf2 = dict(kind='CF96km', loc_horiz_km=96)
+cf3 = dict(kind='CF48km', loc_horiz_km=48)
+cf4 = dict(kind='CF24km', loc_horiz_km=24)
+cf5 = dict(kind='CF12km', loc_horiz_km=12)
 
-for t0 in assim_times_start:
+assimilate_these_observations = [cf5]
 
-    id = None
-    cfg = Config(name=t0.strftime('exp_nat250_VIS_obs12_loc20_oe2_inf4-0.5_%H%M'), 
-        model_dx = 2000,
-        ensemble_size = ensemble_size,
-        dart_nml = dart_nml,
-        geo_em_forecast = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
-        time = t0,
-        
-        use_existing_obsseq = False,
-        nature_wrfout_pattern = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/*/1/wrfout_d01_%Y-%m-%d_%H_%M_%S',
-        #geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
-        geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.250m_1600x1600',
-        
-        update_vars = ['U', 'V', 'W', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'QSNOW', 'PSFC'],
-        input_profile = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/2008-07-30_08:00/1/input_sounding',
-        nature_exp_verif = 'nat_250m_blockavg2km',
-        assimilate_these_observations = [vis,],
-        **cluster_defaults, # type: ignore
+if assimilate_cloudfractions:
+    cluster_defaults.update(
+        dir_dart_src = '/jetfs/home/lkugler/data/compile/DART/DART-v10.9.2-andrea/models/wrf/work/',
     )
+else:
+    for obs in assimilate_these_observations:
+        if obs['kind'] == 'MSG_4_SEVIRI_BDRF':
+            cluster_defaults.update(
+                dir_dart_src = '/jetfs/home/lkugler/data/compile/DART/DART-10.8.3_10pct/models/wrf/work/',
+                )
+
 
-    w = WorkFlows(cfg)
-    w.prepare_WRFrundir(cfg)
-    #id = w.run_ideal(cfg, depends_on=id)
+
+time0 = dt.datetime(2008, 7, 30, 11)
+
+id = None
+cfg = Config(name='exp_nat250_VIS_SO12_loc12_oe2_inf4-0.5', 
+    model_dx = 2000,
+    ensemble_size = ensemble_size,
+    dart_nml = dart_nml,
+    geo_em_forecast = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
+    time = time0,
+    
+    assimilate_these_observations = assimilate_these_observations,
+    assimilate_cloudfractions = assimilate_cloudfractions,
+    CF_config=CF_config,
+    
+    assimilate_existing_obsseq = False,
+    #nature_wrfout_pattern = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/*/1/wrfout_d01_%Y-%m-%d_%H_%M_%S',
+    #geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
+    # geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.250m_1600x1600',
     
-    # assimilate at these times
-    time0 = cfg.time
-    assim_times = pd.date_range(time0, time0 + dt.timedelta(hours=1), freq=timedelta_btw_assim)
-    last_assim_time = assim_times[-1]
+    update_vars = ['U', 'V', 'W', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'QSNOW', 'PSFC'],
+    #input_profile = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/2008-07-30_08:00/1/input_sounding',
+    verify_against = 'nat_250m_blockavg2km',
+    **cluster_defaults)
 
-    # loop over assimilations
-    for i, t in enumerate(assim_times):
-        
-        if i == 0:
-            if t == dt.datetime(2008, 7, 30, 11):
-                prior_init_time = dt.datetime(2008, 7, 30, 8)
-            else:
-                prior_init_time = t - dt.timedelta(minutes=15)
-                
-            cfg.update(time = t,
-                prior_init_time = prior_init_time,
-                prior_valid_time = t,
-                prior_path_exp = '/jetfs/home/lkugler/data/sim_archive/exp_nat250m_noDA/',)
-        else:
-            cfg.update(time = t,
-                prior_init_time = assim_times[i-1],
-                prior_valid_time = t,
-                prior_path_exp = cfg.dir_archive,)
-                    
-                    
-        id = w.assimilate(cfg, depends_on=id)
-
-        # 1) Set posterior = prior
-        id = w.prepare_IC_from_prior(cfg, depends_on=id)
-
-        # 2) Update posterior += updates from assimilation
-        id = w.update_IC_from_DA(cfg, depends_on=id)
-
-        if t == last_assim_time:
-            restart_interval = 9999
-            timedelta_integrate = dt.timedelta(hours=4)
+
+w = WorkFlows(cfg)
+w.prepare_WRFrundir(cfg)
+#id = w.run_ideal(cfg, depends_on=id)
+
+# assimilate at these times
+assim_times = pd.date_range(time0, time0 + dt.timedelta(hours=4), freq=timedelta_btw_assim)
+last_assim_time = assim_times[-1]
+
+# loop over assimilations
+for i, t in enumerate(assim_times):
+    
+    if i == 0:
+        if t == dt.datetime(2008, 7, 30, 11):
+            prior_init_time = dt.datetime(2008, 7, 30, 8)
         else:
-            restart_interval = timedelta_btw_assim.total_seconds()/60  # in minutes
-            timedelta_integrate = dt.timedelta(minutes=15)
+            prior_init_time = t - dt.timedelta(minutes=15)
+            
+        cfg.update(time = t,
+            prior_init_time = prior_init_time,
+            prior_valid_time = t,
+            prior_path_exp = '/jetfs/home/lkugler/data/sim_archive/exp_nat250m_noDA/',)
+    else:
+        cfg.update(time = t,
+            prior_init_time = assim_times[i-1],
+            prior_valid_time = t,
+            prior_path_exp = cfg.dir_archive,)
+                
+
+    id = w.assimilate(cfg, depends_on=id)
+
+    # 1) Set posterior = prior
+    id = w.prepare_IC_from_prior(cfg, depends_on=id)
+
+    # 2) Update posterior += updates from assimilation
+    id = w.update_IC_from_DA(cfg, depends_on=id)
+    
+    if t != last_assim_time:
+        # integrate until next assimilation
+        timedelta_integrate = dt.timedelta(minutes=15)
+        restart_interval = timedelta_btw_assim.total_seconds()/60  # in minutes
 
         cfg.update( WRF_start=t, 
                     WRF_end=t+timedelta_integrate, 
@@ -98,4 +113,24 @@ for t0 in assim_times_start:
         )
 
         # 3) Run WRF ensemble
-        id = w.run_WRF(cfg, depends_on=id)
+        id = w.run_WRF(cfg, depends_on=id)
\ No newline at end of file
+        id = w.run_RTTOV(cfg, depends_on=id)
+    
+    if t.minute == 0 and i != 0:
+        # full hour but not first one
+        # make long forecasts without restart files
+        timedelta_integrate = dt.timedelta(hours=4)
+        restart_interval = 9999
+        
+        cfg.update( WRF_start=t, 
+                    WRF_end=t+timedelta_integrate, 
+                    restart=True, 
+                    restart_interval=restart_interval,
+                    hist_interval_s=300,
+        )
+        # only verification depends on the output
+        # but we only have one run directory, so we need to wait until the other forecast is done
+        id = w.run_WRF(cfg, depends_on=id)
+        id = w.run_RTTOV(cfg, depends_on=id)
+    
+w.verify(cfg, depends_on=id)
\ No newline at end of file

tests/test_assim.py

@@ -32,7 +32,7 @@ def test_dartwrf_cycled_da():
         geo_em_forecast = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
         time = dt.datetime(2008, 7, 30, 11),
         
-        use_existing_obsseq = False,
+        assimilate_existing_obsseq = False,
         nature_wrfout_pattern = '/jetfs/home/lkugler/data/sim_archive/nat_250m_1600x1600x100/*/1/wrfout_d01_%Y-%m-%d_%H_%M_%S',
         #geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.2km_200x200',
         geo_em_nature = '/jetfs/home/lkugler/data/sim_archive/geo_em.d01.nc.250m_1600x1600',