diff --git a/config/cfg.py b/config/cfg.py index 0636712223129221510872e7ce30b960dcd371e9..b643e1bf2670580b399293f6d8f8f2e2286ed4f2 100755 --- a/config/cfg.py +++ b/config/cfg.py @@ -9,7 +9,7 @@ class ExperimentConfiguration(object): exp = ExperimentConfiguration() -exp.expname = "exp_v1.17_Pwbub-1_WV_obs20_loc5" +exp.expname = "exp_v1.17_P1-1_WV_10z" exp.model_dx = 2000 exp.n_ens = 40 exp.n_nodes = 10 @@ -17,19 +17,19 @@ exp.n_nodes = 10 # 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) +exp.cov_loc_vert_km_horiz_km = (2, 20) -n_obs = 256 #121: 30km, 256:16x16 (20km); 961: 10km resoltn # radar: n_obs for each observation height level +n_obs = 961 #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=7.5) + cov_loc_radius_km=20) 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=5) + error_generate=1., error_assimilate=False, + cov_loc_radius_km=20) ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]', kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs, @@ -39,8 +39,8 @@ 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=[7000,], #np.arange(1000, 15001, 1000), - cov_loc_radius_km=5) + heights=np.arange(2000, 14001, 2000), + cov_loc_radius_km=20) t2m = dict(plotname='SYNOP Temperature', plotunits='[K]', kind='SYNOP_TEMPERATURE', n_obs=n_obs, @@ -53,9 +53,9 @@ psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]', cov_loc_radius_km=32) -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'] +exp.observations = [wv73] #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'] # directory paths depend on the name of the experiment cluster.expname = exp.expname diff --git a/config/clusters.py b/config/clusters.py index a51c84b49f8c5e14f2398a9b7537e95654842b4d..8d2709972303374e701dcfcd5238407532ce882f 100755 --- a/config/clusters.py +++ b/config/clusters.py @@ -25,6 +25,7 @@ class ClusterConfig(object): vsc = ClusterConfig() vsc.name = 'vsc' vsc.python = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/python' +vsc.python_enstools = '/home/fs71386/lkugler/miniconda3/envs/enstools/bin/python' vsc.ncks = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/ncks' vsc.tmpfiledir = '/gpfs/data/fs71386/lkugler' vsc.userdir = '/home/fs71386/lkugler' @@ -34,10 +35,9 @@ vsc.dart_srcdir = '/home/fs71386/lkugler/DART/DART-9.11.9/models/wrf/work' vsc.dartrundir = '/gpfs/data/fs71386/lkugler/run_DART' vsc.scriptsdir = '/home/fs71386/lkugler/DART-WRF/scripts/' -vsc.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.16_Pwbub_nat/2008-07-30_09:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' -#vsc.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/wrfout_d01_%Y-%m-%d_%H:%M:%S' +vsc.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.17_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' #vsc.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/from_LMU/raso.nat.<iens>.wrfprof' -vsc.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.nat.001.wrfprof' +vsc.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.fc.<iens>.wrfprof' vsc.ideal = vsc.userdir+'/compile/bin/ideal-v4.2.2_v1.16.exe' vsc.wrfexe = vsc.userdir+'/compile/bin/wrf-v4.3_v1.16.exe' @@ -45,7 +45,7 @@ vsc.namelist = vsc.scriptsdir+'/../templates/namelist.input' vsc.run_WRF = '/home/fs71386/lkugler/DART-WRF/scripts/run_ens.vsc.sh' vsc.slurm_cfg = {"account": "p71386", "partition": "mem_0384", "qos": "p71386_0384", - "ntasks": "1", "nodes": "1", "ntasks-per-node": "48", "ntasks-per-core": "1", + "nodes": "1", "ntasks": "1", "ntasks-per-node": "48", "ntasks-per-core": "1", "mail-type": "FAIL", "mail-user": "lukas.kugler@univie.ac.at"} diff --git a/generate_nature.py b/generate_nature.py index ca9db9a68124f944fe44115aeaa0187d4eab74d1..bf942d0b69fc1b1aa1ed3414b177ccc335eba803 100755 --- a/generate_nature.py +++ b/generate_nature.py @@ -29,32 +29,22 @@ print('starting osse') backup_scripts() id = None -init_time = dt.datetime(2008, 7, 30, 9) +init_time = dt.datetime(2008, 7, 30, 6) id = prepare_wrfinput(init_time) # create initial conditions -# get initial conditions from archive -integration_end_time = dt.datetime(2008, 7, 30, 12) -#exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.16_P1_40mem' -#id = update_wrfinput_from_archive(integration_end_time, init_time, exppath_arch, depends_on=id) +if False: + init_time = dt.datetime(2008, 7, 30, 8) + # get initial conditions from archive + integration_end_time = dt.datetime(2008, 7, 30, 9) + exppath_arch = cluster.archivedir #'/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.16_P1_40mem' + id = update_wrfinput_from_archive(integration_end_time, init_time, exppath_arch, depends_on=id) + + #id = wrfinput_insert_wbubble(depends_on=id) + -id = wrfinput_insert_wbubble(depends_on=id) +begin = dt.datetime(2008, 7, 30, 6) +end = dt.datetime(2008, 7, 31, 0) -begin = dt.datetime(2008, 7, 30, 9, 0) -end = dt.datetime(2008, 7, 30, 12, 0) - -# whole forecast timespan -hist_interval = 5 -radt = 5 -s = my_Slurm("namelist", cfg_update=dict(time="2")) -id = s.run(' '.join([cluster.python, - cluster.scriptsdir+'/prepare_namelist.py', - begin.strftime('%Y-%m-%d_%H:%M'), - end.strftime('%Y-%m-%d_%H:%M'), - str(hist_interval), str(radt),]), - depends_on=[id]) - -s = my_Slurm("EnsWRF", cfg_update={"nodes": "1", "array": "1-"+str(exp.n_nodes), - "mem-per-cpu": "2G", "mail-type": "BEGIN,FAIL,END"}) -cmd = script_to_str(cluster.run_WRF).replace('<expname>', exp.expname) -id = s.run(cmd, depends_on=[id]) +id = run_ENS(begin=begin, end=end, first_minute=False, depends_on=id) +create_satimages(begin, depends_on=id) diff --git a/scheduler.py b/scheduler.py index 7c70e6e3cca5c274468f95209d7a1d981f70f548..b0011fdffc882a512724a113dd88f779309934a5 100755 --- a/scheduler.py +++ b/scheduler.py @@ -207,7 +207,7 @@ def assimilate(assim_time, prior_init_time, def create_satimages(init_time, depends_on=None): - s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "30"}) + s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "30", "nodes": "1"}) s.run(cluster.python+' /home/fs71386/lkugler/RTTOV-WRF/run_init.py '+cluster.archivedir +init_time.strftime('/%Y-%m-%d_%H:%M/'), depends_on=[depends_on]) @@ -219,20 +219,31 @@ def mailme(depends_on=None): def gen_obsseq(depends_on=None): s = my_Slurm("obsseq_netcdf", cfg_update={"time": "10", "mail-type": "FAIL,END"}) - s.run(cluster.python+' '+cluster.scripts_rundir+'/obsseq_to_netcdf.py', + id = s.run(cluster.python+' '+cluster.scripts_rundir+'/obsseq_to_netcdf.py', + depends_on=[depends_on]) + return id + +def verify(depends_on=None): + s = my_Slurm("verify", cfg_update={"time": "180", "mail-type": "FAIL,END", + "ntasks": "96", "ntasks-per-node": "96", "ntasks-per-core": "2"}) + s.run(cluster.python_enstools+' '+cluster.userdir+'/osse_analysis/analyze_fc.py '+exp.expname+' has_node', depends_on=[depends_on]) def copy_to_jet(depends_on=None): - s = my_Slurm("rsync-jet", cfg_update={"time": "1", "mail-type": "FAIL"}) - s.run('bash rsync -avh '+cluster.archivedir+' lkugler@jet01.img.univie.ac.at:/jetfs/home/lkugler/data_jetfs/sim_archive/ &', + Slurm('rsync-jet', slurm_kwargs={"time": "30", + "account": "p71386", "partition": "mem_0384", "qos": "p71386_0384", + "ntasks": "1", "mem": "5gb", + "mail-type": "FAIL", "mail-user": "lukas.kugler@univie.ac.at"}, + log_dir=log_dir, scripts_dir=slurm_scripts_dir, + ).run("bash -c 'nohup rsync -avh "+cluster.archivedir+" lkugler@jet01.img.univie.ac.at:/jetfs/home/lkugler/data/sim_archive/ &'", depends_on=[depends_on]) ################################ if __name__ == "__main__": print('starting osse') - timedelta_integrate = dt.timedelta(minutes=10) - timedelta_btw_assim = dt.timedelta(minutes=10) + timedelta_integrate = dt.timedelta(minutes=15) + timedelta_btw_assim = dt.timedelta(minutes=15) backup_scripts() id = None @@ -257,10 +268,10 @@ if __name__ == "__main__": id = prepare_wrfinput(init_time) # create initial conditions # get initial conditions from archive - integration_end_time = dt.datetime(2008, 7, 30, 13) + integration_end_time = dt.datetime(2008, 7, 30, 10) exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.16_P1_40mem' - #id = update_wrfinput_from_archive(integration_end_time, init_time, exppath_arch, depends_on=id) - id = wrfinput_insert_wbubble(depends_on=id) + id = update_wrfinput_from_archive(integration_end_time, init_time, exppath_arch, depends_on=id) + #id = wrfinput_insert_wbubble(depends_on=id) prior_path_exp = exppath_arch # for next assimilation # values for assimilation @@ -268,7 +279,7 @@ if __name__ == "__main__": assim_time = integration_end_time prior_init_time = init_time - while time <= dt.datetime(2008, 7, 30, 14): + while time <= dt.datetime(2008, 7, 30, 10, 30): id = assimilate(assim_time, prior_init_time, @@ -280,8 +291,8 @@ if __name__ == "__main__": this_forecast_init = assim_time # start integration from here timedelta_integrate = timedelta_btw_assim - if this_forecast_init.minute == 0 and this_forecast_init.hour != dt.datetime(2008, 7, 30, 13): # longer forecast every full hour - timedelta_integrate = dt.timedelta(minutes=3*60) + if this_forecast_init.minute in [30, 0]: # longer forecast every full hour + timedelta_integrate = dt.timedelta(hours=2) this_forecast_end = assim_time + timedelta_integrate @@ -298,5 +309,6 @@ if __name__ == "__main__": assim_time = time prior_init_time = time - timedelta_btw_assim - gen_obsseq(id) # mailme(id) - copy_to_jet(id) + id = gen_obsseq(id) + verify(id) + # copy to jet diff --git a/scripts/apply_obs_op_dart.py b/scripts/apply_obs_op_dart.py index 4563598a47d544826bbfc2b69c11c49a9c929a0c..87f028ed1ff435db58a3f42182fe98e610f2df21 100755 --- a/scripts/apply_obs_op_dart.py +++ b/scripts/apply_obs_op_dart.py @@ -33,7 +33,7 @@ if __name__ == '__main__': savedir = cluster.archivedir+'/obs_seq_final_1min/' aso.set_DART_nml(just_prior_values=True) - osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=False) + osq.create_obsseqin_alltypes(time, exp.observations, 0.) # prepare dummy nature (this Hx is irrelevant) os.chdir(cluster.dartrundir) diff --git a/scripts/assim_synth_obs.py b/scripts/assim_synth_obs.py index e39ad55526388ecde8d7e634984f1fd37bcc65eb..a4a93114ab88de996cb84646c84f497b9edf0d52 100755 --- a/scripts/assim_synth_obs.py +++ b/scripts/assim_synth_obs.py @@ -97,7 +97,8 @@ 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 ', float(previous_err_var.strip()), 'new error', new_err_obs_i) + if debug: + print(line.strip(), '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 @@ -411,12 +412,9 @@ if __name__ == "__main__": n_obs_z = len(obscfg.get('heights', [1,])) n_obs_3d = n_obs * n_obs_z - parametrized = obscfg.get('sat_channel') == 6 - - if not parametrized: - err_assim = np.zeros(n_obs_3d) + obscfg['error_assimilate'] - - else: # error parametrization for WV73 + do_error_parametrization = ((obscfg['error_assimilate'] == False) and (obscfg.get('sat_channel') == 6)) + + if do_error_parametrization: # get observations for sat 6 osq.create_obsseqin_alltypes(time, [obscfg,], np.zeros(n_obs_3d)) run_perfect_model_obs() @@ -426,6 +424,8 @@ if __name__ == "__main__": Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory err_assim = calc_obserr_WV73(Hx_nat, Hx_prior) + else: + err_assim = np.zeros(n_obs_3d) + obscfg['error_assimilate'] error_assimilate.extend(err_assim) # the obs-error we assume for assimilating observations @@ -437,7 +437,7 @@ if __name__ == "__main__": for i, obscfg in enumerate(exp.observations): error_generate.extend(np.zeros(n_obs_3d) + obscfg['error_generate']) - osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=error_generate) + osq.create_obsseqin_alltypes(time, exp.observations, error_generate) set_DART_nml() diff --git a/scripts/create_obsseq.py b/scripts/create_obsseq.py index 1d06ac70aae766dab48f5863648e553c370e8a94..754de0e344099ab811d084a254d8e0a1893370c1 100755 --- a/scripts/create_obsseq.py +++ b/scripts/create_obsseq.py @@ -140,6 +140,7 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True, n_obs_x = int(np.sqrt(n_obs)) nx = len(ds.south_north) # number of gridpoints in one direction model_dx_km = exp.model_dx/1000 + print('assuming', model_dx_km, 'km model grid') omit_covloc_radius_on_boundary = True if omit_covloc_radius_on_boundary: # in order to avoid an increment step on the boundary @@ -382,8 +383,7 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors, archive_obs_coord Args: time_dt (dt.datetime): time of observation list_obscfg (list of dict) - obs_errors (list of float, False): contains observation errors, one for each observation - if False: use zero error + obs_errors (np.array): contains observation errors, one for each observation archive_obs_coords (bool, str): False or str (filepath where `obs_seq.in` will be saved) """ print('creating obs_seq.in:') @@ -411,13 +411,7 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors, archive_obs_coord coords = append_hgt_to_coords(coords, vert_coords) n_obs_3d_thistype = len(coords) - # define obs error - if obs_errors == False: - obs_errors = np.zeros(n_obs_3d_thistype) - assert len(obs_errors) == n_obs_3d_thistype, 'len(obs_errors) == n_obs_3d_thistype' - obserr_std = obs_errors #np.zeros(n_obs_3d_thistype) - #if obs_errors: - # obserr_std += obscfg[obs_errors] + obserr_std = np.zeros(n_obs_3d_thistype) + obs_errors sat_info = write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt) @@ -455,7 +449,7 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors, archive_obs_coord if __name__ == '__main__': # for testing time_dt = dt.datetime(2008, 7, 30, 9, 0) - n_obs = 900 # radar: n_obs for each observation height level + n_obs = 22801 # 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, @@ -490,10 +484,11 @@ if __name__ == '__main__': #create_obsseq_in(time_dt, radar, archive_obs_coords=False) #'./coords_stage1.pkl') - create_obsseqin_alltypes(time_dt, [vis, wv73], obs_errors='error_generate', archive_obs_coords='./obs_coords.pkl') + create_obsseqin_alltypes(time_dt, [wv73], obs_errors=False, archive_obs_coords=False) #'./obs_coords.pkl') - error_assimilate = 5.*np.ones(n_obs*len(radar['heights'])) - import assim_synth_obs as aso - #aso.replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate) + if False: + error_assimilate = 5.*np.ones(n_obs*len(radar['heights'])) + import assim_synth_obs as aso + aso.replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate) diff --git a/templates/input.eval.nml b/templates/input.eval.nml index df66d8ea49826c9a97038fead9e9a0af10450b37..92212cad31ee54584af062f76d80425f4c47fe1a 100644 --- a/templates/input.eval.nml +++ b/templates/input.eval.nml @@ -174,12 +174,12 @@ # &obs_def_radar_mod_nml - apply_ref_limit_to_obs = .false., - reflectivity_limit_obs = -10.0, - lowest_reflectivity_obs = -10.0, - apply_ref_limit_to_fwd_op = .false., - reflectivity_limit_fwd_op = -10.0, - lowest_reflectivity_fwd_op = -10.0, + apply_ref_limit_to_obs = .true., + reflectivity_limit_obs = 5.0, + lowest_reflectivity_obs = 5.0, + apply_ref_limit_to_fwd_op = .true., + reflectivity_limit_fwd_op = 5.0, + lowest_reflectivity_fwd_op = 5.0, max_radial_vel_obs = 1000000, allow_wet_graupel = .false., microphysics_type = 5 , diff --git a/templates/input.nml b/templates/input.nml index 356d7258e25220b20d55fb62530cb22e2690bacc..1281697e1998c3501b206f03f3ea7ee3cf32a481 100644 --- a/templates/input.nml +++ b/templates/input.nml @@ -176,12 +176,12 @@ # &obs_def_radar_mod_nml - apply_ref_limit_to_obs = .false., - reflectivity_limit_obs = -10.0, - lowest_reflectivity_obs = -10.0, - apply_ref_limit_to_fwd_op = .false., - reflectivity_limit_fwd_op = -10.0, - lowest_reflectivity_fwd_op = -10.0, + apply_ref_limit_to_obs = .true., + reflectivity_limit_obs = 5.0, + lowest_reflectivity_obs = 5.0, + apply_ref_limit_to_fwd_op = .true., + reflectivity_limit_fwd_op = 5.0, + lowest_reflectivity_fwd_op = 5.0, max_radial_vel_obs = 1000000, allow_wet_graupel = .false., microphysics_type = 5 , diff --git a/templates/input.prioronly.nml b/templates/input.prioronly.nml index 19471bb280abc0e238bc62657d674c5338aa5d10..994e6e33cf499d236705e3e2dee8c57dd7c17d59 100644 --- a/templates/input.prioronly.nml +++ b/templates/input.prioronly.nml @@ -170,12 +170,12 @@ # &obs_def_radar_mod_nml - apply_ref_limit_to_obs = .false., - reflectivity_limit_obs = -10.0, - lowest_reflectivity_obs = -10.0, - apply_ref_limit_to_fwd_op = .false., - reflectivity_limit_fwd_op = -10.0, - lowest_reflectivity_fwd_op = -10.0, + apply_ref_limit_to_obs = .true., + reflectivity_limit_obs = 5.0, + lowest_reflectivity_obs = 5.0, + apply_ref_limit_to_fwd_op = .true., + reflectivity_limit_fwd_op = 5.0, + lowest_reflectivity_fwd_op = 5.0, max_radial_vel_obs = 1000000, allow_wet_graupel = .false., microphysics_type = 5 , diff --git a/templates/namelist.input b/templates/namelist.input index b488b18c87a33df6d74b41939a2ebdd30ce1004e..78cabae47630847d7572dc18c453ad134eef5cc2 100644 --- a/templates/namelist.input +++ b/templates/namelist.input @@ -11,6 +11,7 @@ end_hour = <HH2>, 00, 00, end_minute = <MM2>, 120, 120, end_second = 0, 00, 00, + interval_seconds = 86400 history_interval = <hist_interval>, 15, 15, frames_per_outfile = 1, 1, 1, history_outname = '<archivedir>/wrfout_d<domain>_<date>' @@ -30,7 +31,6 @@ starting_time_step = 8, max_time_step = 16, min_time_step = 6, - step_to_output_time = .true., time_step = 8, time_step_fract_num = 0, time_step_fract_den = 1,