diff --git a/assimilation_code/modules/assimilation/filter_mod.f90 b/assimilation_code/modules/assimilation/filter_mod.f90 index 4db5c98a95..1df96a2f9c 100644 --- a/assimilation_code/modules/assimilation/filter_mod.f90 +++ b/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) diff --git a/models/wrf_hydro/create_identity_streamflow_obs.f90 b/models/wrf_hydro/create_identity_streamflow_obs.f90 index aeaa661996..2479f21e31 100644 --- a/models/wrf_hydro/create_identity_streamflow_obs.f90 +++ b/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 diff --git a/models/wrf_hydro/model_mod.f90 b/models/wrf_hydro/model_mod.f90 index 573c8a23dc..343871c021 100644 --- a/models/wrf_hydro/model_mod.f90 +++ b/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 !------------------------------------------------------------------