Merge branch 'hybrid_hydrodart_exp' into hybrid

assimilation_code/modules/assimilation/filter_mod.f90

@@ -666,11 +666,19 @@ subroutine filter_main()
 
 ! for now, assume that we only allow cycling if single_file_out is true.
 ! code in this call needs to know how to initialize the output files.
-call initialize_file_information(num_state_ens_copies ,                     &
+if (do_hybrid) then 
+   call initialize_file_information(num_state_ens_copies ,                  &
                                  file_info_input      , file_info_mean_sd,  &
                                  file_info_forecast   , file_info_preassim, &
                                  file_info_postassim  , file_info_analysis, &
                                  file_info_output     , file_info_hybrid)
+else
+   call initialize_file_information(num_state_ens_copies ,                  &
+                                 file_info_input      , file_info_mean_sd,  &
+                                 file_info_forecast   , file_info_preassim, &
+                                 file_info_postassim  , file_info_analysis, &
+                                 file_info_output     )
+endif
 
 call check_file_info_variable_shape(file_info_output, state_ens_handle)
 

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
@@ -209,9 +209,18 @@ 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)
+! MEG 6/29/23 [Commented out 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)
 
+print *, gages_list_file
+print *, 'n_wanted_gages: ', n_wanted_gages
+
 num_new_obs = estimate_total_obs_count(input_files, nfiles)
 
 inquire(file=output_file, exist=file_exist)
@@ -316,12 +325,21 @@ program create_identity_streamflow_obs
       if (indx == 0) cycle OBSLOOP
 
       ! relate the physical location to the dart state vector index
-      dart_index = linkloc_to_dart(lat(indx), lon(indx))
+      dart_index = linkloc_to_dart(lat(indx), lon(indx)) 
+
+      ! desired gauges get the provided obs_err
+      ! remaining ones are for verification purposes
+      if (ANY(desired_gages == station_strings(n))) then  
+        oerr = max(discharge(n)*obs_fraction_for_error, MIN_OBS_ERR_STD)
+      else 
+        oerr = MAX_OBS_ERR_STD
+      endif
 
+      ! MEG: 6/29/23 [Commented out, oerr is conditionally computed now]
       ! 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)
+      !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 

models/wrf_hydro/model_mod.f90

@@ -26,7 +26,8 @@ module model_mod
 
 use  netcdf_utilities_mod, only : nc_check, nc_add_global_attribute, &
                                   nc_synchronize_file, nc_end_define_mode, &
-                                  nc_add_global_creation_time, nc_begin_define_mode
+                                  nc_add_global_creation_time, nc_begin_define_mode, &
+                                  nc_get_dimension_size, nc_open_file_readonly
 
 use obs_def_utilities_mod, only : track_status
 
@@ -536,6 +537,7 @@ function read_model_time(filename)
 integer :: DimID, VarID, strlen, ntimes
 logical :: isLsmFile, isClimFile
 integer :: ncid, io
+integer :: c_link
 
 io = nf90_open(filename, NF90_NOWRITE, ncid)
 call nc_check(io,routine,'open',filename)
@@ -587,6 +589,21 @@ function read_model_time(filename)
    allocate(datestring(ntimes))
    datestring(1) = '1980-01-01_00:00:00'
 
+   ! Also check if the state in the climatology is consistent 
+   ! with the state in the restarts
+   ncid   = nc_open_file_readonly(filename, routine)
+   c_link = nc_get_dimension_size(ncid, 'links', routine)
+
+   if ( c_link /= n_link ) then
+      write(string1,'(A)')'The size of the state in the climatology files is not consistent with the current domain size.'
+      write(string2, *   )'number of links: ', c_link, &
+                      ' from "'//trim(filename)//'"'
+      write(string3,*)'number of links: ',int(n_link,i8), &
+                      ' from "'//get_hydro_domain_filename()//'"'
+      call error_handler(E_ERR, routine, string1, &
+                 source, revision, revdate, text2=string2, text3=string3)
+   endif
+
 else ! Get the time from the hydro or parameter file
 
    io = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'Restart_Time', len=strlen)
@@ -612,6 +629,7 @@ function read_model_time(filename)
 
 deallocate(datestring)
 
+
 end function read_model_time
 
 !------------------------------------------------------------------