Skip to content
Snippets Groups Projects
Commit 0da787c8 authored by lkugler's avatar lkugler
Browse files

bugfix

parent 66349851
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
scripts/assim_synth_obs.py
+ 55
− 44
View file @ 0da787c8
...@@ -74,6 +74,7 @@ def replace_errors_obsseqout(f, new_errors): ...@@ -74,6 +74,7 @@ def replace_errors_obsseqout(f, new_errors):
new_errors (np.array) : standard deviation, new_errors (np.array) : standard deviation,
shape must match the number of observations shape must match the number of observations
""" """
debug = True
obsseq = open(f, 'r').readlines() obsseq = open(f, 'r').readlines()
# find number of lines between two ' OBS ' lines: # find number of lines between two ' OBS ' lines:
...@@ -96,7 +97,7 @@ def replace_errors_obsseqout(f, new_errors): ...@@ -96,7 +97,7 @@ def replace_errors_obsseqout(f, new_errors):
previous_err_var = obsseq[line_error_obs_i] previous_err_var = obsseq[line_error_obs_i]
new_err_obs_i = new_errors[i_obs]**2 # variance in obs_seq.out new_err_obs_i = new_errors[i_obs]**2 # variance in obs_seq.out
#print('previous err var ', previous_err_var, 'new error', new_err_obs_i) if debug: print('previous err var ', previous_err_var, 'new error', new_err_obs_i)
obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n' obsseq[line_error_obs_i] = ' '+str(new_err_obs_i)+' \n'
i_obs += 1 # next iteration i_obs += 1 # next iteration
...@@ -189,9 +190,10 @@ def obs_operator_ensemble(istage): ...@@ -189,9 +190,10 @@ def obs_operator_ensemble(istage):
# DART may need a wrfinput file as well, which serves as a template for dimension sizes # DART may need a wrfinput file as well, which serves as a template for dimension sizes
symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01') symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
# I dont think this is necessary, we do this already in pre_assim.py
# add geodata, if istage>0, wrfout is DART output (has coords) # add geodata, if istage>0, wrfout is DART output (has coords)
if istage == 0: #if istage == 0:
wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01') # wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
# run perfect_model obs (forward operator) # run perfect_model obs (forward operator)
os.system('mpirun -np 12 ./perfect_model_obs > /dev/null') os.system('mpirun -np 12 ./perfect_model_obs > /dev/null')
...@@ -221,6 +223,7 @@ def link_nature_to_dart_truth(time): ...@@ -221,6 +223,7 @@ def link_nature_to_dart_truth(time):
cluster.dartrundir+'/wrfout_d01') cluster.dartrundir+'/wrfout_d01')
# DART may need a wrfinput file as well, which serves as a template for dimension sizes # DART may need a wrfinput file as well, which serves as a template for dimension sizes
symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01') symlink(cluster.dartrundir+'/wrfout_d01', cluster.dartrundir+'/wrfinput_d01')
print('linked', time.strftime(cluster.nature_wrfout), 'to', cluster.dartrundir+'/wrfout_d01')
def prepare_nature_dart(time): def prepare_nature_dart(time):
...@@ -254,7 +257,7 @@ def run_perfect_model_obs(): ...@@ -254,7 +257,7 @@ def run_perfect_model_obs():
os.system('mpirun -np 12 ./perfect_model_obs > log.perfect_model_obs') os.system('mpirun -np 12 ./perfect_model_obs > log.perfect_model_obs')
if not os.path.exists(cluster.dartrundir+'/obs_seq.out'): if not os.path.exists(cluster.dartrundir+'/obs_seq.out'):
raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir, raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir,
'\n look for '+cluster.dartrundir+'log.perfect_model_obs') '\n look for '+cluster.dartrundir+'/log.perfect_model_obs')
def assimilate(nproc=96): def assimilate(nproc=96):
print('time now', dt.datetime.now()) print('time now', dt.datetime.now())
...@@ -279,6 +282,8 @@ def archive_diagnostics(archive_dir, time): ...@@ -279,6 +282,8 @@ def archive_diagnostics(archive_dir, time):
# warnings.warn(str(e)) # warnings.warn(str(e))
def recycle_output(): def recycle_output():
"""Use output of assimilation (./filter) as input for another assimilation (with ./filter)
Specifically, this copies the state fields from filter_restart_d01.000x to the wrfout files in advance_temp folders"""
update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'QRAIN', 'U10', 'V10', 'T2', 'Q2', 'TSK', 'PSFC', 'CLDFRA'] update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'QRAIN', 'U10', 'V10', 'T2', 'Q2', 'TSK', 'PSFC', 'CLDFRA']
updates = ','.join(update_vars) updates = ','.join(update_vars)
...@@ -344,6 +349,9 @@ if __name__ == "__main__": ...@@ -344,6 +349,9 @@ if __name__ == "__main__":
Assumptions: Assumptions:
- x_ensemble is already linked for DART to advance_temp<iens>/wrfout_d01 - x_ensemble is already linked for DART to advance_temp<iens>/wrfout_d01
Example call:
python assim.py 2008-08-07_12:00
""" """
time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M') time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
...@@ -360,20 +368,21 @@ if __name__ == "__main__": ...@@ -360,20 +368,21 @@ if __name__ == "__main__":
################################################ ################################################
print(' 1) get the assimilation errors in a single vector ') print(' 1) get the assimilation errors in a single vector ')
error_assimilate = [] error_generate = []
for i, obscfg in enumerate(exp.observations): for i, obscfg in enumerate(exp.observations):
n_obs = obscfg['n_obs'] n_obs = obscfg['n_obs']
n_obs_z = obscfg.get('heights', 1) n_obs_z = len(obscfg.get('heights', [1,]))
n_obs_3d = n_obs * n_obs_z n_obs_3d = n_obs * n_obs_z
parametrized = obscfg.get('sat_channel') == 6 parametrized = obscfg.get('sat_channel') == 6
if not parametrized: if not parametrized:
err_this_type = np.zeros(n_obs_3d) + obscfg['error_generate'] err_this_type = np.zeros(n_obs_3d) + obscfg['error_generate']
else: # error parametrization for WV73 else: # error parametrization for WV73
# get observations for sat 6 # get observations for sat 6
osq.create_obsseqin_alltypes(time, [obscfg,], obs_errors='error_generate') osq.create_obsseqin_alltypes(time, [obscfg,], np.zeros(n_obs_3d))
run_perfect_model_obs() run_perfect_model_obs()
# depends on obs_seq.out produced before by run_perfect_model_obs() # depends on obs_seq.out produced before by run_perfect_model_obs()
...@@ -382,11 +391,11 @@ if __name__ == "__main__": ...@@ -382,11 +391,11 @@ if __name__ == "__main__":
Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
err_this_type = calc_obserr_WV73(Hx_nat, Hx_prior) err_this_type = calc_obserr_WV73(Hx_nat, Hx_prior)
error_assimilate.extend(err_this_type) error_generate.extend(err_this_type) # the obs-error we assume for generating observations
################################################ ################################################
print(' 2) generate observations ') print(' 2) generate observations ')
osq.create_obsseqin_alltypes(time, exp.observations, obs_errors='error_generate', osq.create_obsseqin_alltypes(time, exp.observations, obs_errors=error_generate,
archive_obs_coords=archive_stage+'/obs_coords.pkl') archive_obs_coords=archive_stage+'/obs_coords.pkl')
first_obstype = exp.observations[0] first_obstype = exp.observations[0]
...@@ -397,58 +406,60 @@ if __name__ == "__main__": ...@@ -397,58 +406,60 @@ if __name__ == "__main__":
################################################ ################################################
print(' 3) assimilate with adapted observation-errors ') print(' 3) assimilate with adapted observation-errors ')
error_assimilate = np.zeros(n_obs_3d) + obscfg['error_assimilate'] # the obs-error we assume for assimilation
replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate) replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
t = time.time() t = time_module.time()
assimilate() assimilate()
print('filter took', time_module.time()-t, 'seconds')
dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage) dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage)
archive_diagnostics(dir_obsseq, time) archive_diagnostics(dir_obsseq, time)
archive_output(archive_stage) archive_output(archive_stage)
sys.exit() # below is the code for separate assimilation of different obs types #sys.exit() # below is the code for separate assimilation of different obs types
for istage, obscfg in enumerate(exp.observations): #for istage, obscfg in enumerate(exp.observations):
print('running observation stage', istage, obscfg) # print('running observation stage', istage, obscfg)
archive_stage = archive_time + '/assim_stage'+str(istage)+'/' # archive_stage = archive_time + '/assim_stage'+str(istage)+'/'
n_obs = obscfg['n_obs'] # n_obs = obscfg['n_obs']
n_obs_3d = n_obs * len(obscfg['heights']) # n_obs_3d = n_obs * len(obscfg['heights'])
sat_channel = obscfg.get('sat_channel', False) # sat_channel = obscfg.get('sat_channel', False)
# debug option: overwrite time in prior files # # debug option: overwrite time in prior files
# for iens in range(1,41): # # for iens in range(1,41):
# os.system('ncks -A -v Times '+cluster.dartrundir+'/wrfout_d01 '+cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01') # # os.system('ncks -A -v Times '+cluster.dartrundir+'/wrfout_d01 '+cluster.dartrundir+'/advance_temp'+str(iens)+'/wrfout_d01')
if error_assimilate == False: # if error_assimilate == False:
# use error parametrization for assimilation-obs.errors # # use error parametrization for assimilation-obs.errors
if sat_channel != 6: # if sat_channel != 6:
raise NotImplementedError('sat channel '+str(sat_channel)) # raise NotImplementedError('sat channel '+str(sat_channel))
# depends on obs_seq.out produced before by run_perfect_model_obs() # # depends on obs_seq.out produced before by run_perfect_model_obs()
Hx_nat, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out') # Hx_nat, _ = read_truth_obs_obsseq(cluster.dartrundir+'/obs_seq.out')
Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory # Hx_prior = obs_operator_ensemble(istage) # files are already linked to DART directory
error_assimilate = calc_obserr_WV73(Hx_nat, Hx_prior) # error_assimilate = calc_obserr_WV73(Hx_nat, Hx_prior)
else: # else:
# overwrite error values in obs_seq.out # # overwrite error values in obs_seq.out
error_assimilate = np.zeros(n_obs_3d) + error_assimilate # ensure shape # error_assimilate = np.zeros(n_obs_3d) + error_assimilate # ensure shape
replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate) # replace_errors_obsseqout(cluster.dartrundir+'/obs_seq.out', error_assimilate)
assimilate() # assimilate()
dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage) # dir_obsseq = cluster.archivedir+'/obs_seq_final/assim_stage'+str(istage)
archive_diagnostics(dir_obsseq, time) # archive_diagnostics(dir_obsseq, time)
if istage < n_stages-1: # if istage < n_stages-1:
# recirculation: filter output -> input # # recirculation: filter output -> input
archive_output(archive_stage) # archive_output(archive_stage)
recycle_output() # recycle_output()
elif istage == n_stages-1: # elif istage == n_stages-1:
# last assimilation, continue integration now # # last assimilation, continue integration now
archive_output(archive_stage) # archive_output(archive_stage)
else: # else:
RuntimeError('this should never occur?!') # RuntimeError('this should never occur?!')
scripts/create_obsseq.py
+ 10
− 6
View file @ 0da787c8
...@@ -143,7 +143,7 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True, ...@@ -143,7 +143,7 @@ def calc_obs_locations(n_obs, coords_from_domaincenter=True,
omit_covloc_radius_on_boundary = True omit_covloc_radius_on_boundary = True
if omit_covloc_radius_on_boundary: # in order to avoid an increment step on the boundary if omit_covloc_radius_on_boundary: # in order to avoid an increment step on the boundary
skip_km = cov_loc_radius_km*1.5 skip_km = 50 # cov_loc_radius_km*1.5
skip_gridpoints = int(skip_km/model_dx_km) # skip this many gridpoints on each side skip_gridpoints = int(skip_km/model_dx_km) # skip this many gridpoints on each side
gridpoints_left = nx - 2*skip_gridpoints gridpoints_left = nx - 2*skip_gridpoints
...@@ -376,13 +376,14 @@ def create_obsseq_in_separate_obs(time_dt, obscfg, obs_errors=False, ...@@ -376,13 +376,14 @@ def create_obsseq_in_separate_obs(time_dt, obscfg, obs_errors=False,
write_file(txt, output_path=cluster.dartrundir+'/obs_seq.in') write_file(txt, output_path=cluster.dartrundir+'/obs_seq.in')
def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors='error_assimilate', archive_obs_coords=False): def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors, archive_obs_coords=False):
"""Create obs_seq.in with multiple obs types in one file """Create obs_seq.in with multiple obs types in one file
Args: Args:
time_dt (dt.datetime): time of observation time_dt (dt.datetime): time of observation
list_obscfg (list of dict) list_obscfg (list of dict)
obs_errors (False or str): Key to obsdict, which field contains the error values obs_errors (list of float, False): contains observation errors, one for each observation
if False: use zero error
archive_obs_coords (bool, str): False or str (filepath where `obs_seq.in` will be saved) archive_obs_coords (bool, str): False or str (filepath where `obs_seq.in` will be saved)
""" """
print('creating obs_seq.in:') print('creating obs_seq.in:')
...@@ -411,9 +412,12 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors='error_assimilate' ...@@ -411,9 +412,12 @@ def create_obsseqin_alltypes(time_dt, list_obscfg, obs_errors='error_assimilate'
n_obs_3d_thistype = len(coords) n_obs_3d_thistype = len(coords)
# define obs error # define obs error
obserr_std = np.zeros(n_obs_3d_thistype) if obs_errors == False:
if obs_errors: obs_errors = np.zeros(n_obs_3d_thistype)
obserr_std += obscfg[obs_errors] 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]
sat_info = write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt) sat_info = write_sat_angle_appendix(sat_channel, lat0, lon0, time_dt)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment