diff --git a/config/cfg.py b/config/cfg.py index 9bb0369803d1f40bc64889090102e95319ff073c..ab01c007582c717ee32ef221128cbc0f9d199c5b 100755 --- a/config/cfg.py +++ b/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, ] diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py index 91b39dff9328115ca2cf8d2e2141eafe064158e8..d8035347ff696da9869fc71de4465e594994747a 100755 --- a/scripts/create_obsseq.py +++ b/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 diff --git a/scripts/gen_synth_obs.py b/scripts/gen_synth_obs.py index ec21526354673e6e61f9e6bb91eee83dfc53dc23..c40a7c73b69bbc7f9c8ffb0c6f7ef198496888b0 100755 --- a/scripts/gen_synth_obs.py +++ b/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', diff --git a/templates/input.nml b/templates/input.nml index a016646b7447825c68a5835d437e74c52a48a2b5..2e593ef49d260b916fee15b9d76bb3315e605882 100644 --- a/templates/input.nml +++ b/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.,