Skip to content
Snippets Groups Projects
Commit 68fc2efe authored by lkugler's avatar lkugler
Browse files

update

parent ad19784d
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 92 additions and 95 deletions
config/cfg.py
+ 11
− 11
View file @ 68fc2efe
......@@ -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
config/clusters.py
+ 4
− 4
View file @ 68fc2efe
......@@ -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"}
......
generate_nature.py
+ 14
− 24
View file @ 68fc2efe
......@@ -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
if False:
init_time = dt.datetime(2008, 7, 30, 8)
# 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)
id = wrfinput_insert_wbubble(depends_on=id)
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])
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)
begin = dt.datetime(2008, 7, 30, 6)
end = dt.datetime(2008, 7, 31, 0)
id = run_ENS(begin=begin, end=end, first_minute=False, depends_on=id)
create_satimages(begin, depends_on=id)
scheduler.py
+ 26
− 14
View file @ 68fc2efe
......@@ -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
scripts/apply_obs_op_dart.py
+ 1
− 1
View file @ 68fc2efe
......@@ -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)
......
scripts/assim_synth_obs.py
+ 8
− 8
View file @ 68fc2efe
......@@ -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
do_error_parametrization = ((obscfg['error_assimilate'] == False) and (obscfg.get('sat_channel') == 6))
if not parametrized:
err_assim = np.zeros(n_obs_3d) + obscfg['error_assimilate']
else: # error parametrization for WV73
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()
......
scripts/create_obsseq.py
+ 9
− 14
View file @ 68fc2efe
......@@ -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')
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)
aso.replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
templates/input.eval.nml
+ 6
− 6
View file @ 68fc2efe
......@@ -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 ,
......
templates/input.nml
+ 6
− 6
View file @ 68fc2efe
......@@ -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 ,
......
templates/input.prioronly.nml
+ 6
− 6
View file @ 68fc2efe
......@@ -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 ,
......
templates/namelist.input
+ 1
− 1
View file @ 68fc2efe
......@@ -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,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment