different cutoffs for each type

config/cfg.py

@@ -9,22 +9,27 @@ class ExperimentConfiguration(object):
 
 
 exp = ExperimentConfiguration()
-exp.expname = "exp_v1.16_Pwbub-1_40mem"
+exp.expname = "exp_v1.17_Pwbub-1_WV_obs20_loc5"
 exp.model_dx = 2000
 exp.n_ens = 40
 exp.n_nodes = 10
 
-n_obs = 121  #961 900: 10km resoltn # radar: n_obs for each observation height level
+# localize vertically, if it has a vertical position
+# needs a horizontal scale too, to calculate the vertical normalization
+# since you can not specify different vertical localizations for diff. variables
+exp.cov_loc_vert_km_horiz_km = (3, 5)  
+
+n_obs = 256  #121: 30km, 256:16x16 (20km); 961: 10km resoltn # radar: n_obs for each observation height level
 
 vis = dict(plotname='VIS 0.6µm',  plotunits='[1]',
            kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs, 
            error_generate=0.03, error_assimilate=0.06,
-           cov_loc_radius_km=32)
+           cov_loc_radius_km=7.5)
 
 wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
             kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs, 
-            error_generate=1., error_assimilate=False,
-            cov_loc_radius_km=32)
+            error_generate=1., error_assimilate=False, # 
+            cov_loc_radius_km=5)
 
 ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
              kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs, 
@@ -34,21 +39,21 @@ ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
 radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
              kind='RADAR_REFLECTIVITY', n_obs=n_obs, 
              error_generate=2.5, error_assimilate=5.,
-             heights=np.arange(1000, 15001, 1000),
-             cov_loc_radius_km=32, cov_loc_vert_km=4)
+             heights=[7000,], #np.arange(1000, 15001, 1000),
+             cov_loc_radius_km=5)
 
 t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
            kind='SYNOP_TEMPERATURE', n_obs=n_obs, 
            error_generate=0.1, error_assimilate=1.,
-           cov_loc_radius_km=20, cov_loc_vert_km=3)
+           cov_loc_radius_km=20)
 
 psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
             kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, 
             error_generate=50., error_assimilate=100.,
-            cov_loc_radius_km=32, cov_loc_vert_km=5)
+            cov_loc_radius_km=32)
 
 
-exp.observations = [] #wv73, vis] # 108, wv73, vis]
+exp.observations = [radar] # 108, wv73, vis]
 exp.update_vars = ['T', 'QVAPOR', 'QCLOUD', 'QICE','CLDFRA']
 #exp.update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'TSK', 'CLDFRA']
 

scripts/apply_obs_op_dart.py

@@ -48,37 +48,34 @@ if __name__ == '__main__':
 
     # only the prior state values are of interest in this file
     # observation and truth is wrong in this file (dummy)
-    istage = 0
-    archive_stage = savedir+'/assim_stage'+str(istage)
-    aso.archive_diagnostics(archive_stage, time)
+    aso.archive_osq_final(time, posterior_1min=True)
 
-    sys.exit()  # multi stage below
+    # sys.exit()  # multi stage below
+    # n_stages = len(exp.observations)
+    # for istage, obscfg in enumerate(exp.observations):
 
-    n_stages = len(exp.observations)
-    for istage, obscfg in enumerate(exp.observations):
+    #     n_obs = obscfg['n_obs']
+    #     sat_channel = obscfg.get('sat_channel', False)
+    #     obscfg['folder_obs_coords'] = False
 
-        n_obs = obscfg['n_obs']
-        sat_channel = obscfg.get('sat_channel', False)
-        obscfg['folder_obs_coords'] = False
+    #     aso.set_DART_nml(cov_loc_radius_km=obscfg['cov_loc_radius_km'],
+    #                      cov_loc_vert_km=obscfg.get('cov_loc_vert_km', False), 
+    #                      just_prior_values=True)
 
-        aso.set_DART_nml(cov_loc_radius_km=obscfg['cov_loc_radius_km'],
-                         cov_loc_vert_km=obscfg.get('cov_loc_vert_km', False), 
-                         just_prior_values=True)
+    #     osq.create_obsseq_in(time, obscfg)
 
-        osq.create_obsseq_in(time, obscfg)
 
 
+    #     # prepare dummy nature (this Hx is irrelevant)
+    #     os.chdir(cluster.dartrundir)
+    #     os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
+    #     wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
+    #                        cluster.dartrundir+'/wrfout_d01')
+    #     aso.run_perfect_model_obs()
+    #     aso.assimilate(nproc=96)
 
-        # prepare dummy nature (this Hx is irrelevant)
-        os.chdir(cluster.dartrundir)
-        os.system('cp ./advance_temp1/wrfout_d01 ./wrfout_d01')
-        wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
-                           cluster.dartrundir+'/wrfout_d01')
-        aso.run_perfect_model_obs()
-        aso.assimilate(nproc=96)
-
-        # only the prior state values are of interest in this file
-        # observation and truth is wrong in this file (dummy)
-        archive_stage = savedir+'/assim_stage'+str(istage)
-        aso.archive_diagnostics(archive_stage, time)
+    #     # only the prior state values are of interest in this file
+    #     # observation and truth is wrong in this file (dummy)
+    #     archive_stage = savedir+'/assim_stage'+str(istage)
+    #     aso.archive_diagnostics(archive_stage, time)
 

scripts/assim_synth_obs.py

@@ -97,7 +97,7 @@ def replace_errors_obsseqout(f, new_errors):
 
             previous_err_var = obsseq[line_error_obs_i]
             new_err_obs_i = new_errors[i_obs]**2  # variance in obs_seq.out
-            if debug: print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
+            if debug: print('previous err var ', float(previous_err_var.strip()), 'new error', new_err_obs_i)
             obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n'
 
             i_obs += 1  # next iteration
@@ -110,6 +110,7 @@ def replace_errors_obsseqout(f, new_errors):
 
 def set_DART_nml_singleobstype(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
                  just_prior_values=False):
+    """this function is outdated and will not work"""
     cov_loc_radian = cov_loc_radius_km/earth_radius_km
     
     if just_prior_values:
@@ -145,8 +146,19 @@ def set_DART_nml_singleobstype(sat_channel=False, cov_loc_radius_km=32, cov_loc_
         rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
         append_file(cluster.dartrundir+'/input.nml', rttov_nml)
 
-def set_DART_nml(cov_loc_radius_km=32, cov_loc_vert_km=False, just_prior_values=False):
-    cov_loc_radian = cov_loc_radius_km/earth_radius_km
+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:
         template = cluster.scriptsdir+'/../templates/input.eval.nml'
@@ -156,15 +168,20 @@ def set_DART_nml(cov_loc_radius_km=32, cov_loc_vert_km=False, just_prior_values=
 
     # options keys are replaced in input.nml with values
     options = {'<n_ens>': str(int(exp.n_ens)),
-               '<cov_loc_radian>': str(cov_loc_radian)}
+               '<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),
+               }
 
-    if cov_loc_vert_km:
-        vert_norm_rad = earth_radius_km*cov_loc_vert_km/cov_loc_radius_km*1000
-        options['<horiz_dist_only>'] = '.false.'
-        options['<vert_norm_hgt>'] = str(vert_norm_rad)
-    else:
+    # if cov_loc_vert_km:
     options['<horiz_dist_only>'] = '.true.'
-        options['<vert_norm_hgt>'] = '50000.0'  # dummy value
+    
+    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)
+    # else:
+    #     options['<horiz_dist_only>'] = '.true.'
+    #     options['<vert_norm_hgt>'] = '50000.0'  # dummy value
 
     for key, value in options.items():
         sed_inplace(cluster.dartrundir+'/input.nml', key, value)
@@ -295,42 +312,56 @@ def recycle_output():
 
 ############### archiving
 
-def archive_assimilation(time)
+def archive_osq_final(time, posterior_1min=False):
+    """Save obs_seq.final file for later. 
+    time (dt.datetime) : time of sampling values from files
+    posterior_1min (bool) : False if usual assimilation
+    """
 
-    print('archive obs space diagnostics')
+    if posterior_1min:
+        archive_dir = cluster.archivedir+'/obs_seq_final_1min/'
+    else:
         archive_dir = cluster.archivedir+'/obs_seq_final/'
     mkdir(archive_dir)
     fout = archive_dir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final')
     copy(cluster.dartrundir+'/obs_seq.final', fout)
     print(fout, 'saved.')
 
-    # try:  # what are regression diagnostics?!
-    #     print('archive regression diagnostics')
+    # try:  
+    #     print('archive regression diagnostics') # what are regression diagnostics?!
     #     copy(cluster.dartrundir+'/reg_diagnostics', archive_dir+'/reg_diagnostics')
     # except Exception as e:
     #     warnings.warn(str(e))
 
+def archive_filteroutput(time):
     print('archiving output')
-    archive_assim = cluster.archivedir + '/assim_stage0/'
+    archive_assim = cluster.archivedir + time.strftime('/%Y-%m-%d_%H:%M/assim_stage0/')
     mkdir(archive_assim)
     copy(cluster.dartrundir+'/input.nml', archive_assim+'/input.nml')
 
     for iens in range(1, exp.n_ens+1):  # single members
-        filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
-        copy(filter_out, archive_assim+'/filter_restart_d01.'+str(iens).zfill(4))
+        copy(cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4), 
+             archive_assim+'/filter_restart_d01.'+str(iens).zfill(4))
+
+    try:  # not necessary for next forecast run
+        for iens in range(1, exp.n_ens+1):
+            copy(cluster.dartrundir+'/postassim_member_'+str(iens).zfill(4)+'.nc', 
+                archive_assim+'/postassim_member_'+str(iens).zfill(4)+'.nc')
 
         for f in ['output_mean.nc', 'output_sd.nc']:  # copy mean and sd to archive
             copy(cluster.dartrundir+'/'+f, archive_assim+'/'+f)
 
-def archive_obs_generation(time):
+    except Exception as e:
+        warnings.warn(str(e))
 
+
+def archive_osq_out(time):
     dir_obsseq = cluster.archivedir+'/obs_seq_out/'
     os.makedirs(dir_obsseq, exist_ok=True)
     copy(cluster.dartrundir+'/obs_seq.out', dir_obsseq+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.out'))
 
 
 if __name__ == "__main__":
-
     """Assimilate observations (different obs types)
     as defined in config/cfg.py
     for a certain timestamp (argument) of the nature run (defined in config/clusters.py)
@@ -339,7 +370,7 @@ if __name__ == "__main__":
     for each assimilation stage (one obs_seq.in and e.g. one observation type):
     1) create obs_seq.in with obs-errors
     2) prepare nature run for DART
-    3) create obs from nature (obs_seq.out)
+    3) create obs from nature (obs_seq.out) using pre-defined obs-errors
     4) Assimilate
       - adapt obs errors for assimilation
         - calculate assim obs error from parametrization
@@ -408,12 +439,10 @@ if __name__ == "__main__":
 
     osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=error_generate)
 
-    first_obstype = exp.observations[0]  # TODO: different for each observation type
-    set_DART_nml(cov_loc_radius_km=first_obstype['cov_loc_radius_km'],
-                 cov_loc_vert_km=first_obstype.get('cov_loc_vert_km', False))
+    set_DART_nml()
 
     run_perfect_model_obs()  # actually create observations that are used to assimilate
-    archive_obs_generation(time)
+    archive_osq_out(time)
 
     ################################################
     print(' 3) assimilate with observation-errors for assimilation')
@@ -423,4 +452,5 @@ if __name__ == "__main__":
     assimilate()
     print('filter took', time_module.time()-t, 'seconds')
     
-    archive_assimilation(time)
+    archive_filteroutput(time)
+    archive_osq_final(time)

templates/input.nml

@@ -51,7 +51,7 @@
    output_timestamps        = .false.,
    trace_execution          = .false.,
 
-   stages_to_write          = 'output'
+   stages_to_write          = 'postassim', 'output'
    output_members           = .true.
    output_mean              = .true.
    output_sd                = .true.
@@ -114,6 +114,8 @@
    convert_all_state_verticals_first = .true.
    convert_all_obs_verticals_first   = .true.
    print_every_nth_obs               = 0,
+   special_localization_obs_types  = <list_obstypes>,
+   special_localization_cutoffs    = <list_cutoffs>
    /
 
 &cov_cutoff_nml