minor fixes

config/cfg.py

@@ -17,23 +17,30 @@ exp.n_nodes = 10
 
 n_obs = 1600  # radar: n_obs for each observation height level
 
-vis = dict(kind='MSG_4_SEVIRI_BDRF',
+vis = dict(plotname='VIS 0.6µm',  plotunits='[1]',
+           kind='MSG_4_SEVIRI_BDRF',
            sat_channel=1, n_obs=n_obs, err_std=0.03,
            cov_loc_radius_km=20)
-wv73 = dict(kind='MSG_4_SEVIRI_TB',
+wv73 = dict(plotname='Brightness temperature WV 7.3µm',  plotunits='[K]',
+            kind='MSG_4_SEVIRI_TB',
             sat_channel=6, n_obs=n_obs, err_std=False,
             cov_loc_radius_km=20)
-ir108 = dict(kind='MSG_4_SEVIRI_TB',
+ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
+             kind='MSG_4_SEVIRI_TB',
              sat_channel=9, n_obs=n_obs, err_std=5.,
              cov_loc_radius_km=20)
 
-radar = dict(kind='RADAR', n_obs=n_obs, err_std=5.,
+radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
+             kind='RADAR', 
+             n_obs=n_obs, err_std=5.,
              heights=np.arange(1000, 15001, 1000),
              cov_loc_radius_km=10, cov_loc_vert_km=2)
 
-t2m = dict(kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=1.0, 
+t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
+           kind='SYNOP_TEMPERATURE', n_obs=n_obs, err_std=1.0, 
            cov_loc_radius_km=32, cov_loc_vert_km=1)
-psfc = dict(kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
+psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
+            kind='SYNOP_SURFACE_PRESSURE', n_obs=n_obs, err_std=50.,
             cov_loc_radius_km=32, cov_loc_vert_km=5)
 
 

scheduler.py

@@ -58,7 +58,7 @@ def backup_scripts():
 
     # reproducibility
     for f in ['scheduler.py', 'config/clusters.py', 'config/cfg.py']:
-        func(current+'/../'+f,  main_a+'/scheduler.py', shutil.copy)
+        func(current+'/../'+f,  main_a+f, shutil.copy)
         
     for f in os.listdir(current):
         func(os.path.join(current, f), main_a+'/', shutil.copy)
@@ -230,7 +230,7 @@ timedelta_btw_assim = dt.timedelta(minutes=30)
 backup_scripts()
 id = None
 
-start_from_existing_state = True
+start_from_existing_state = False
 is_new_run = not start_from_existing_state
 
 if is_new_run:

scripts/apply_obs_op_dart.py

@@ -27,12 +27,10 @@ if __name__ == '__main__':
     pre_assim.run(time, prev_forecast_init, exppath_firstguess)
 
     savedir = cluster.archivedir()+'/obs_seq_final_1min/'
-    mkdir(savedir)
 
     n_stages = len(exp.observations)
     for istage, obscfg in enumerate(exp.observations):
 
-        kind = obscfg['kind']
         n_obs = obscfg['n_obs']
         sat_channel = obscfg.get('sat_channel', False)
         obscfg['folder_obs_coords'] = False
@@ -50,9 +48,10 @@ if __name__ == '__main__':
         wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc',
                            cluster.dartrundir+'/wrfout_d01')
         aso.run_perfect_model_obs()
-        aso.assimilate(nproc=48)
+        aso.assimilate(nproc=96)
 
-        archive_stage = savedir+kind
         # only the prior state values are of interest in this file
-        aso.archive_diagnostics(archive_stage, time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final'))
+        # observation and truth is wrong in this file (dummy)
+        archive_stage = savedir+'/assim_stage'+str(istage)
+        aso.archive_diagnostics(archive_stage, time)
 

scripts/assim_synth_obs.py

@@ -106,7 +106,7 @@ def set_DART_nml(sat_channel=False, cov_loc_radius_km=32, cov_loc_vert_km=False,
         rttov_nml = cluster.scriptsdir+'/../templates/obs_def_rttov.IR.nml'
         append_file(cluster.dartrundir+'/input.nml', rttov_nml)
 
-def obs_operator_ensemble():
+def obs_operator_ensemble(istage):
     # assumes that prior ensemble is already linked to advance_temp<i>/wrfout_d01
     print('running obs operator on ensemble forecast')
     os.chdir(cluster.dartrundir)
@@ -122,7 +122,8 @@ def obs_operator_ensemble():
             # 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')
             
-            # add geodata
+            # add geodata, if istage>0, wrfout is DART output (has coords)
+            if istage == 0:
                 wrfout_add_geo.run(cluster.dartrundir+'/geo_em.d01.nc', cluster.dartrundir+'/wrfout_d01')
 
             # run perfect_model obs (forward operator)
@@ -182,7 +183,7 @@ def run_perfect_model_obs():
     try_remove(cluster.dartrundir+'/obs_seq.out')
     if not os.path.exists(cluster.dartrundir+'/obs_seq.in'):
         raise RuntimeError('obs_seq.in does not exist in '+cluster.dartrundir)
-    os.system('mpirun -np 12 ./perfect_model_obs')
+    os.system('mpirun -np 12 ./perfect_model_obs > log.perfect_model_obs')
 
 def assimilate(nproc=96):
     print('running filter')
@@ -190,13 +191,14 @@ def assimilate(nproc=96):
     try_remove(cluster.dartrundir+'/obs_seq.final')
     if not os.path.exists(cluster.dartrundir+'/obs_seq.out'):
         raise RuntimeError('obs_seq.out does not exist in '+cluster.dartrundir)
-    os.system('mpirun -genv I_MPI_PIN_PROCESSOR_LIST=0-'+str(int(nproc)-1)+' -np '+str(int(nproc))+' ./filter')
+    os.system('mpirun -genv I_MPI_PIN_PROCESSOR_LIST=0-'+str(int(nproc)-1)+' -np '+str(int(nproc))+' ./filter > log.filter')
 
 def archive_diagnostics(archive_dir, time):
     print('archive obs space diagnostics')
     mkdir(archive_dir)
-    copy(cluster.dartrundir+'/obs_seq.final', 
-         archive_dir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final'))
+    fout = archive_dir+time.strftime('/%Y-%m-%d_%H:%M_obs_seq.final')
+    copy(cluster.dartrundir+'/obs_seq.final', fout)
+    print(fout, 'saved.')
 
     # try:  # what are regression diagnostics?!
     #     print('archive regression diagnostics')
@@ -258,8 +260,7 @@ if __name__ == "__main__":
     n_stages = len(exp.observations)
     for istage, obscfg in enumerate(exp.observations):
 
-        kind = obscfg['kind']
-        archive_stage = archive_time + '/assim_stage'+str(istage)+'_'+kind
+        archive_stage = archive_time + '/assim_stage'+str(istage)
         n_obs = obscfg['n_obs']
         sat_channel = obscfg.get('sat_channel', False)
         obscfg['folder_obs_coords'] = archive_stage+'/obs_coords.pkl'
@@ -276,7 +277,7 @@ if __name__ == "__main__":
             osq.create_obsseq_in(time, obscfg, zero_error=True)  # zero error to get truth vals
 
             Hx_nat = obs_operator_nature(time) 
-            Hx_prior = obs_operator_ensemble()  # files are already linked to DART directory
+            Hx_prior = obs_operator_ensemble(istage)  # files are already linked to DART directory
             
             obscfg['err_std'] = calc_obserr_WV73(Hx_nat, Hx_prior)
         else:
@@ -287,7 +288,7 @@ if __name__ == "__main__":
         run_perfect_model_obs()
 
         assimilate()
-        dir_obsseq = cluster.archivedir()+'/obs_seq_final/assim_stage'+str(istage)+'_'+kind
+        dir_obsseq = cluster.archivedir()+'/obs_seq_final/assim_stage'+str(istage)
         archive_diagnostics(dir_obsseq, time)
 
         if istage < n_stages-1:

scripts/run_obs_diag.py

@@ -45,11 +45,11 @@ def run_obsdiag(filepaths, f_out='./obsdiag.nc'):
         os.system('./obs_diag >& obs_diag.log')  # caution, this overwrites obs_seq_to_netcdf
 
         # move output to archive
-        #outdir = outdir #'/'.join(folder_obs_seq_final.split('/')[:-1])
+        #outdir = os.path.dirname(f_out)  #'/'.join(folder_obs_seq_final.split('/')[:-1])
         if obserr_iszero == '.true.':
-            fout = outdir+'/'+f_out[:-3]+'_wrt_truth.nc'   
+            fout = f_out[:-3]+'_wrt_truth.nc'   
         elif obserr_iszero == '.false.':
-            fout = outdir+'/'+f_out[:-3]+'_wrt_obs.nc' 
+            fout = f_out[:-3]+'_wrt_obs.nc' 
         shutil.move(rundir_program+'/obs_diag_output.nc', fout)
         print(fout, 'saved.')
 
@@ -59,6 +59,7 @@ def run_obs_seq_to_netcdf(filepaths, f_out='./obs_epoch.nc'):
     write_input_filelist(filepaths)
     print('------ running obs_seq_to_netcdf program')
     shutil.copy(cluster.dart_srcdir+'/obs_seq_to_netcdf-bak', rundir_program+'/obs_seq_to_netcdf')
+    os.chdir(rundir_program)
     os.system('./obs_seq_to_netcdf  >& obs_seq_to_netcdf.log')  # caution, overwrites its own binary?!
     shutil.move(rundir_program+'/obs_epoch_001.nc', f_out)
     print(f_out, 'saved.')

scripts/wrfout_add_geo.py

@@ -23,9 +23,9 @@ def run(geo_data_file, wrfout_file):
 
     for old, new in zip(fields_old, fields_new):
         print('moving old field', old, 'into new field', new)
-        print(geo_ds.variables[old][:].shape, wrfinp_ds.variables[new][:].shape)
+        #print(geo_ds.variables[old][:].shape, wrfinp_ds.variables[new][:].shape)
         wrfinp_ds.variables[new][:] = geo_ds.variables[old][:]
-        print(wrfinp_ds.variables[new][:])
+        #print(wrfinp_ds.variables[new][:])
 
     wrfinp_ds.close()
     geo_ds.close()

templates/input.nml

@@ -284,7 +284,7 @@
    layout = 1,
    tasks_per_node = 96
    communication_configuration = 1
-   debug                       = .true.
+   debug                       = .false.
    /
 
 &obs_def_gps_nml