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

minimal changes

parent d37e5153
Branches
No related tags found
No related merge requests found
import numpy as np
from config import clusters # from . = problem in archivedir from config import clusters # from . = problem in archivedir
cluster = clusters.jet # change cluster configuration here cluster = clusters.jet # change cluster configuration here
...@@ -7,10 +6,10 @@ class ExperimentConfiguration(object): ...@@ -7,10 +6,10 @@ class ExperimentConfiguration(object):
pass pass
exp = ExperimentConfiguration() exp = ExperimentConfiguration()
exp.expname = "exp_v1.22_P3_wbub7_WV73_obs10_loc20" exp.expname = "test_jet" #"exp_v1.22_P3_wbub7_WV62_obs10_loc20_oe1"
exp.model_dx = 2000 exp.model_dx = 2000
exp.n_ens = 40 exp.n_ens = 40
exp.n_nodes = 10 exp.n_nodes = 40
exp.filter_kind = 1 exp.filter_kind = 1
exp.inflation = True exp.inflation = True
...@@ -20,13 +19,13 @@ exp.cov_loc_vert_km_horiz_km = False # (3, 20) ...@@ -20,13 +19,13 @@ exp.cov_loc_vert_km_horiz_km = False # (3, 20)
exp.superob_km = False # False or int (spatial averaging of observations) exp.superob_km = False # False or int (spatial averaging of observations)
exp.use_existing_obsseq = False # False or pathname (use precomputed obs_seq.out files) exp.use_existing_obsseq = False # False or pathname (use precomputed obs_seq.out files)
#exp.use_existing_obsseq = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.21_P2_rr_REFL_obs10_loc20_oe1/obs_seq_out/2008-07-30_%H:%M_obs_seq.out' #exp.use_existing_obsseq = '/jetfs/home/lkugler/data/sim_archive/exp_v1.21_P3_wbub7_VIS_obs10_loc20/obs_seq_out/2008-07-30_%H:%M_obs_seq.out'
#exp.use_existing_obsseq = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.21_P3_wbub7_REFL2D_obs10_loc20_oe5/obs_seq_out/2008-07-30_%H:%M_obs_seq.out' #exp.use_existing_obsseq = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.21_P3_wbub7_REFL2D_obs10_loc20_oe5/obs_seq_out/2008-07-30_%H:%M_obs_seq.out'
#exp.use_existing_obsseq = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.21_P2_rr_VIS_obs20_loc4/obs_seq_out/2008-07-30_%H:%M_obs_seq.out' #exp.use_existing_obsseq = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.21_P2_rr_VIS_obs20_loc4/obs_seq_out/2008-07-30_%H:%M_obs_seq.out'
#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P5+su_nat2/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' #exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P5+su_nat2/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' exp.nature_wrfout = '/jetfs/home/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_Pwbub5_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' #exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_Pwbub5_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.18_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' #exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.18_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P4_nat/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S' #exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P4_nat/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
...@@ -35,7 +34,7 @@ exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_n ...@@ -35,7 +34,7 @@ exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_n
#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.nat.<iens>.wrfprof' #exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.nat.<iens>.wrfprof'
#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.fc.<iens>.wrfprof' #exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.fc.<iens>.wrfprof'
#exp.input_profile = '/home/fs71386/lkugler/data/initial_profiles/wrf/ens/large_mean_error/raso.nat.<iens>.wrfprof' #exp.input_profile = '/home/fs71386/lkugler/data/initial_profiles/wrf/ens/large_mean_error/raso.nat.<iens>.wrfprof'
exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-03-31/raso.fc.<iens>.wrfprof' exp.input_profile = '/jetfs/home/lkugler/data/initial_profiles/wrf/ens/2022-03-31/raso.fc.<iens>.wrfprof'
#exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-03-31/raso.nat.<iens>.wrfprof' #exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-03-31/raso.nat.<iens>.wrfprof'
#exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-05-18/raso.fc.<iens>.wrfprof' #exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-05-18/raso.fc.<iens>.wrfprof'
...@@ -48,18 +47,18 @@ n_obs = 961 # 22500: 2km, 5776: 4km, 121: 30km, 256:16x16 (20km); 961: 10km res ...@@ -48,18 +47,18 @@ n_obs = 961 # 22500: 2km, 5776: 4km, 121: 30km, 256:16x16 (20km); 961: 10km res
vis = dict(plotname='VIS 0.6µm', plotunits='[1]', vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs, kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs,
error_generate=0.03, error_assimilate=0.06, error_generate=0.03, error_assimilate=0.03,
cov_loc_radius_km=20) cov_loc_radius_km=20)
wv62 = dict(plotname='Brightness temperature WV 6.2µm', plotunits='[K]', wv62 = dict(plotname='Brightness temperature WV 6.2µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=5, n_obs=n_obs, kind='MSG_4_SEVIRI_TB', sat_channel=5, n_obs=n_obs,
error_generate=1., error_assimilate=False, error_generate=1., error_assimilate=1.,
cov_loc_radius_km=20) cov_loc_radius_km=20)
wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]', wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs, kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs,
error_generate=1., error_assimilate=False, error_generate=1., error_assimilate=1.,
cov_loc_radius_km=20) cov_loc_radius_km=10)
ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]', ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs, kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs,
...@@ -69,19 +68,19 @@ ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]', ...@@ -69,19 +68,19 @@ ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
radar = dict(plotname='Radar reflectivity', plotunits='[dBz]', radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
kind='RADAR_REFLECTIVITY', n_obs=n_obs, kind='RADAR_REFLECTIVITY', n_obs=n_obs,
error_generate=2.5, error_assimilate=5, error_generate=2.5, error_assimilate=5,
heights=[2000,], #np.arange(2000, 14001, 2000), heights=range(2000, 14001, 2000),
cov_loc_radius_km=20) cov_loc_radius_km=20)
t = dict(plotname='Temperature', plotunits='[K]', t = dict(plotname='Temperature', plotunits='[K]',
kind='RADIOSONDE_TEMPERATURE', n_obs=n_obs, kind='RADIOSONDE_TEMPERATURE', n_obs=n_obs,
error_generate=0.2, error_assimilate=0.4, error_generate=0.2, error_assimilate=0.4,
heights=[5000,], #np.arange(1000, 17001, 2000), heights=[5000,], #range(1000, 17001, 2000),
cov_loc_radius_km=4) cov_loc_radius_km=4)
q = dict(plotname='Specific humidity', plotunits='[kg/kg]', q = dict(plotname='Specific humidity', plotunits='[kg/kg]',
kind='RADIOSONDE_SPECIFIC_HUMIDITY', n_obs=n_obs, kind='RADIOSONDE_SPECIFIC_HUMIDITY', n_obs=n_obs,
error_generate=0., error_assimilate=5*1e-5, error_generate=0., error_assimilate=5*1e-5,
heights=[1000], #np.arange(1000, 17001, 2000), heights=[1000], #range(1000, 17001, 2000),
cov_loc_radius_km=0.1) cov_loc_radius_km=0.1)
t2m = dict(plotname='SYNOP Temperature', plotunits='[K]', t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
...@@ -95,7 +94,7 @@ psfc = dict(plotname='SYNOP Pressure', plotunits='[Pa]', ...@@ -95,7 +94,7 @@ psfc = dict(plotname='SYNOP Pressure', plotunits='[Pa]',
cov_loc_radius_km=32) cov_loc_radius_km=32)
exp.observations = [wv73] exp.observations = [wv62]
exp.update_vars = ['U', 'V', 'W', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC'] exp.update_vars = ['U', 'V', 'W', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC']
#exp.update_vars = ['U', 'V', 'W', 'T', 'PH', 'MU', 'QVAPOR', 'PSFC'] #exp.update_vars = ['U', 'V', 'W', 'T', 'PH', 'MU', 'QVAPOR', 'PSFC']
......
...@@ -94,8 +94,8 @@ jet = ClusterConfig() ...@@ -94,8 +94,8 @@ jet = ClusterConfig()
jet.name = 'jet' jet.name = 'jet'
# binaries # binaries
jet.python = '/jetfs/home/lkugler/miniconda3/bin/python' jet.python = '/jetfs/home/lkugler/miniconda3/envs/DART/bin/python'
jet.python_verif = '/jetfs/home/lkugler/miniconda3/envs/ensdask/bin/python' jet.python_verif = '/jetfs/home/lkugler/miniconda3/envs/enstools/bin/python'
jet.ncks = '/jetfs/spack/opt/spack/linux-rhel8-skylake_avx512/intel-20.0.2/nco-4.9.3-dhlqiyog7howjmaleyfhm6lkt7ra37xf/bin/ncks' jet.ncks = '/jetfs/spack/opt/spack/linux-rhel8-skylake_avx512/intel-20.0.2/nco-4.9.3-dhlqiyog7howjmaleyfhm6lkt7ra37xf/bin/ncks'
jet.ideal = '/jetfs/home/lkugler/bin/ideal-v4.3_v1.22.exe' jet.ideal = '/jetfs/home/lkugler/bin/ideal-v4.3_v1.22.exe'
jet.wrfexe = '/jetfs/home/lkugler/bin/wrf-v4.3_v1.22.exe' jet.wrfexe = '/jetfs/home/lkugler/bin/wrf-v4.3_v1.22.exe'
...@@ -108,14 +108,15 @@ jet.archive_base = '/jetfs/home/lkugler/data/sim_archive/' ...@@ -108,14 +108,15 @@ jet.archive_base = '/jetfs/home/lkugler/data/sim_archive/'
# paths used as input # paths used as input
jet.srcdir = '/jetfs/home/lkugler/data/compile/WRF-4.3/run' jet.srcdir = '/jetfs/home/lkugler/data/compile/WRF-4.3/run'
jet.dart_srcdir = '/jetfs/home/lkugler/data/compile/DART-9.11.9/DART/models/wrf/work' jet.dart_srcdir = '/jetfs/home/lkugler/data/compile/DART/DART-10.5.3/models/wrf/work'
jet.rttov_srcdir = '/jetfs/home/lkugler/data/compile/RTTOV13/rtcoef_rttov13/' jet.rttov_srcdir = '/jetfs/home/lkugler/data/compile/RTTOV13/rtcoef_rttov13/'
jet.scriptsdir = '/jetfs/home/lkugler/DART-WRF/dartwrf/' jet.scriptsdir = '/jetfs/home/lkugler/DART-WRF/dartwrf/'
jet.geo_em = '/jetfs/home/lkugler/data/geo_em.d01.nc'
# templates/run scripts # templates/run scripts
jet.namelist = jet.scriptsdir+'/../templates/namelist.input' jet.namelist = jet.scriptsdir+'/../templates/namelist.input'
jet.run_WRF = '/jetfs/home/lkugler/DART-WRF/dartwrf/run_ens.jet.sh' jet.run_WRF = '/jetfs/home/lkugler/DART-WRF/dartwrf/run_ens.jet.sh'
jet.slurm_cfg = {"account": "", "partition": "compute", jet.slurm_cfg = {"account": "lkugler", "partition": "compute",
"nodes": "1", "ntasks": "1", "ntasks-per-node": "40", "ntasks-per-core": "1", "ntasks": "1", "ntasks-per-core": "1", "mem": "50G",
"mail-type": "FAIL", "mail-user": "lukas.kugler@univie.ac.at"} "mail-type": "FAIL", "mail-user": "lukas.kugler@univie.ac.at"}
...@@ -76,9 +76,7 @@ def link_nature_to_dart_truth(time): ...@@ -76,9 +76,7 @@ def link_nature_to_dart_truth(time):
def prepare_nature_dart(time): def prepare_nature_dart(time):
print("linking nature to DART & georeferencing") print("linking nature to DART & georeferencing")
link_nature_to_dart_truth(time) link_nature_to_dart_truth(time)
wrfout_add_geo.run( wrfout_add_geo.run(cluster.geo_em, cluster.dartrundir + "/wrfout_d01")
cluster.dartrundir + "/../geo_em.d01.nc", cluster.dartrundir + "/wrfout_d01"
)
def prepare_prior_ensemble(assim_time, prior_init_time, prior_valid_time, prior_path_exp): def prepare_prior_ensemble(assim_time, prior_init_time, prior_valid_time, prior_path_exp):
...@@ -114,7 +112,7 @@ def prepare_prior_ensemble(assim_time, prior_init_time, prior_valid_time, prior_ ...@@ -114,7 +112,7 @@ def prepare_prior_ensemble(assim_time, prior_init_time, prior_valid_time, prior_
cluster.dartrundir+"/wrfout_d01 "+ wrfout_dart) cluster.dartrundir+"/wrfout_d01 "+ wrfout_dart)
# this seems to be necessary (else wrong level selection) # this seems to be necessary (else wrong level selection)
wrfout_add_geo.run(cluster.dartrundir + "/../geo_em.d01.nc", wrfout_dart) wrfout_add_geo.run(cluster.geo_em, wrfout_dart)
write_list_of_inputfiles_prior() write_list_of_inputfiles_prior()
write_list_of_outputfiles() write_list_of_outputfiles()
...@@ -168,7 +166,7 @@ def run_perfect_model_obs(nproc=12, verbose=True): ...@@ -168,7 +166,7 @@ def run_perfect_model_obs(nproc=12, verbose=True):
"obs_seq.out does not exist in " + cluster.dartrundir, "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 filter(nproc=48): def filter(nproc=12):
print("time now", dt.datetime.now()) print("time now", dt.datetime.now())
print("running filter") print("running filter")
os.chdir(cluster.dartrundir) os.chdir(cluster.dartrundir)
...@@ -223,15 +221,15 @@ def get_parametrized_error(obscfg, osf_prior): ...@@ -223,15 +221,15 @@ def get_parametrized_error(obscfg, osf_prior):
"""Calculate the parametrized error for an ObsConfig (one obs type) """Calculate the parametrized error for an ObsConfig (one obs type)
Args Args
obscfg () obscfg (object): configuration of observations
df_obs (obsseq.ObsRecord): contains observations from obs_seq.out osf_prior (obsseq.ObsRecord): contains truth and prior values from obs_seq.final
(output of ./final in evaluate-mode (no posterior))
Returns Returns
np.array observation error std-dev for assimilation np.array observation error std-dev for assimilation
""" """
df_obs = osf_prior Hx_prior = osf_prior.get_prior_Hx().T
Hx_prior = df_obs.get_prior_Hx().T Hx_truth = osf_prior.get_truth_Hx()
Hx_truth = df_obs.get_truth_Hx()
# compute the obs error for assimilation on the averaged grid # compute the obs error for assimilation on the averaged grid
# since the assimilation is done on the averaged grid # since the assimilation is done on the averaged grid
...@@ -246,6 +244,17 @@ def get_parametrized_error(obscfg, osf_prior): ...@@ -246,6 +244,17 @@ def get_parametrized_error(obscfg, osf_prior):
def set_obserr_assimilate_in_obsseqout(oso, osf_prior, outfile="./obs_seq.out"): def set_obserr_assimilate_in_obsseqout(oso, osf_prior, outfile="./obs_seq.out"):
""""Overwrite existing variance values in obs_seq.out files
Args:
oso (ObsSeq) : python representation of obs_seq.out file,
will be modified and written to file
osf_prior (ObsSeq): python representation of obs_seq.final (output of filter in evaluate-mode without posterior)
contains prior values; used for parameterized errors
Returns:
None (writes to file)
"""
for obscfg in exp.observations: for obscfg in exp.observations:
kind_str = obscfg['kind'] kind_str = obscfg['kind']
kind = osq.obs_kind_nrs[kind_str] kind = osq.obs_kind_nrs[kind_str]
...@@ -358,10 +367,11 @@ def generate_obsseq_out(time): ...@@ -358,10 +367,11 @@ def generate_obsseq_out(time):
def ensure_physical_vis(oso): # set reflectance < surface albedo to surface albedo def ensure_physical_vis(oso): # set reflectance < surface albedo to surface albedo
print(" 2.2) removing obs below surface albedo ") print(" 2.2) removing obs below surface albedo ")
if_vis_obs = oso.df['kind'].values == 262 clearsky_albedo = 0.2928
if_obs_below_surface_albedo = oso.df['observations'].values < 0.2928
oso.df.loc[if_vis_obs & if_obs_below_surface_albedo, ('observations')] = 0.2928 if_vis_obs = oso.df['kind'].values == 262
if_obs_below_surface_albedo = oso.df['observations'].values < clearsky_albedo
oso.df.loc[if_vis_obs & if_obs_below_surface_albedo, ('observations')] = clearsky_albedo
oso.to_dart(f=cluster.dartrundir + "/obs_seq.out") oso.to_dart(f=cluster.dartrundir + "/obs_seq.out")
return oso return oso
...@@ -384,7 +394,7 @@ def generate_obsseq_out(time): ...@@ -384,7 +394,7 @@ def generate_obsseq_out(time):
dir_obsseq=cluster.archivedir + "/obs_seq_out/" dir_obsseq=cluster.archivedir + "/obs_seq_out/"
os.makedirs(dir_obsseq, exist_ok=True) os.makedirs(dir_obsseq, exist_ok=True)
osq.create_obsseqin_alltypes(time, exp.observations) osq.create_obs_seq_in(time, exp.observations)
run_perfect_model_obs() # generate observation, draws from gaussian run_perfect_model_obs() # generate observation, draws from gaussian
print(" 2.1) obs preprocessing") print(" 2.1) obs preprocessing")
...@@ -411,15 +421,18 @@ def get_obsseq_out(time): ...@@ -411,15 +421,18 @@ def get_obsseq_out(time):
print(f_obsseq, 'copied to', cluster.dartrundir+'/obs_seq.out') print(f_obsseq, 'copied to', cluster.dartrundir+'/obs_seq.out')
oso = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.out") oso = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.out")
else: else:
# decision to NOT use existing obs_seq.out file
# did we already generate observations? # did we already generate observations?
f_oso_thisexp = cluster.archivedir+'/obs_seq_out/'+time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out") # f_oso_thisexp = cluster.archivedir+'/obs_seq_out/'+time.strftime("/%Y-%m-%d_%H:%M_obs_seq.out")
if os.path.isfile(f_oso_thisexp): # if os.path.isfile(f_oso_thisexp):
# oso exists # # oso exists
copy(f_oso_thisexp, cluster.dartrundir+'/obs_seq.out') # copy(f_oso_thisexp, cluster.dartrundir+'/obs_seq.out')
print('copied existing obsseqout from', f_oso_thisexp) # print('copied existing obsseqout from', f_oso_thisexp)
oso = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.out") # oso = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.out")
else: # else:
# we need to generate observations
# generate observations with new observation noise
oso = generate_obsseq_out(time) oso = generate_obsseq_out(time)
return oso return oso
...@@ -453,7 +466,7 @@ if __name__ == "__main__": ...@@ -453,7 +466,7 @@ if __name__ == "__main__":
options = [] options = []
if len(sys.argv) >4: if len(sys.argv) >4:
options = sys.argv[5:] options = sys.argv[5:]
nproc = 6 if 'headnode' in options else 48 nproc = 6 if 'headnode' in options else 12
archive_time = cluster.archivedir + time.strftime("/%Y-%m-%d_%H:%M/") archive_time = cluster.archivedir + time.strftime("/%Y-%m-%d_%H:%M/")
os.makedirs(cluster.dartrundir, exist_ok=True) # create directory to run DART in os.makedirs(cluster.dartrundir, exist_ok=True) # create directory to run DART in
......
...@@ -58,8 +58,8 @@ def OE_model_harnisch_WV62(ci): ...@@ -58,8 +58,8 @@ def OE_model_harnisch_WV62(ci):
# Kelvin, fit of Fig 7a, Harnisch 2016 # Kelvin, fit of Fig 7a, Harnisch 2016
x_ci = [0, 2.5, 4.5, 5.5, 7.5] # average cloud impact [K] x_ci = [0, 2.5, 4.5, 5.5, 7.5] # average cloud impact [K]
y_oe = [1.2, 3, 5, 6, 6.5] # adjusted observation error [K] y_oe = [1.2, 3, 5, 6, 6.5] # adjusted observation error [K]
oe_73_linear = interp1d(x_ci, y_oe, assume_sorted=True) oe_linear = interp1d(x_ci, y_oe, assume_sorted=True)
return oe_73_linear(ci) return oe_linear(ci)
else: # assign highest observation error else: # assign highest observation error
return 6.5 return 6.5
...@@ -69,8 +69,10 @@ def OE_model_harnisch_WV73(ci): ...@@ -69,8 +69,10 @@ def OE_model_harnisch_WV73(ci):
# Kelvin, fit of Fig 7b, Harnisch 2016 # Kelvin, fit of Fig 7b, Harnisch 2016
x_ci = [0, 5, 10.5, 13, 16] # average cloud impact [K] x_ci = [0, 5, 10.5, 13, 16] # average cloud impact [K]
y_oe = [1, 4.5, 10, 12, 13] # adjusted observation error [K] y_oe = [1, 4.5, 10, 12, 13] # adjusted observation error [K]
oe_73_linear = interp1d(x_ci, y_oe, assume_sorted=True)
return oe_73_linear(ci) #y_oe = [1.2, 3, 5, 6, 6.5] # OE for WV62 !!!!
oe_linear = interp1d(x_ci, y_oe, assume_sorted=True)
return oe_linear(ci)
else: # assign highest observation error else: # assign highest observation error
return 13.0 return 13.0
...@@ -143,3 +143,96 @@ def read_prior_obs(f_obs_prior): ...@@ -143,3 +143,96 @@ def read_prior_obs(f_obs_prior):
OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens))) OBSs.append(dict(observed=observed, truth=truth, prior_ens=np.array(prior_ens)))
return OBSs return OBSs
def create_obsseq_in_separate_obs(time_dt, obscfg, obs_errors=False,
archive_obs_coords=False):
"""Create obs_seq.in of one obstype
Args:
time_dt (dt.datetime): time of observation
obscfg (dict)
obs_errors (int, np.array) : values of observation errors (standard deviations)
e.g. 0 = use zero error
archive_obs_coords (str, False): path to folder
channel_id (int): SEVIRI channel number
see https://nwp-saf.eumetsat.int/downloads/rtcoef_rttov12/ir_srf/rtcoef_msg_4_seviri_srf.html
coords (list of 2-tuples with (lat,lon))
obserr_std (np.array): shape (n_obs,), one value for each observation,
gaussian error with this std-dev is added to the truth at observation time
archive_obs_coords (bool, str): False or str (filepath where `obs_seq.in` will be saved)
"""
n_obs = obscfg['n_obs']
coords = calc_obs_locations(n_obs,
coords_from_domaincenter=False,
distance_between_obs_km=obscfg.get('distance_between_obs_km', False),
fpath_coords=archive_obs_coords)
kind = obscfg['kind']
sat_channel = obscfg.get('sat_channel', False)
# determine vertical coordinates
if not sat_channel:
if 'SYNOP' in kind:
vert_coord_sys = "-1" # meters AGL
vert_coords = [2, ]
else:
vert_coord_sys = "3" # meters AMSL
vert_coords = obscfg['heights']
else:
vert_coord_sys = "-2" # undefined height
vert_coords = ["-888888.0000000000", ]
coords = append_hgt_to_coords(coords, vert_coords)
n_obs_3d = len(coords)
# define obs error
obserr_std = np.zeros(n_obs_3d)
if obs_errors:
obserr_std += obs_errors
# other stuff for obsseq.in
obs_kind_nr = obs_kind_nrs[kind]
line_obstypedef = ' '+obs_kind_nr+' '+kind
time_dt = add_timezone_UTC(time_dt)
dart_date_day, secs_thatday = get_dart_date(time_dt)
print('secs, days:', secs_thatday, dart_date_day)
if sat_channel:
sun_az = str(get_azimuth(lat0, lon0, time_dt))
sun_zen = str(90. - get_altitude(lat0, lon0, time_dt))
print('sunzen', sun_zen, 'sunazi', sun_az)
appendix = """visir
"""+sat_az+""" """+sat_zen+""" """+sun_az+"""
"""+sun_zen+"""
12 4 21 """+str(sat_channel)+"""
-888888.000000000
1"""
else:
appendix = ''
txt = preamble(n_obs_3d, line_obstypedef)
for i_obs in range(n_obs_3d):
last = False
if i_obs == int(n_obs_3d)-1:
last = True # last_observation
txt += write_section(dict(i=i_obs+1,
kind_nr=obs_kind_nr,
dart_date_day=dart_date_day,
secs_thatday=secs_thatday,
lon=coords[i_obs][1],
lat=coords[i_obs][0],
vert_coord=coords[i_obs][2],
vert_coord_sys=vert_coord_sys,
obserr_var=obserr_std[i_obs]**2,
appendix=appendix),
last=last)
write_file(txt, output_path=cluster.dartrundir+'/obs_seq.in')
...@@ -5,8 +5,8 @@ import subprocess ...@@ -5,8 +5,8 @@ import subprocess
def shell(args): def shell(args):
print(args) print(args)
subprocess.run(args.split(' ')) #, shell=True) #, stderr=subprocess.STDOUT) #subprocess.run(args.split(' ')) #, shell=True) #, stderr=subprocess.STDOUT)
#os.system(args) os.system(args)
def print(*args): def print(*args):
__builtin__.print(*args, flush=True) __builtin__.print(*args, flush=True)
......
name: DART
channels:
- conda-forge
- defaults
dependencies:
- cdo
- docopt
- hdf5
- ipython
- matplotlib
- nco
- netcdf-fortran
- netcdf4
- numpy
- pandas
- pip
- python
- python-dateutil
- pytz
- proj4
- scipy
- setuptools
- wrf-python
- xarray
- yaml
- pysolar
...@@ -75,8 +75,8 @@ do ...@@ -75,8 +75,8 @@ do
mv $rundir/rsl.out.0000 $rundir/rsl.out.input mv $rundir/rsl.out.0000 $rundir/rsl.out.input
done done
""" """
s = my_Slurm("ideal", cfg_update={"ntasks": str(exp.n_ens), "nodes": "1", s = my_Slurm("ideal", cfg_update={"ntasks": str(exp.n_ens),
"time": "10", "mem-per-cpu": "2G"}) "time": "10", "mem": "100G"})
id = s.run(cmd, depends_on=[depends_on]) id = s.run(cmd, depends_on=[depends_on])
return id return id
...@@ -141,13 +141,13 @@ def run_ENS(begin, end, depends_on=None, first_minute=True, ...@@ -141,13 +141,13 @@ def run_ENS(begin, end, depends_on=None, first_minute=True,
if output_restart_interval: if output_restart_interval:
args.append('--restart_interval='+str(int(float(output_restart_interval)))) args.append('--restart_interval='+str(int(float(output_restart_interval))))
s = my_Slurm("preWRF2", cfg_update=dict(time="2")) s = my_Slurm("preWRF", cfg_update=dict(time="2"))
id = s.run(' '.join(args), depends_on=[id]) id = s.run(' '.join(args), depends_on=[id])
time_in_simulation_hours = (end-begin).total_seconds()/3600 time_in_simulation_hours = (end-begin).total_seconds()/3600
runtime_wallclock_mins_expected = int(8+time_in_simulation_hours*9.5) # usually below 9 min/hour runtime_wallclock_mins_expected = int(8+time_in_simulation_hours*9.5) # usually below 9 min/hour
s = my_Slurm("runWRF2", cfg_update={"nodes": "1", "array": "1-"+str(exp.n_nodes), s = my_Slurm("WRF", cfg_update={"array": "1-"+str(exp.n_nodes), "ntasks": "10", "nodes": "1",
"time": str(runtime_wallclock_mins_expected), "mem-per-cpu": "2G"}) "time": str(runtime_wallclock_mins_expected), "mem": "140G"})
cmd = script_to_str(cluster.run_WRF).replace('<exp.expname>', exp.expname cmd = script_to_str(cluster.run_WRF).replace('<exp.expname>', exp.expname
).replace('<cluster.wrf_rundir_base>', cluster.wrf_rundir_base) ).replace('<cluster.wrf_rundir_base>', cluster.wrf_rundir_base)
id = s.run(cmd, depends_on=[id]) id = s.run(cmd, depends_on=[id])
...@@ -166,8 +166,8 @@ def assimilate(assim_time, prior_init_time, prior_valid_time, prior_path_exp, ...@@ -166,8 +166,8 @@ def assimilate(assim_time, prior_init_time, prior_valid_time, prior_path_exp,
if not os.path.exists(prior_path_exp): if not os.path.exists(prior_path_exp):
raise IOError('prior_path_exp does not exist: '+prior_path_exp) raise IOError('prior_path_exp does not exist: '+prior_path_exp)
id = my_Slurm("Assim", cfg_update={"nodes": "1", "ntasks": "96", "time": "60", id = my_Slurm("Assim", cfg_update={"ntasks": "12", "time": "60",
"mem": "300G", "ntasks-per-node": "96", "ntasks-per-core": "2"} "mem": "200G", "ntasks-per-node": "12", "ntasks-per-core": "2"}
).run(cluster.python+' '+cluster.scripts_rundir+'/assim_synth_obs.py ' ).run(cluster.python+' '+cluster.scripts_rundir+'/assim_synth_obs.py '
+assim_time.strftime('%Y-%m-%d_%H:%M ') +assim_time.strftime('%Y-%m-%d_%H:%M ')
+prior_init_time.strftime('%Y-%m-%d_%H:%M ') +prior_init_time.strftime('%Y-%m-%d_%H:%M ')
...@@ -200,8 +200,8 @@ def update_IC_from_DA(assim_time, depends_on=None): ...@@ -200,8 +200,8 @@ def update_IC_from_DA(assim_time, depends_on=None):
def create_satimages(init_time, depends_on=None): def create_satimages(init_time, depends_on=None):
s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "80", "nodes": "1"}) s = my_Slurm("RTTOV", cfg_update={"ntasks": "12", "time": "80", "mem": "200G"})
id = s.run('/home/fs71386/lkugler/RTTOV-WRF/withmodules /home/fs71386/lkugler/RTTOV-WRF/run_init.py '+cluster.archivedir id = s.run(cluster.python_verif+' ~/RTTOV-WRF/run_init.py '+cluster.archivedir
+init_time.strftime('/%Y-%m-%d_%H:%M/'), +init_time.strftime('/%Y-%m-%d_%H:%M/'),
depends_on=[depends_on]) depends_on=[depends_on])
return id return id
...@@ -221,19 +221,19 @@ def gen_obsseq(depends_on=None): ...@@ -221,19 +221,19 @@ def gen_obsseq(depends_on=None):
def verify_sat(depends_on=None): def verify_sat(depends_on=None):
s = my_Slurm("verif-SAT-"+exp.expname, cfg_update={"time": "120", "mail-type": "FAIL,END", "ntasks": "1", s = my_Slurm("verif-SAT-"+exp.expname, cfg_update={"time": "60", "mail-type": "FAIL,END", "ntasks": "20",
"ntasks-per-node": "1", "ntasks-per-core": "1"}) "ntasks-per-node": "20", "ntasks-per-core": "1", "mem": "100G",})
cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node sat verif1d FSS BS' cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node sat verif1d FSS BS'
s.run(cmd, depends_on=[depends_on]) s.run(cmd, depends_on=[depends_on])
def verify_wrf(depends_on=None): def verify_wrf(depends_on=None):
s = my_Slurm("verif-WRF-"+exp.expname, cfg_update={"time": "180", "mail-type": "FAIL,END", "ntasks": "1", s = my_Slurm("verif-WRF-"+exp.expname, cfg_update={"time": "120", "mail-type": "FAIL,END", "ntasks": "20",
"ntasks-per-node": "1", "ntasks-per-core": "1"}) "ntasks-per-node": "20", "ntasks-per-core": "1", "mem": "250G"})
cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node wrf verif1d FSS BS' cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node wrf verif1d FSS BS'
s.run(cmd, depends_on=[depends_on]) s.run(cmd, depends_on=[depends_on])
def verify_fast(depends_on=None): def verify_fast(depends_on=None):
s = my_Slurm("verif-fast-"+exp.expname, cfg_update={"time": "30", "mail-type": "FAIL", "ntasks": "1", s = my_Slurm("verif-fast-"+exp.expname, cfg_update={"time": "10", "mail-type": "FAIL", "ntasks": "1",
"ntasks-per-node": "1", "ntasks-per-core": "1"}) "ntasks-per-node": "1", "ntasks-per-core": "1"})
cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_fast/plot_single_exp.py '+exp.expname cmd = cluster.python_verif+' /jetfs/home/lkugler/osse_analysis/plot_fast/plot_single_exp.py '+exp.expname
s.run(cmd, depends_on=[depends_on]) s.run(cmd, depends_on=[depends_on])
...@@ -316,7 +316,7 @@ if __name__ == "__main__": ...@@ -316,7 +316,7 @@ if __name__ == "__main__":
# update time variables # update time variables
prior_init_time = time - timedelta_btw_assim prior_init_time = time - timedelta_btw_assim
id = gen_obsseq(id) #id = gen_obsseq(id)
verify_sat(id_sat) verify_sat(id_sat)
verify_wrf(id) verify_wrf(id)
verify_fast(id) verify_fast(id)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment