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

scripts into ifmain

parent 6eb26f5b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dartwrf/create_wbubble_wrfinput.py
+ 65
− 64
View file @ f3947c60
...@@ -11,78 +11,79 @@ dx_km = 2 ...@@ -11,78 +11,79 @@ dx_km = 2
cr = 15 # km horizontal relaxation distance cr = 15 # km horizontal relaxation distance
cz = 2000 # meters vertical relaxation distance cz = 2000 # meters vertical relaxation distance
args = sys.argv[1:] if __name__ == "__main__":
args = sys.argv[1:]
if 'plot' in args: if 'plot' in args:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(5,5)) fig, ax = plt.subplots(figsize=(5,5))
perturbations = True if 'perturb' in args else False perturbations = True if 'perturb' in args else False
print('update state in wrfinput wrfout file to DART background file') print('update state in wrfinput wrfout file to DART background file')
print('perturb wbubble = ', perturbations) print('perturb wbubble = ', perturbations)
for iens in range(1, exp.n_ens+1): for iens in range(1, exp.n_ens+1):
print('iens', iens) print('iens', iens)
wrfin = cluster.wrf_rundir(iens)+'/wrfinput_d01' wrfin = cluster.wrf_rundir(iens)+'/wrfinput_d01'
# insert wbubble # insert wbubble
with nc.Dataset(wrfin, 'r+') as ds: with nc.Dataset(wrfin, 'r+') as ds:
t = ds.variables['T'][:] t = ds.variables['T'][:]
nt, nz, ny, nx = t.shape nt, nz, ny, nx = t.shape
z = ds.variables['PHB'][0,:,100,100]/9.81 z = ds.variables['PHB'][0,:,100,100]/9.81
z = np.array(z) z = np.array(z)
z = z - z[0] z = z - z[0]
z = (z[1:]+z[:-1])/2 z = (z[1:]+z[:-1])/2
z = z[:, np.newaxis, np.newaxis] z = z[:, np.newaxis, np.newaxis]
if perturbations: if perturbations:
cx = (70 + 60*np.random.rand())*dx_km cx = (70 + 60*np.random.rand())*dx_km
cy = (70 + 60*np.random.rand())*dx_km cy = (70 + 60*np.random.rand())*dx_km
else: else:
cx = 100*dx_km cx = 100*dx_km
cy = 100*dx_km cy = 100*dx_km
print('wbubble center is', cx, cy) print('wbubble center is', cx, cy)
x = np.arange(0,nx)*dx_km x = np.arange(0,nx)*dx_km
y = np.arange(0,ny)*dx_km y = np.arange(0,ny)*dx_km
dx = x-cx dx = x-cx
dy = y-cy dy = y-cy
xx, yy = np.meshgrid(dx, dy) xx, yy = np.meshgrid(dx, dy)
dr = np.sqrt(xx**2 + yy**2)[np.newaxis, :, :] dr = np.sqrt(xx**2 + yy**2)[np.newaxis, :, :]
amplitude = 3 amplitude = 3
if perturbations: if perturbations:
amplitude += np.random.rand()*2 - 1 # Uniformly random +/- 1 K amplitude += np.random.rand()*2 - 1 # Uniformly random +/- 1 K
increment = amplitude*np.exp(-(dr/cr)**2)*np.exp(-(z/cz)**2) increment = amplitude*np.exp(-(dr/cr)**2)*np.exp(-(z/cz)**2)
# now perturbations are centered in domain # now perturbations are centered in domain
# to shift it in direction NW # to shift it in direction NW
shift = 50 shift = 50
incr2 = increment.copy()*0. incr2 = increment.copy()*0.
incr2[:, :-shift, :-shift] = increment[:, shift:, shift:] # main area incr2[:, :-shift, :-shift] = increment[:, shift:, shift:] # main area
incr2[:, -shift:, :-shift] = increment[:, :shift, shift:] # lower part incr2[:, -shift:, :-shift] = increment[:, :shift, shift:] # lower part
incr2[:, :-shift, -shift:] = increment[:, shift:, :shift] # right part incr2[:, :-shift, -shift:] = increment[:, shift:, :shift] # right part
incr2[:, -shift:, -shift:] = increment[:, :shift, :shift] # bottom right corner incr2[:, -shift:, -shift:] = increment[:, :shift, :shift] # bottom right corner
if 'plot' in args: if 'plot' in args:
pdata = incr2[0,:,:] # select lowest level pdata = incr2[0,:,:] # select lowest level
pdata = pdata - pdata.min() pdata = pdata - pdata.min()
c = next(ax._get_lines.prop_cycler)['color'] c = next(ax._get_lines.prop_cycler)['color']
ax.contour(pdata, levels=[2.,], linewidths=1, colors=[c], zorder=3) ax.contour(pdata, levels=[2.,], linewidths=1, colors=[c], zorder=3)
ax.contour(pdata, levels=[1.,], linewidths=0.5, colors=[c], zorder=3) ax.contour(pdata, levels=[1.,], linewidths=0.5, colors=[c], zorder=3)
ds.variables['T'][0,...] += incr2 ds.variables['T'][0,...] += incr2
ds.variables['THM'][0,...] += incr2 ds.variables['THM'][0,...] += incr2
print(wrfin, 'wbubble inserted.') print(wrfin, 'wbubble inserted.')
if 'plot' in args: if 'plot' in args:
ax.set_aspect(1) ax.set_aspect(1)
ax.plot([], [], color='k', lw=1, label='+2 K perturbation') ax.plot([], [], color='k', lw=1, label='+2 K perturbation')
ax.plot([], [], color='k', lw=0.5, label='+1 K perturbation') ax.plot([], [], color='k', lw=0.5, label='+1 K perturbation')
ax.legend() ax.legend()
fout = '/home/fs71386/lkugler/data/analysis/'+exp.expname+'/shiftedbubbles.png' fout = '/home/fs71386/lkugler/data/analysis/'+exp.expname+'/shiftedbubbles.png'
os.makedirs(os.path.dirname(fout), exist_ok=True) os.makedirs(os.path.dirname(fout), exist_ok=True)
fig.savefig(fout, dpi=200) fig.savefig(fout, dpi=200)
print(fout, 'saved.') print(fout, 'saved.')
dartwrf/link_dart_rttov.py
+ 26
− 24
View file @ f3947c60
...@@ -5,34 +5,36 @@ from utils import symlink, copy_scp_srvx8, copy, sed_inplace ...@@ -5,34 +5,36 @@ from utils import symlink, copy_scp_srvx8, copy, sed_inplace
joinp = os.path.join joinp = os.path.join
# DART executables if __name__ == "__main__":
bins = ['perfect_model_obs', 'filter', 'obs_diag', 'obs_seq_to_netcdf']
for b in bins:
symlink(joinp(cluster.dart_srcdir, b),
joinp(cluster.dartrundir, b))
print(joinp(cluster.dartrundir, b), 'created')
rttov_files = ['rttov13pred54L/rtcoef_msg_4_seviri_o3.dat', # DART executables
#'mfasis_lut/rttov_mfasis_cld_msg_4_seviri_deff.dat', bins = ['perfect_model_obs', 'filter', 'obs_diag', 'obs_seq_to_netcdf']
'mfasis_lut/rttov_mfasis_cld_msg_4_seviri_deff.H5', for b in bins:
'cldaer_visir/sccldcoef_msg_4_seviri.dat'] symlink(joinp(cluster.dart_srcdir, b),
joinp(cluster.dartrundir, b))
print(joinp(cluster.dartrundir, b), 'created')
for f_src in rttov_files: rttov_files = ['rttov13pred54L/rtcoef_msg_4_seviri_o3.dat',
destname = os.path.basename(f_src) #'mfasis_lut/rttov_mfasis_cld_msg_4_seviri_deff.dat',
if 'rtcoef' in f_src: 'mfasis_lut/rttov_mfasis_cld_msg_4_seviri_deff.H5',
destname = 'rtcoef_msg_4_seviri.dat' 'cldaer_visir/sccldcoef_msg_4_seviri.dat']
symlink(cluster.rttov_srcdir + f_src, for f_src in rttov_files:
cluster.dartrundir+'/'+destname) destname = os.path.basename(f_src)
if 'rtcoef' in f_src:
destname = 'rtcoef_msg_4_seviri.dat'
################## symlink(cluster.rttov_srcdir + f_src,
symlink(cluster.dartrundir+'/rttov_mfasis_cld_msg_4_seviri_deff.H5', cluster.dartrundir+'/'+destname)
cluster.dartrundir+'/rttov_mfasis_cld_msg_4_seviri.H5')
symlink(cluster.dart_srcdir+'/../../../observations/forward_operators/rttov_sensor_db.csv', ##################
cluster.dartrundir+'/rttov_sensor_db.csv') symlink(cluster.dartrundir+'/rttov_mfasis_cld_msg_4_seviri_deff.H5',
cluster.dartrundir+'/rttov_mfasis_cld_msg_4_seviri.H5')
symlink(cluster.dart_srcdir+'/../../../assimilation_code/programs/gen_sampling_err_table/work/sampling_error_correction_table.nc', symlink(cluster.dart_srcdir+'/../../../observations/forward_operators/rttov_sensor_db.csv',
cluster.dartrundir+'/sampling_error_correction_table.nc') cluster.dartrundir+'/rttov_sensor_db.csv')
print('prepared DART & RTTOV links in', cluster.dartrundir) symlink(cluster.dart_srcdir+'/../../../assimilation_code/programs/gen_sampling_err_table/work/sampling_error_correction_table.nc',
cluster.dartrundir+'/sampling_error_correction_table.nc')
print('prepared DART & RTTOV links in', cluster.dartrundir)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment