hydro/Routing/Reservoirs/module_reservoir_read_timeslice_data.F90

   1 ! This module reads USGS (U.S. Geological Survey) or USACE
   2 ! (U.S. Army Corps of Engineers) timeslice files to get gage discharge
   3 ! values that will be used by reservoirs. An observation lookback
   4 ! period is passed in to determine how far back in time from the
   5 ! current model time the module will look for timeslice files.
   6 ! The observation resolution determines the time increments the
   7 ! module will look back. For instance, a standard lookback period
   8 ! would be 18 hours with an observation resolution of 15 minutes,
   9 ! where a model current time of 8:00 PM would search for timeslice
  10 ! files at every 15 minute increment between 2:00 AM and 8:00 PM
  11 ! that day. The module will first search for the most recent
  12 ! timeslice files and grab the discharge for a particular
  13 ! lake/reservoir if the gage quality standard is met at that time.
  14 ! If a gage discharge is missing or if the gage quality standard
  15 ! is not met for any particular lake/reservoir in the given
  16 ! timeslice file, the module will continue to look back at every
  17 ! observation resolution increment until either all
  18 ! lakes/reservoirs have a good quality discharge or the end of
  19 ! the lookback period is reached. The total lookback seconds
  20 ! from current model time that the discharge is read will also
  21 ! be returned.
  22
  23 module module_reservoir_read_timeslice_data
  24     use module_reservoir_utilities, only: get_timeslice_array_dimension, &
  25                                           read_timeslice_gage_ids, &
  26                                           read_timeslice_netcdf_real_1D_variables, &
  27                                           read_persistence_netcdf_gage_ids_all, &
  28                                           geth_newdate, &
  29                                           handle_err
  30     use netcdf
  31     implicit none
  32
  33     ! The timeslice_data is a singleton in that the object is instantiated only once per
  34     ! processor by the first reservoir on that processor that needs it. The object is then used
  35     ! by every other reservoir on that processor that needs it.
  36     type :: timeslice_data_type
  37
  38         character(len=19)                             :: start_date
  39         character(len=19)                             :: current_date
  40         character(len=256)                            :: timeslice_path
  41         character(len=15)                             :: data_source
  42         integer                                       :: number_of_reservoir_gages
  43         integer                                       :: observation_lookback_hours
  44         integer                                       :: current_lookback_seconds
  45         integer                                       :: last_update_time_seconds
  46         integer,            allocatable, dimension(:) :: gage_lookback_seconds
  47         character(len=15),  allocatable, dimension(:) :: reservoir_gage_ids
  48         real,               allocatable, dimension(:) :: reservoir_gage_discharges
  49         character(len=256), allocatable, dimension(:) :: timeslice_file_names
  50         character(len=15),  allocatable, dimension(:) :: gage_id
  51         real,               allocatable, dimension(:) :: gage_discharge
  52         real,               allocatable, dimension(:) :: gage_quality
  53         logical                                       :: initialized = .FALSE.
  54
  55     contains
  56
  57         procedure :: init => timeslice_data_init
  58         procedure :: destroy => timeslice_data_destroy
  59         procedure :: setup_read_timeslice => setup_read_timeslice
  60         procedure :: read_timeslice_file => read_timeslice_file
  61
  62     end type timeslice_data_type
  63
  64     type (timeslice_data_type), target :: usgs_timeslice_data, usace_timeslice_data
  65
  66     integer, parameter :: observation_resolution_minutes = 15
  67     real,    parameter :: gage_quality_threshold = 0.9
  68
  69 contains
  70
  71     ! Timeslice Data Type Constructor
  72     subroutine timeslice_data_init(this, start_date, timeslice_path, &
  73         reservoir_parameter_file, data_source, observation_lookback_hours)
  74         implicit none
  75         class(timeslice_data_type), intent(inout) :: this ! object being initialized
  76         character(len=19),  intent(in)                   :: start_date
  77         character(len=256), intent(in)                   :: timeslice_path
  78         character(len=*),   intent(in)                   :: reservoir_parameter_file
  79         character(len=*),   intent(in)                   :: data_source
  80         integer,            intent(in)                   :: observation_lookback_hours
  81         integer                                          :: ncid, var_id
  82         integer                                          :: status  ! status of reading NetCDF
  83         integer                                          :: timeslice_gage_index
  84         character(len=15), allocatable, dimension(:)     :: gage_ids_string_array
  85         character(len=15)                                :: gage_id_string, gage_id_string_trimmed
  86         integer                                          :: gage_id_integer
  87         integer                                          :: char_index, number_counter, temp_val
  88         character(len=1)                                 :: temp_char
  89
  90         ! Return from subroutine if this singleton is already initialized
  91         if (this%initialized) return
  92         this%initialized = .true.
  93
  94         this%data_source = ADJUSTL(trim(data_source))
  95
  96         status = nf90_open(path = reservoir_parameter_file, mode = nf90_nowrite, ncid = ncid)
  97         if (status /= nf90_noerr) call handle_err(status, "Could not open reservoir parameter file " &
  98         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
  99
 100         call read_persistence_netcdf_gage_ids_all(ncid, ADJUSTL(trim(this%data_source)) // "_gage_id", &
 101         reservoir_parameter_file, var_id, this%number_of_reservoir_gages)
 102
 103         allocate (this%reservoir_gage_ids(this%number_of_reservoir_gages))
 104         allocate (gage_ids_string_array(this%number_of_reservoir_gages))
 105
 106         status = nf90_get_var(ncid, var_id, gage_ids_string_array)
 107         if (status /= nf90_noerr) call handle_err(status, "Error reading " // ADJUSTL(trim(this%data_source)) &
 108          // "_gage_id from " // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 109
 110         status = nf90_close(ncid)
 111         if (status /= nf90_noerr) call handle_err(status, "Could not close reservoir parameter file " &
 112         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 113
 114         ! Convert gage ids to integers
 115         do timeslice_gage_index = 1, this%number_of_reservoir_gages
 116             gage_id_string = gage_ids_string_array(timeslice_gage_index)
 117
 118             do char_index = 1, 15
 119                 temp_val = ichar(gage_id_string(char_index:char_index))
 120
 121                 ! Check if null character, then pad with a space
 122                 if (temp_val == 0) then
 123                     gage_id_string(char_index:char_index) = ' '
 124                 end if
 125             end do
 126
 127             gage_id_string_trimmed = ADJUSTL(trim(gage_id_string))
 128
 129             this%reservoir_gage_ids(timeslice_gage_index) = gage_id_string_trimmed
 130
 131         end do
 132
 133         if(allocated(gage_ids_string_array)) deallocate(gage_ids_string_array)
 134
 135         this%start_date = start_date
 136         this%current_date = start_date
 137         this%timeslice_path = timeslice_path
 138         this%observation_lookback_hours = observation_lookback_hours
 139
 140         this%last_update_time_seconds = -999
 141
 142         ! Allocate lookback, dishcharge, and timeslice name arrays
 143         allocate (this%gage_lookback_seconds(this%number_of_reservoir_gages))
 144         allocate (this%reservoir_gage_discharges(this%number_of_reservoir_gages))
 145         allocate(this%timeslice_file_names(this%number_of_reservoir_gages))
 146
 147     end subroutine timeslice_data_init
 148
 149
 150     ! Timeslice Data Type Destructor
 151     subroutine timeslice_data_destroy(this)
 152
 153         implicit none
 154         class(timeslice_data_type), intent(inout) :: this ! object being destroyed
 155
 156     end subroutine timeslice_data_destroy
 157
 158
 159     ! Set up list of timeslice files to read
 160     ! The timeslices are read only once on each processor at each timeslice update time
 161     subroutine setup_read_timeslice(this, update_interval_seconds, current_reservoir_time_seconds, &
 162         reservoir_gage_id, gage_lookback_seconds, reservoir_gage_discharge, reservoir_timeslice_file_name)
 163         implicit none
 164
 165         class(timeslice_data_type), intent(inout) :: this
 166         integer, intent(in)  :: update_interval_seconds
 167         integer, intent(in)  :: current_reservoir_time_seconds
 168         character(len=*), intent(in)  :: reservoir_gage_id
 169         integer, intent(out) :: gage_lookback_seconds
 170         real, intent(out)    :: reservoir_gage_discharge
 171         character(len=256), intent(out)   :: reservoir_timeslice_file_name
 172         character(len=256)   :: timeslice_file_name
 173         character(len=256)   :: timeslice_file_name_full
 174         character(len=19)    :: old_date, new_date
 175         character(len=2)     :: observation_resolution_string
 176         integer :: total_timeslice_time_periods, timeslice_file_index, reservoir_subset_index, lake_index
 177         integer :: observation_minute, observation_resolution_seconds
 178         logical :: file_exists
 179
 180         ! This check ensures that the timeslice is read only once for each processor at each timeslice update time
 181         if (current_reservoir_time_seconds .ne. this%last_update_time_seconds) then
 182
 183             ! Checks if not the first timestep, then update date.
 184             if (this%last_update_time_seconds .ne. -999) then
 185                 call geth_newdate(this%current_date, this%current_date, update_interval_seconds)
 186             end if
 187
 188             this%last_update_time_seconds = current_reservoir_time_seconds
 189
 190             ! Set gage discharges of subset lakes to default -1.0
 191             this%reservoir_gage_discharges = -1.0
 192
 193             this%gage_lookback_seconds = 0
 194             this%current_lookback_seconds = 0
 195
 196             ! Set timeslice file names to empty strings by default
 197             this%timeslice_file_names = ""
 198
 199             total_timeslice_time_periods = this%observation_lookback_hours * (60 / observation_resolution_minutes)
 200
 201             old_date = this%current_date
 202
 203             read (old_date(15:16), *) observation_minute
 204
 205             ! Match minutes to proper observation resolution interval
 206             observation_minute = observation_minute / observation_resolution_minutes * observation_resolution_minutes
 207
 208             ! Formatting for writing minutes
 209             if (observation_minute < 10) then
 210                 old_date(15:16) = "00"
 211             else
 212                 write(old_date(15:16), "(I2)") observation_minute
 213             end if
 214
 215             ! Zero out seconds
 216             old_date(18:19) = "00"
 217
 218             write(observation_resolution_string, "(I2)") observation_resolution_minutes
 219
 220             ! Negative for going back in time
 221             observation_resolution_seconds = observation_resolution_minutes * 60 * (-1)
 222
 223             ! Loop through the total timeslice periods to look for and read timeslice files
 224             do timeslice_file_index = 1, total_timeslice_time_periods
 225
 226                 ! Use below for timeslice files with colons in the time
 227                 old_date = old_date(:13) // ':' // old_date(15:16) // ':' // old_date(18:)
 228
 229                 ! If any gages that still have discharges equal to -1.0, read the previous
 230                 ! time timeslice file available.
 231                 if ( any(this%reservoir_gage_discharges == -1.0)) then
 232
 233                     ! Construct timeslice filename
 234                     timeslice_file_name = old_date // "." // &
 235                             observation_resolution_string // 'min.' // ADJUSTL(trim(this%data_source)) // 'TimeSlice.ncdf'
 236
 237                     ! Add path to timeslice filename
 238                     timeslice_file_name_full = trim(this%timeslice_path) // "/" // ADJUSTL(trim(timeslice_file_name))
 239
 240                     ! Check if file exists
 241                     inquire(FILE = timeslice_file_name_full, EXIST = file_exists)
 242                     if (file_exists) then
 243                         ! Call subroutine to read a particular timeslice file
 244                         call this%read_timeslice_file(timeslice_file_name_full, timeslice_file_name)
 245                     end if
 246
 247                     ! Call subroutine to get the date from one observation resolution back in time
 248                     call geth_newdate(new_date, old_date, observation_resolution_seconds)
 249                     old_date = new_date
 250
 251                     ! Lookback seconds to account for how long back a timeslice file was read
 252                     this%current_lookback_seconds = this%current_lookback_seconds - observation_resolution_seconds
 253
 254                 else
 255                     exit
 256                 end if
 257             end do
 258
 259         end if
 260
 261         ! Loop through gages to return gage lookback seconds and gage discharge for a given reservoir
 262         do reservoir_subset_index = 1, this%number_of_reservoir_gages
 263             if (this%reservoir_gage_ids(reservoir_subset_index) .eq. reservoir_gage_id) then
 264                 gage_lookback_seconds = this%gage_lookback_seconds(reservoir_subset_index)
 265                 reservoir_gage_discharge = this%reservoir_gage_discharges(reservoir_subset_index)
 266                 reservoir_timeslice_file_name = this%timeslice_file_names(reservoir_subset_index)
 267             end if
 268         end do
 269
 270     end subroutine setup_read_timeslice
 271
 272
 273     ! Read given timeslice file to get gage discharges
 274     subroutine read_timeslice_file(this, timeslice_file_full, timeslice_file)
 275         implicit none
 276         class(timeslice_data_type), intent(inout) :: this
 277         character(len=256), intent(in) :: timeslice_file_full
 278         character(len=256), intent(in) :: timeslice_file
 279         integer :: reservoir_gage_index, timeslice_gage_index, number_of_gages
 280         integer :: ncid, status
 281
 282         ! Open Timeslice NetCDF file
 283         status = nf90_open(path = trim(timeslice_file_full), mode = nf90_nowrite, ncid = ncid)
 284         if (status /= nf90_noerr) call handle_err(status, "Could not open timeslice file " &
 285         // trim(ADJUSTL(timeslice_file_full)) // ".")
 286
 287         ! Get dimension of gage arrays
 288         call get_timeslice_array_dimension(ncid, 'discharge', timeslice_file_full, number_of_gages)
 289
 290         ! Allocate gage info arrays
 291         allocate(this%gage_id(number_of_gages))
 292         allocate(this%gage_discharge(number_of_gages))
 293         allocate(this%gage_quality(number_of_gages))
 294
 295         ! Get gage ids
 296         call read_timeslice_gage_ids(ncid, 'stationId', timeslice_file_full, this%gage_id)
 297
 298         ! Get gage discharges
 299         call read_timeslice_netcdf_real_1D_variables(ncid, 'discharge', timeslice_file_full, this%gage_discharge)
 300
 301         ! Get gage qualities
 302         call read_timeslice_netcdf_real_1D_variables(ncid, 'discharge_quality', timeslice_file_full, this%gage_quality)
 303
 304         ! Normalize gage qualities to 1
 305         this%gage_quality = this%gage_quality/100
 306
 307         ! First, quality check on the quality flag
 308         where(this%gage_quality .lt. 0 .or. this%gage_quality .gt. 1) this%gage_quality=0
 309
 310         ! Currently using line below that would set quality to 0 if there is a negative
 311         ! discharge but there may actually be a deficit in flow, i.e. no flow available
 312         ! at all that the transducer has given us a negative pressure. We do not account
 313         ! for that case here but may need to consider this later.
 314         where(this%gage_discharge .lt. 0.000) this%gage_quality=0
 315
 316         ! Currently not using line below instead of above line that sets quality to 0 for negative discharge
 317         !where(this%gage_discharge .le. 0.000) this%gage_discharge = 0.0
 318
 319         !! peak flow on MS river *2
 320         !! http://nwis.waterdata.usgs.gov/nwis/peak?site_no=07374000&agency_cd=USGS&format=html
 321         !! baton rouge 1945: 1,473,000cfs=41,711cms, multiply it roughly by 2
 322         where(this%gage_discharge .ge. 90000.0) this%gage_quality=0
 323
 324         ! Loop through reservoir gage array to look for gages that still have discharges set to -1.0.
 325         ! Match gage ids from the timeslice file gage id arrays to gage ids in the reservoir gage subset
 326         ! array. If the matched gage is good quality, then set that reservoir gage subset discharge to
 327         ! the discharge read from the timeslice file along with the total lookback seconds from current
 328         ! model time that the gage discharge is read from a timeslice file.
 329         do reservoir_gage_index = 1, this%number_of_reservoir_gages
 330             if (this%reservoir_gage_discharges(reservoir_gage_index) == -1.0) then
 331                 do timeslice_gage_index = 1, number_of_gages
 332                     if (ADJUSTL(trim(this%gage_id(timeslice_gage_index))) &
 333                     .eq. ADJUSTL(trim(this%reservoir_gage_ids(reservoir_gage_index)))) then
 334                         if (this%gage_quality(timeslice_gage_index) .gt. gage_quality_threshold) then
 335                             this%reservoir_gage_discharges(reservoir_gage_index) = this%gage_discharge(timeslice_gage_index)
 336                             this%gage_lookback_seconds(reservoir_gage_index) = this%current_lookback_seconds
 337                             this%timeslice_file_names(reservoir_gage_index) = ADJUSTL(trim(timeslice_file))
 338                         end if
 339                         exit
 340                     end if
 341                 end do
 342             end if
 343         end do
 344
 345     ! Close timeslice NetCDF file
 346     status = nf90_close(ncid)
 347     if (status /= nf90_noerr) call handle_err(status, "Could not close timeslice file " // trim(ADJUSTL(timeslice_file_full)) // ".")
 348
 349     ! Deallocate gage arrays
 350     if(allocated(this%gage_id)) deallocate(this%gage_id)
 351     if(allocated(this%gage_discharge)) deallocate(this%gage_discharge)
 352     if(allocated(this%gage_quality)) deallocate(this%gage_quality)
 353
 354     end subroutine read_timeslice_file
 355
 356 end module module_reservoir_read_timeslice_data