new obstypes, cloud dep. error, other updates

config/cfg.py

@@ -7,16 +7,17 @@ class ExperimentConfiguration(object):
         pass
 
 exp = ExperimentConfiguration()
-exp.expname = "exp_v1.11_LMU_filter_domainobs"
+exp.expname = "exp_v1.12_LMU_WV73_cde"
 exp.model_dx = 2000
 exp.timestep = 10
-exp.n_ens = 40
-exp.n_nodes = 10
-exp.n_obs = 100
-
-exp.sat_channels = [1,]
-exp.error_variance = 0.0009
+exp.n_ens = 20
+exp.n_nodes = 5
+exp.n_obs = 100  # radar: n_obs for each observation height level
 
+exp.sat_channels = [6,]
+exp.sat_err = 0.03
+exp.radar_err = 5.
+exp.distance_between_obs_meters = 10000
 
 # directory paths depend on the name of the experiment
 cluster.expname = exp.expname

config/clusters.py

@@ -17,7 +17,7 @@ class ClusterConfig(object):
 
 vsc = ClusterConfig()
 vsc.name = 'vsc'
-vsc.python = '/home/fs71386/lkugler/miniconda3/bin/python'
+vsc.python = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/python'
 vsc.ncks = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/ncks'
 vsc.userdir = '/home/fs71386/lkugler'
 vsc.srcdir = '/home/fs71386/lkugler/compile/WRF/WRF-4.2.1/run'

scripts/archive_assim.py

@@ -34,15 +34,17 @@ try:
     for iens in range(1, exp.n_ens+1):
         savedir = cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/')+str(iens)
         mkdir(savedir)
-        
+
         copy(cluster.dartrundir+'/input.nml',
              cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M/input.nml'))
 
         filter_in = cluster.dartrundir+'/preassim_member_'+str(iens).zfill(4)+'.nc'
         filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
 
-        copy(filter_in, savedir+time.strftime('/%Y-%m-%d_%H:%M_prior'))
-        copy(filter_out, savedir+time.strftime('/%Y-%m-%d_%H:%M_posterior'))
+        # files not necessary as wrfouts exist
+        #copy(filter_in, savedir+time.strftime('/%Y-%m-%d_%H:%M_prior'))
+        #copy(filter_out, savedir+time.strftime('/%Y-%m-%d_%H:%M_posterior'))
+
 
     # copy mean and sd to archive
     print('copy preassim, postassim mean and sd')

scripts/gen_synth_obs.py

+import os, sys, shutil
+import datetime as dt
+import numpy as np
+from scipy.interpolate import interp1d
+from config.cfg import exp, cluster
+from utils import symlink, copy, sed_inplace, append_file
+import create_obsseq
+
+# fit of Fig 7, Harnisch 2016
+x_ci = [0,   5, 10.5, 13, 16]
+y_oe = [1, 4.5,   10, 12, 13]  # Kelvin
+oe_73_linear = interp1d(x_ci, y_oe, assume_sorted=True)
+
+def oe_73(ci):
+    if ci < 13:
+        return oe_73_linear(ci)
+    else:
+        return 16.
+
+def cloudimpact_73(bt_mod, bt_obs):
+    """
+    follows Harnisch 2016
+    """
+    biascor_obs = 0
+    bt_lim = 255  # Kelvin for 7.3 micron WV channel
+
+    ci_obs = max(0, bt_lim-(bt_obs - biascor_obs))
+    ci_mod = max(0, bt_lim-bt_mod)
+    ci = (ci_obs+ci_mod)/2
+    return ci
+
+def read_prior_obs(f_obs_prior):
+    """
+    docstring
+    """
+    obsseq = open(f_obs_prior, 'r').readlines()
+    OBSs = []
+    # read observations from obs_seq.final
+    for i, line in enumerate(obsseq):
+        if ' OBS ' in line:
+            observed = float(obsseq[i+1].strip())
+            truth = float(obsseq[i+2].strip())
+            prior_ensmean = float(obsseq[i+3].strip())
+            prior_enssd = float(obsseq[i+4].strip())
+            prior_ens = []
+            for j in range(5, 5+exp.n_ens):
+                prior_ens.append(float(obsseq[i+j].strip()))
+
+            OBSs.append(dict(observed=observed, truth=truth,
+                            prior_ens=np.array(prior_ens)))
+    return OBSs
+
+def edit_obserr_in_obsseq(fpath_obsseqin, OEs):
+    """
+    overwrite observation errors in a obs_seq.out file
+    according to the values in OEs
+    """
+    # write to txt (write whole obs_seq.out again)
+    obsseq = open(fpath_obsseqin, 'r').readlines()
+    obsseq_new = obsseq.copy()
+    i_obs = 0
+    for i, line in enumerate(obsseq):
+        if 'kind\n' in line:
+            i_line_oe = i+9  # 9 for satellite obs
+            obsseq_new[i_line_oe] = ' '+str(OEs[i_obs])+'   \n'
+            i_obs += 1
+
+    os.rename(fpath_obsseqin, fpath_obsseqin+'-bak')  # backup
+    # write cloud dependent errors (actually whole file)
+    with open(fpath_obsseqin, 'w') as f:
+        for line in obsseq_new:
+            f.write(line)
+
+def set_input_nml(sat_channel=False, just_prior_values=False):
+    """descr"""
+    if just_prior_values:
+        template = cluster.scriptsdir+'/../templates/input.prioronly.nml'
+    else:
+        template = cluster.scriptsdir+'/../templates/input.nml'
+    copy(template, cluster.dartrundir+'/input.nml')
+    sed_inplace(cluster.dartrundir+'/input.nml', '<n_ens>', str(int(exp.n_ens)))
+
+    # input.nml for RTTOV
+    if sat_channel > 0:
+        if sat_channel in [1, 2, 3, 12]:
+            rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.VIS.nml'
+        else:
+            rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
+        append_file(cluster.dartrundir+'/input.nml', rttov_nml)
+
+if __name__ == "__main__":
+
+    time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
+    print(dt.datetime.now())
+
+    # error_var = (5.)**2
+    # kind = 'RADAR'
+
+    # error_var = (0.5)**2
+    # kind = 'RASO_T'
+
+    #obs.generic_obs(kind, time, exp.n_obs, exp.radar_err**2, exp.distance_between_obs_meters,
+    #                output_path=cluster.dartrundir,
+    #                fpath_obs_locations=cluster.archivedir()+time.strftime('/%Y-%m-%d_%H:%M')
+    #                +'/obs_coords.pkl')
+
+
+    for channel_id in exp.sat_channels:
+
+        n_obs = 100
+        channel_id = 6
+        distance_between_obs_meters = 10000
+
+        error_var = (1.)**2
+        create_obsseq.sat(time, channel_id, n_obs, error_var, 
+                          distance_between_obs_meters,
+                          output_path=cluster.dartrundir, 
+                          fpath_obs_locations='./domain.pkl')
+        print(dt.datetime.now())
+
+        if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
+            raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
+
+        # generate observations
+        set_input_nml(sat_channel=channel_id, just_prior_values=True)
+        os.chdir(cluster.dartrundir)
+        t = dt.datetime.now()
+        os.system('mpirun -np 12 ./perfect_model_obs')
+        print('1st perfect_model_obs', (dt.datetime.now()-t).total_seconds())
+
+        if channel_id == 6:
+            # run ./filter to have prior observation estimates from model state
+            set_input_nml(sat_channel=channel_id, just_prior_values=True)
+            t = dt.datetime.now()
+            os.system('mv obs_seq.out obs_seq_all.out; mpirun -np 20 ./filter')
+            print('1st filter', (dt.datetime.now()-t).total_seconds())
+
+            # find the observation error for a pair of (H(x_nature), H(x_background))
+            f_obs_prior = cluster.dartrundir+'/obs_seq.final'
+            OBSs = read_prior_obs(f_obs_prior)
+
+            # compute the observation error necessary
+            # to achieve a certain operational FGD distribution
+            OEs = []
+            for obs in OBSs:
+                bt_y = obs['truth']
+                bt_x_ens = obs['prior_ens']
+                CIs = [cloudimpact_73(bt_x, bt_y) for bt_x in bt_x_ens]
+
+                oe_nature = oe_73(np.mean(CIs))
+                OEs.append(oe_nature)
+
+            # write obs_seq.out
+            fpath_obsseqout = cluster.dartrundir+'/obs_seq.in'
+            edit_obserr_in_obsseq(fpath_obsseqout, OEs)
+            print('after editing oe', dt.datetime.now())
+
+            # generate actual observations (with correct error)
+            os.chdir(cluster.dartrundir)
+            t = dt.datetime.now()
+            os.system('mpirun -np 12 ./perfect_model_obs')
+            print('real obs gen', (dt.datetime.now()-t).total_seconds())
+
+        # correct input.nml for actual assimilation later on
+        set_input_nml(sat_channel=channel_id, just_prior_values=False)
+
+    # FIXME: missing is the concatenation of
+    os.chdir(cluster.dartrundir)
+    os.system('cat obs_seq.out >> obs_seq_all.out')
+    print(dt.datetime.now())
+    # FIXME: what if different channels in one obs_seq.out -> need different input.nml for different channels

scripts/pre_assim.py

@@ -16,7 +16,7 @@ for iens in range(1, exp.n_ens+1):
 
     print('link wrfout file to DART background file')
     wrfout_run = exppath_firstguess+background_init_time.strftime('/%Y-%m-%d_%H:%M/')  \
-                 +str(iens)+'/'+assim_time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
+                 +str(iens)+assim_time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
     dart_ensdir = cluster.dartrundir+'/advance_temp'+str(iens)
     wrfout_dart = dart_ensdir+'/wrfout_d01'
 
@@ -49,5 +49,5 @@ os.system('rm -rf '+cluster.dartrundir+'/output_mean*')
 os.system('rm -rf '+cluster.dartrundir+'/output_sd*')
 os.system('rm -rf '+cluster.dartrundir+'/perfect_output_*')
 
-print('replace measurement error with obs error for assimilation')  # FIXME !!! temporary only
-sed_inplace(cluster.dartrundir+'/obs_seq.out', '9.000000000000000E-004', '0.04')
+#print('replace measurement error with obs error for assimilation')  # not needed anymore
+#sed_inplace(cluster.dartrundir+'/obs_seq.out', '9.000000000000000E-004', '0.04')

scripts/prepare_namelist.py

-import os, sys, shutil
+import os, sys, shutil, warnings
 import datetime as dt
 
 sys.path.append(os.getcwd())
 from config.cfg import exp, cluster
-from utils import sed_inplace, copy, symlink
+from utils import sed_inplace, copy, symlink, mkdir
 
 def run(cluster, iens, begin, end):
     rundir = cluster.wrf_rundir(iens)
@@ -26,6 +26,17 @@ def run(cluster, iens, begin, end):
                  '<HH2>': '%H', '<MM2>': '%M'}.items():
         sed_inplace(rundir+'/namelist.input', k, end.strftime(v))
 
+    #########################
+    try:
+        print('copy wrfinput of this run to archive')
+        wrfin_old = rundir+'/wrfinput_d01'
+        init_dir = cluster.archivedir()+begin.strftime('/%Y-%m-%d_%H:%M/')+str(iens)
+        os.makedirs(init_dir, exist_ok=True)
+        wrfin_arch = init_dir+'/wrfinput_d01'
+        copy(wrfin_old, wrfin_arch)
+    except Exception as e:
+        warnings.warn(str(e))
+
 
 if __name__ == '__main__':
     begin = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')

scripts/run_obs_diag.py

@@ -28,19 +28,25 @@ def run(folder_obs_seq_final):
     print('running obs_diag program')
     os.chdir(rundir_program)
     symlink(cluster.dart_srcdir+'/obs_diag', rundir_program+'/obs_diag')
-    os.system('./obs_diag >& obs_diag.log')
+    try:
+        os.remove(rundir_program+'/obs_seq_to_netcdf')
+    except:
+        pass
+    os.system('./obs_diag >& obs_diag.log')  # caution, this overwrites obs_seq_to_netcdf
 
     outdir = '/'.join(folder_obs_seq_final.split('/')[:-1])
     print('moving output to', outdir+'/obs_diag_output.nc')
     copy(rundir_program+'/obs_diag_output.nc', outdir+'/obs_diag_output.nc')
 
     print('running obs_seq_to_netcdf program')
+    shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', cluster.dart_srcdir+'/obs_seq_to_netcdf')
     symlink(cluster.dart_srcdir+'/obs_seq_to_netcdf', rundir_program+'/obs_seq_to_netcdf')
-    os.system('./obs_seq_to_netcdf  >& obs_seq_to_netcdf.log')
-    print('moving output to', outdir+'/obs_seq_output.nc')
-    copy(rundir_program+'/obs_diag_output.nc', outdir+'/obs_seq_output.nc')
+    os.system('./obs_seq_to_netcdf  >& obs_seq_to_netcdf.log')  # caution, overwrites its own binary?!
+    print('moving output to', outdir+'/obs_seq...')
+    os.system('mv '+rundir_program+'/obs_epoch_*.nc '+outdir+'/')
 
 
 if __name__ == '__main__':
-    folder_obs_seq_final = '/home/fs71386/lkugler/data/sim_archive/exp_v1.11_LMU_filter_domainobs/obs_seq_final/'
+    #folder_obs_seq_final = '/home/fs71386/lkugler/data/sim_archive/exp_v1.11_LMU_filter2/obs_seq_final/'
+    folder_obs_seq_final = str(sys.argv[1])
     run(folder_obs_seq_final)

scripts/update_wrfinput_from_filteroutput.py

 import os, sys, warnings
 import datetime as dt
+import netCDF4 as nc
 
 from config.cfg import exp, cluster
 from utils import symlink, copy_scp_srvx8, copy, mkdir, mkdir_srvx8, clean_wrfdir
 
 time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
+background_init_time = dt.datetime.strptime(sys.argv[2], '%Y-%m-%d_%H:%M')
+exppath_firstguess = str(sys.argv[3])
 
 #if cluster.name != 'srvx8':
 #    copy = copy_scp_srvx8
 #    mkdir = mkdir_srvx8
+cycle_vars = ['U', 'V', 'P', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW',
+              'QGRAUP', 'QNICE', 'QNRAIN', 'U10', 'V10', 'T2', 'Q2', 'PSFC', 'TSLB',
+              'SMOIS', 'TSK']
 
-update_vars = ['Times', 'U', 'V', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
+update_vars = ['Times', 'U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
+
+# variables which are updated need not to be cycled
+for var in update_vars:
+    if var in cycle_vars:
+        cycle_vars.remove(var)
+
+cycles = ','.join(cycle_vars)
+updates = ','.join(update_vars)
 
 print('move output to WRF dir as new initial conditions')
 for iens in range(1, exp.n_ens+1):
     clean_wrfdir(cluster.wrf_rundir(iens))
+    prior_wrf = exppath_firstguess + background_init_time.strftime('/%Y-%m-%d_%H:%M/') \
+              +str(iens)+time.strftime('/wrfout_d01_%Y-%m-%d_%H:%M:%S')
     filter_out = cluster.dartrundir+'/filter_restart_d01.'+str(iens).zfill(4)
     wrf_ic = cluster.wrf_rundir(iens) + '/wrfinput_d01'
 
-    # overwrite variables in wrfinput file
-    vars = ','.join(update_vars)
-    print('updating', vars, 'in', wrf_ic)
-    os.system(cluster.ncks+' -A -v '+vars+' '+filter_out+' '+wrf_ic)
+    # cycles variables from wrfout (prior state)
+    print('cycling', cycles, 'into', wrf_ic, 'from', prior_wrf)
+    os.system(cluster.ncks+' -A -v '+cycles+' '+prior_wrf+' '+wrf_ic)
+
+    print('updating', updates, 'in', wrf_ic, 'from', filter_out)
+    os.system(cluster.ncks+' -A -v '+updates+' '+filter_out+' '+wrf_ic)
+
+    print('writing T into THM of wrfinput')
+    thm_in = nc.Dataset(filter_out, 'r').variables['T'][:]
+    dsout = nc.Dataset(wrf_ic, 'r+')
+    dsout.variables['THM'][:] = thm_in
+    dsout.close()
 
     # clean up
     #try:

templates/input.nml

@@ -155,7 +155,8 @@
                                 'ACARS_V_WIND_COMPONENT',
                                 'ACARS_TEMPERATURE',
                                 'GPSRO_REFRACTIVITY',
-                                'LAND_SFC_ALTIMETER',
+                                'DOPPLER_RADIAL_VELOCITY',
+				'RADAR_REFLECTIVITY',
                                 'MSG_4_SEVIRI_RADIANCE',
 				'MSG_4_SEVIRI_TB',
 				'MSG_4_SEVIRI_BDRF'
@@ -177,7 +178,7 @@
    lowest_reflectivity_fwd_op  =  -10.0,
    max_radial_vel_obs          =   1000000,
    allow_wet_graupel           =   .false.,
-   microphysics_type           =       2  ,
+   microphysics_type           =       5  ,
    allow_dbztowt_conv          =   .false.,
    dielectric_factor           =  0.224,
    n0_rain                     =  8.0e6,
@@ -218,6 +219,7 @@
                          'PSFC',  'QTY_SURFACE_PRESSURE',     'TYPE_PSFC', 'UPDATE','999',
                          'T2',    'QTY_2M_TEMPERATURE',       'TYPE_T',    'UPDATE','999',
                          'TSK',   'QTY_SKIN_TEMPERATURE',     'TYPE_T',    'UPDATE','999',
+                         'REFL_10CM','QTY_RADAR_REFLECTIVITY','TYPE_REFL', 'UPDATE','999',
    wrf_state_bounds    = 'QVAPOR','0.0','NULL','CLAMP',
                          'QRAIN', '0.0','NULL','CLAMP',
                          'QCLOUD','0.0','NULL','CLAMP',
@@ -298,10 +300,10 @@
 &obs_diag_nml
    obs_sequence_name  = '',
    obs_sequence_list  = 'obsdiag_inputlist.txt',
-   first_bin_center   =  2008, 7,30,11, 0, 0 ,
-   last_bin_center    =  2008, 7,30,16, 0, 0 ,
-   bin_separation     =     0, 0, 0, 0,15, 0 ,
-   bin_width          =     0, 0, 0, 0,15, 0 ,
+   first_bin_center   =  2008, 7,30,10, 0, 0 ,
+   last_bin_center    =  2008, 7,30,18, 0, 0 ,
+   bin_separation     =     0, 0, 0, 0,30, 0 ,
+   bin_width          =     0, 0, 0, 0,30, 0 ,
    time_to_skip       =     0, 0, 0, 0, 0, 0 ,
    max_num_bins       = 1000,
    trusted_obs        = 'null',
@@ -314,7 +316,7 @@
    print_mismatched_locs = .false.,
    create_rank_histogram = .true.,
    outliers_in_histogram = .true.,
-   use_zero_error_obs    = .false.,
+   use_zero_error_obs    = .true.,
    verbose               = .false.
    /
 
@@ -332,7 +334,7 @@
 &obs_seq_to_netcdf_nml  
    obs_sequence_name = '',
    obs_sequence_list = 'obsdiag_inputlist.txt',
-   append_to_netcdf  = .false.,
+   append_to_netcdf  = .true.,
    lonlim1    =    0.0,
    lonlim2    =  360.0,
    latlim1    =  -90.0,

templates/namelist.input

@@ -90,6 +90,7 @@
  tke_adv_opt                         = 1,      1,      1,
  non_hydrostatic                     = .true., .true., .true.,
  mix_full_fields                     = .true., .true., .true.,
+ use_theta_m			     = 0,
  /
 
  &bdy_control