Merge pull request #638 from mgharamti/hydroDART

Recent WRF-Hydro Developments

Author: mjs2369

models/wrf_hydro/create_identity_streamflow_obs.f90

@@ -66,7 +66,7 @@ program create_identity_streamflow_obs
 integer,  parameter :: NUM_COPIES      = 1      ! number of copies in sequence
 integer,  parameter :: NUM_QC          = 1      ! number of QC entries
 real(r8), parameter :: MIN_OBS_ERR_STD = 0.1_r8 ! m^3/sec
-real(r8), parameter :: MAX_OBS_ERR_STD = 100000.0_r8 
+real(r8), parameter :: MAX_OBS_ERR_STD = 1000000.0_r8 
 real(r8), parameter :: NORMAL_FLOW     = 10.0_r8
 real(r8), parameter :: contract        = 0.001_r8
 
@@ -104,7 +104,7 @@ program create_identity_streamflow_obs
 real(r8),                       allocatable :: discharge(:)
 
 character(len=IDLength),        allocatable :: desired_gages(:)
-integer         :: n_wanted_gages
+integer         :: n_wanted_gages, n_desired_gages
 real(r8)        :: oerr, qc
 integer         :: oday, osec
 type(obs_type)  :: obs
@@ -127,6 +127,7 @@ program create_identity_streamflow_obs
 character(len=256) :: location_file   = 'location.nc'
 character(len=256) :: gages_list_file = ''
 real(r8)           :: obs_fraction_for_error = 0.01
+logical            :: assimilate_all  = .false. 
 integer            :: debug = 0
 
 namelist / create_identity_streamflow_obs_nml / &
@@ -135,6 +136,7 @@ program create_identity_streamflow_obs
                location_file, &
                gages_list_file, &
                obs_fraction_for_error, &
+               assimilate_all, &
                debug
 
 !-------------------------------------------------------------------------------
@@ -209,7 +211,12 @@ program create_identity_streamflow_obs
 call init_obs(obs,      num_copies=NUM_COPIES, num_qc=NUM_QC)
 call init_obs(prev_obs, num_copies=NUM_COPIES, num_qc=NUM_QC)
 
-n_wanted_gages = set_desired_gages(gages_list_file)
+! Collect all the gauges: 
+! - desired ones will have the provided obs_err_sd
+! - remaining gauges are dummy with very large obs_err_sd
+
+n_desired_gages = set_desired_gages(gages_list_file)
+n_wanted_gages  = 0 !set_desired_gages(gages_list_file)
 call find_textfile_dims(input_files, nfiles)
 
 num_new_obs = estimate_total_obs_count(input_files, nfiles)
@@ -308,7 +315,8 @@ program create_identity_streamflow_obs
 
    OBSLOOP: do n = 1, nobs
 
-      if ( discharge(n) < 0.0_r8 ) cycle OBSLOOP
+      ! make sure discharge is physical 
+      if ( discharge(n) < 0.0_r8 .or. discharge(n) /= discharge(n) ) cycle OBSLOOP
 
       ! relate the TimeSlice:station to the RouteLink:gage so we can
       ! determine the location
@@ -318,13 +326,13 @@ program create_identity_streamflow_obs
       ! relate the physical location to the dart state vector index
       dart_index = linkloc_to_dart(lat(indx), lon(indx))
 
-      ! oerr is the observation error standard deviation in this application.
-      ! The observation error variance encoded in the observation file
-      ! will be oerr*oerr
-      oerr = max(discharge(n)*obs_fraction_for_error, MIN_OBS_ERR_STD)
-
-      ! MEG: A fix to not crush the ensemble in a no-flood period (stagnant water).  
-      !if ( discharge(n) < NORMAL_FLOW ) then 
+      ! desired gauges get the provided obs_err
+      ! remaining ones are for verification purposes
+      if (ANY(desired_gages == station_strings(n)) .or. assimilate_all) then  
+        oerr = max(discharge(n)*obs_fraction_for_error, MIN_OBS_ERR_STD)
+      else 
+        oerr = MAX_OBS_ERR_STD
+      endif
          ! don't correct that much, the gauge observations imply that the flow 
          ! in the stream is small. This is not a flood period. Streamflow values
          ! indicate a more or less lake situation rather than a strongly flowing stream. 
@@ -660,9 +668,9 @@ function estimate_total_obs_count(file_list,nfiles) result (num_obs)
 
 ! We need to know how many observations there may be.
 ! Specifying too many is not really a problem.
-! I am adding 20%
+! I am multiplying by 10.
 
-num_obs = 1.2_r8 * nobs * nfiles
+num_obs = 10.0_r8 * nobs * nfiles
 
 end function estimate_total_obs_count
 

models/wrf_hydro/hydro_dart_py/hydrodartpy/core/create_usgs_daily_obs_seq.py

@@ -50,7 +50,7 @@ def parallel_process_day(arg_dict):
 
     the_cmd = exe_cmd.format(
         **{
-            'cmd': './' + the_converter.name,
+            'cmd': './create_identity_streamflow_obs',
             'nproc': 1
         }
     )

models/wrf_hydro/hydro_dart_py/hydrodartpy/core/setup_usgs_daily.py

@@ -27,7 +27,8 @@ def setup_usgs_daily(
     input_dir = usgs_daily_config['input_dir']
     output_dir = usgs_daily_config['output_dir']
     # Output directory: make if DNE
-    output_dir.mkdir(exist_ok=False, parents=True)
+    #output_dir.mkdir(exist_ok=False, parents=True)
+    output_dir.mkdir(exist_ok=True, parents=True)
 
     # converter: identity or regular obs converter?
     # Check that the desired obs converter is in the dart build
@@ -75,7 +76,8 @@ def setup_usgs_daily(
     run_dir = config['experiment']['run_dir']
     m0 = pickle.load(open(run_dir / "member_000/WrfHydroSim.pkl", 'rb'))
     route_link_f = run_dir / 'member_000' / m0.base_hydro_namelist['hydro_nlist']['route_link_f']
-    (output_dir / route_link_f.name).symlink_to(route_link_f)
+    if not route_link_f.is_file():
+        (output_dir / route_link_f.name).symlink_to(route_link_f)
     input_nml[converter_nml]['location_file'] = route_link_f.name
 
     #input.nml input_files: create a list of files in the start and end range.
@@ -101,27 +103,35 @@ def setup_usgs_daily(
     if usgs_daily_config['identity_obs']:
 
         hydro_rst_file = run_dir / 'member_000' / m0.base_hydro_namelist['hydro_nlist']['restart_file']
-        (output_dir / hydro_rst_file.name).symlink_to(hydro_rst_file)
+        if not hydro_rst_file.is_file():
+            (output_dir / hydro_rst_file.name).symlink_to(hydro_rst_file)
         input_nml['model_nml']['domain_order'] = 'hydro'
         input_nml['model_nml']['domain_shapefiles'] = str(hydro_rst_file.name)
 
         f90nml.Namelist(m0.base_hydro_namelist).write(output_dir / 'hydro.namelist', force=True)
         top_level_dir = get_top_level_dir_from_config(config, m0)
-        (output_dir / top_level_dir).symlink_to(config['wrf_hydro']['domain_src'] / top_level_dir)
+        nwm_dir = config['wrf_hydro']['domain_src'] / top_level_dir
+        if not nwm_dir.is_dir():
+            (output_dir / top_level_dir).symlink_to(config['wrf_hydro']['domain_src'] / top_level_dir)
 
     # Now we are done editing it, write the input.nml back out.
-    input_nml.write(output_dir / 'input.nml')
+    out_input = output_dir / 'input.nml'
+    if not out_input.is_file():
+        input_nml.write(output_dir / 'input.nml')
 
     # Symlink the config file into the output_dir so the default yaml file name
     # can be used by create_usgs_daily_obs_seq.
     if config_file is None:
         config_file = sorted(exp_dir.glob('original.*.yaml'))[0]
-    (output_dir / 'config_file.yaml').symlink_to(config_file)
+    if not config_file.is_file():
+        (output_dir / 'config_file.yaml').symlink_to(config_file)
 
     # Stage the file that does the batch processing.
     this_file = pathlib.Path(__file__)
     batcher_base = 'create_usgs_daily_obs_seq.py'
-    (output_dir / batcher_base).symlink_to(this_file.parent / batcher_base)
+    pyscript = this_file.parent / batcher_base
+    if not pyscript.is_file():
+        (output_dir / batcher_base).symlink_to(this_file.parent / batcher_base)
 
     # Setup the scheduled script.
     orig_submit_script = this_file.parent / 'submission_scripts/submit_usgs_daily_obs_converter.sh'
@@ -152,10 +162,11 @@ def setup_usgs_daily(
 
         # Select statement
         # Right now, only single node processing
-        select_stmt = 'select=1:ncpus={ncpus}:mpiprocs={mpiprocs}'.format(
+        select_stmt = 'select=1:ncpus={ncpus}:mpiprocs={mpiprocs}:mem={reqmem}GB'.format(
             **{
                 'ncpus': usgs_sched['ncpus'],
-                'mpiprocs': usgs_sched['mpiprocs']
+                'mpiprocs': usgs_sched['mpiprocs'],
+                'reqmem': usgs_sched['reqmem']
             }
         )
         replace_in_file(this_submit_script, 'PBS_SELECT_TEMPLATE', select_stmt)
@@ -191,6 +202,7 @@ def setup_usgs_daily(
     all_obs_dir = pathlib.PosixPath(config['observation_preparation']['all_obs_dir'])
     all_obs_seq = output_dir.glob('obs_seq.*')
     for oo in all_obs_seq:
-        (all_obs_dir / oo.name).symlink_to(oo)
+        if not oo.is_file():
+            (all_obs_dir / oo.name).symlink_to(oo)
 
     return 0