vert cov loc option

config/cfg.py

@@ -20,7 +20,7 @@ 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,
            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_radius_km=10, cov_loc_vert_km=5)
 
 exp.observations = [vis, ]
 

scripts/create_obsseq.py

@@ -104,11 +104,27 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True,
     assert len(coords) == n_obs, (len(coords), n_obs)
     return coords
 
+
 def write_generic_obsseq(obs_name, obs_kind_nr, error_var, coords,
-                         dart_date_day, secs_thatday, output_path):
+                         dart_date_day, secs_thatday, output_path,
+                         vert_coord_sfc=False):
+    """
+    Args:
+        dart_date_day (str): DART internal time formatted date
+        secs_thatday (str): DART internal time of day (seconds since 0 UTC)
+        vert_coord_sfc (bool):
+            if True, then vertical coordinate is height above ground, i.e. "surface observation"
+            if False, then vertical is hgt_AMSL
+    """
+
+    vert_coord_sys = 3  # meters AMSL
+    if vert_coord_sfc:
+        vert_coord_sys = -1
+
     n_obs_str = str(int(n_obs))
     error_var = str(error_var)
     line_obstypedef = obs_kind_nr+' '+obs_name
+    vert_coord_sys = str(vert_coord_sys)
 
     msg = """
  obs_sequence
@@ -135,7 +151,7 @@ obs_kind_definitions
           -1           """+str(i_obs+1)+"""          -1
 obdef
 loc3d
-     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
+     """+lon_rad+"        "+lat_rad+"        "+hgt_m+"     "+vert_coord_sys+"""
 kind
          """+obs_kind_nr+"""
  """+secs_thatday+"""     """+dart_date_day+"""
@@ -147,7 +163,7 @@ kind
           """+str(i_obs-1)+"""           -1          -1
 obdef
 loc3d
-     """+lon_rad+"""        """+lat_rad+"""        """+hgt_m+"""     3
+     """+lon_rad+"        "+lat_rad+"        "+hgt_m+"     "+vert_coord_sys+"""
 kind
          """+obs_kind_nr+"""
  """+secs_thatday+"""     """+dart_date_day+"""
@@ -163,6 +179,7 @@ kind
         f.write(msg)
         print(fpath, 'saved.')
 
+
 def sat(time_dt, channel_id, coords, error_var, output_path='./'):
     """Create obs_seq.in
 

scripts/gen_synth_obs.py

@@ -72,17 +72,31 @@ def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
         for line in obsseq_new:
             f.write(line)
 
+
 def set_input_nml(sat_channel=False, just_prior_values=False,
-                  cov_loc_radius_km=32):
+                  cov_loc_radius_km=32, cov_loc_vert_km=False):
     """descr"""
+    cov_loc_radian = cov_loc_radius_km/earth_radius_km
+    
     if just_prior_values:
         template = cluster.scriptsdir+'/../templates/input.prioronly.nml'
     else:
         template = cluster.scriptsdir+'/../templates/input.nml'
     copy(template, cluster.dartrundir+'/input.nml')
-    sed_inplace(cluster.dartrundir+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
-    cov_loc_radian = cov_loc_radius_km/earth_radius_km
-    sed_inplace(cluster.dartrundir+'/input.nml', '<cov_loc_radian>', str(cov_loc_radian))
+
+    options = {'<n_ens>': str(int(exp.n_ens)),
+               '<cov_loc_radian>': str(cov_loc_radian)}
+
+    if cov_loc_vert_km:
+        cov_loc_vert_rad = cov_loc_vert_km/cov_loc_radian
+        options['<horiz_dist_only>'] = '.false.'
+        options['<vert_norm_hgt>'] = str(cov_loc_vert_rad)
+    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)
 
     # input.nml for RTTOV
     if sat_channel > 0:
@@ -109,6 +123,7 @@ if __name__ == "__main__":
         sat_channel = obscfg.get('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)
 
         # generate obs_seq.in
         obs_coords = osq.calc_obs_locations(n_obs, coords_from_domaincenter=False, 
@@ -126,7 +141,9 @@ if __name__ == "__main__":
             raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
 
         # generate observations (obs_seq.out)
-        set_input_nml(sat_channel=sat_channel, cov_loc_radius_km=cov_loc)
+        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')
@@ -136,7 +153,8 @@ if __name__ == "__main__":
         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)
+                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())
@@ -168,7 +186,8 @@ if __name__ == "__main__":
 
                 # correct input.nml for actual assimilation later on
                 set_input_nml(sat_channel=sat_channel,
-                              cov_loc_radius_km=cov_loc)
+                              cov_loc_radius_km=cov_loc,
+                              cov_loc_vert_km=cov_loc_vert_km)
 
         # rename according to i_obs
         os.rename(cluster.dartrundir+'/obs_seq.out', 

templates/input.nml

@@ -159,7 +159,8 @@
 				'RADAR_REFLECTIVITY',
                                 'MSG_4_SEVIRI_RADIANCE',
 				'MSG_4_SEVIRI_TB',
-				'MSG_4_SEVIRI_BDRF'
+				'MSG_4_SEVIRI_BDRF',
+                                'LAND_SFC_PRESSURE'
    evaluate_these_obs_types   = 'RADIOSONDE_SPECIFIC_HUMIDITY',
    /
 
@@ -244,9 +245,9 @@
 # an internal structure that speeds up searches.  don't change it
 # based on your grid size.  nlon must be an odd number.
 &location_nml
-   horiz_dist_only             = .true.,
+   horiz_dist_only             = <horiz_dist_only>,
    vert_normalization_pressure = 6666666.7,
-   vert_normalization_height   = 5000000.0,
+   vert_normalization_height   = <vert_norm_hgt>,
    vert_normalization_level    = 2666.7,
    vert_normalization_scale_height = 10.0,
    approximate_distance        = .false.,