hydro/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F

   1 ! This module defines and instantiates objects
   2 ! for a hybrid persistence levelpool type
   3 ! reservoir. The hybrid reservoir type
   4 ! inherits input and output types from the
   5 ! reservoir base module and calls instantiation
   6 ! of these into sub-objects. The hybrid
   7 ! reservoir type also points to types for
   8 ! hybrid properties and state and calls
   9 ! instantiation of these into sub-objects.
  10 ! A pointer to a levelpool reservoir object
  11 ! is also held in state, and this module
  12 ! instantiates that levelpool object. There
  13 ! is also a subroutine to run hybrid reservoir
  14 ! that is derived from the reservoir base
  15 ! type interface to run reservoir.
  16
  17 module module_persistence_levelpool_hybrid
  18
  19     use module_persistence_levelpool_hybrid_properties, only: hybrid_properties_interface
  20     use module_persistence_levelpool_hybrid_state, only: hybrid_state_interface
  21     use module_levelpool, only: levelpool
  22     use module_reservoir_utilities, only: modify_for_projected_storage, warn_negative_inflow, &
  23                                           create_hybrid_diagnostic_log_file, &
  24                                           log_hybrid_diagnostic_data
  25     use module_reservoir, only: reservoir, reservoir_input, reservoir_output
  26     use module_reservoir_read_timeslice_data, only: usgs_timeslice_data, usace_timeslice_data, timeslice_data_type
  27     use module_hydro_stop, only: HYDRO_stop
  28     use iso_fortran_env, only: int64
  29     implicit none
  30
  31     ! Extend/derive hybrid type from the abstract base
  32     ! type for reservoirs.
  33     type, extends(reservoir) :: persistence_levelpool_hybrid
  34
  35         ! Define pointers to sub-types / sub-objects to and
  36         ! held by a level pool reservoir object.
  37         type (hybrid_properties_interface), pointer :: properties => null()
  38         type (hybrid_state_interface), pointer :: state => null()
  39         type (timeslice_data_type), pointer :: timeslice_data => null()
  40
  41         logical :: pointer_allocation_guard = .false.
  42
  43     contains
  44
  45         procedure :: init => hybrid_init
  46         procedure :: destroy => hybrid_destroy
  47         procedure :: run => run_hybrid_reservoir
  48
  49     end type persistence_levelpool_hybrid
  50
  51 #ifndef NCEP_WCOSS
  52     integer, parameter :: log_warning = 6
  53 #else
  54     integer, parameter :: log_warning = 78
  55 #endif
  56
  57 contains
  58
  59     ! Hybrid Constructor
  60     subroutine hybrid_init(this, water_elevation,  &
  61         lake_area, weir_elevation, weir_coeffecient, &
  62         weir_length, dam_length, orifice_elevation, orifice_coefficient, &
  63         orifice_area, lake_max_water_elevation, initial_fractional_depth, &
  64         lake_number, reservoir_type, reservoir_parameter_file, start_date, &
  65         usgs_timeslice_path, usace_timeslice_path, observation_lookback_hours, &
  66         observation_update_time_interval_seconds)
  67         implicit none
  68         class(persistence_levelpool_hybrid), intent(inout) :: this ! object being initialized
  69         real,    intent(inout) :: water_elevation           ! meters AMSL
  70         real,    intent(in)    :: lake_area                         ! area of lake (km^2)
  71         real,    intent(in)    :: weir_elevation            ! bottom of weir elevation (meters AMSL)
  72         real,    intent(in)    :: weir_coeffecient          ! weir coefficient
  73         real,    intent(in)    :: weir_length               ! weir length (meters)
  74         real,    intent(in)    :: dam_length                ! dam length (meters)
  75         real,    intent(in)    :: orifice_elevation         ! orifice elevation (meters AMSL)
  76         real,    intent(in)    :: orifice_coefficient       ! orifice coefficient
  77         real,    intent(in)    :: orifice_area              ! orifice area (meters^2)
  78         real,    intent(in)    :: lake_max_water_elevation  ! max water elevation (meters)
  79         real,    intent(in)    :: initial_fractional_depth  ! initial fraction water depth
  80         integer(kind=int64), intent(in)    :: lake_number               ! lake number
  81         integer, intent(in)    :: reservoir_type            ! reservoir type
  82         character(len=*),   intent(in) :: reservoir_parameter_file
  83         character(len=19),  intent(in) :: start_date
  84         character(len=256), intent(in) :: usgs_timeslice_path
  85         character(len=256), intent(in) :: usace_timeslice_path
  86         integer,            intent(in) :: observation_lookback_hours
  87         integer,            intent(in) :: observation_update_time_interval_seconds
  88         character(len=15)              :: lake_number_string
  89         character(len=15)              :: reservoir_type_string
  90
  91
  92 #ifdef RESERVOIR_D
  93         ! Create diagnostic log file only for development/debugging purposes
  94         call create_hybrid_diagnostic_log_file(lake_number)
  95 #endif
  96
  97         if (this%pointer_allocation_guard .eqv. .false. ) then
  98
  99             ! Call to initialize timeslice data object. This object is a singleton for either
 100             ! USGS (U.S. Geological Survey) or USACE (U.S. Army Corps of Engineers) data source,
 101             ! so if another reservoir has already initialized the object of a particular data
 102             ! source, it will return immediately.
 103             if (reservoir_type == 2) then
 104                 call usgs_timeslice_data%init(start_date, usgs_timeslice_path, &
 105                 reservoir_parameter_file, "usgs", observation_lookback_hours)
 106
 107                 ! Assign pointer to usgs_timeslice_data object
 108                 this%timeslice_data=>usgs_timeslice_data
 109
 110             else if (reservoir_type == 3) then
 111                 call usace_timeslice_data%init(start_date, usace_timeslice_path, &
 112                 reservoir_parameter_file, "usace", observation_lookback_hours)
 113
 114                 ! Assign pointer to usace_timeslice_data object
 115                 this%timeslice_data=>usace_timeslice_data
 116
 117             else
 118                 ! Call hydro_stop if reservoir_type is not USGS or USACE because the
 119                 ! data source is not available.
 120                 write(lake_number_string, "(I15)") lake_number
 121                 write(reservoir_type_string, "(I15)") reservoir_type
 122                 call hydro_stop("ERROR: Incorrect reservoir type for reservoir " // trim(ADJUSTL(lake_number_string)) // &
 123                 ". Expected reservoir type 2 or 3 and instead received reservoir type "  &
 124                 // trim(ADJUSTL(reservoir_type_string)) // "." )
 125             end if
 126
 127             ! try to allocate input
 128             allocate ( this%input )
 129             if ( .not. associated(this%input) ) then
 130                 ! if the input structure could not be created, call hydro_stop.
 131                 write(lake_number_string, "(I15)") lake_number
 132                 call hydro_stop("ERROR: Failure to allocate hybrid input structure for reservoir " &
 133                 // trim(ADJUSTL(lake_number_string)) // ".")
 134             else
 135                 ! initialize the input structure
 136                 call this%input%init()
 137             end if
 138
 139             ! try to allocate output
 140             allocate ( this%output )
 141             if ( .not. associated(this%output) ) then
 142                 ! if the output structure could not be created, call hydro_stop.
 143                 write(lake_number_string, "(I15)") lake_number
 144                 call hydro_stop("ERROR: Failure to allocate hybrid output structure for reservoir " &
 145                 // trim(ADJUSTL(lake_number_string)) // ".")
 146             else
 147                 ! initialize the output structure
 148                 call this%output%init()
 149             end if
 150
 151             ! try to allocate properties
 152             allocate ( this%properties )
 153             if ( .not. associated(this%properties) ) then
 154                 ! if the properties structure could not be created, call hydro_stop.
 155                 write(lake_number_string, "(I15)") lake_number
 156                 call hydro_stop("ERROR: Failure to allocate hybrid properties structure for reservoir " &
 157                 // trim(ADJUSTL(lake_number_string)) // ".")
 158             else
 159                 ! initialize hybrid properties
 160                 call this%properties%init(lake_area, lake_max_water_elevation, orifice_elevation, lake_number, &
 161                 reservoir_type, observation_lookback_hours, observation_update_time_interval_seconds, &
 162                 reservoir_parameter_file)
 163             end if
 164             this%pointer_allocation_guard = .true.
 165
 166             ! try to allocate state
 167             allocate ( this%state )
 168             if ( .not. associated(this%state) ) then
 169                 ! if the state structure could not be created, call hydro_stop.
 170                 write(lake_number_string, "(I15)") lake_number
 171                 call hydro_stop("ERROR: Failure to allocate hybrid state structure for reservoir " &
 172                 // trim(ADJUSTL(lake_number_string)) // ".")
 173             else
 174                 ! initialize hybrid state
 175                 call this%state%init(water_elevation, lake_area, lake_max_water_elevation, orifice_elevation, &
 176                 initial_fractional_depth, reservoir_type)
 177             end if
 178             this%pointer_allocation_guard = .true.
 179
 180             ! Initialize persistence weight index to be out of bounds of the persistence weights array
 181             this%state%persistence_weight_index = SIZE(this%properties%persistence_weighted_coefficients) + 1
 182
 183             ! Output warning if the initial storage is greater than max storage or less than zero.
 184             if (this%state%current_storage > this%properties%max_storage) then
 185                 this%state%current_storage = this%properties%max_storage
 186                 write(lake_number_string, "(I15)") lake_number
 187                 write(log_warning,*) "WARNING: initial storage exceeds max storage for reservoir ", lake_number_string, &
 188                 ". Setting initial storage to max storage."
 189
 190             else if (this%state%current_storage < this%properties%min_storage) then
 191                 this%state%current_storage = this%properties%min_storage
 192                 write(lake_number_string, "(I15)") lake_number
 193                 write(log_warning,*) "WARNING: initial storage is less than zero for reservoir ", lake_number_string, &
 194                 ". Setting initial storage to zero."
 195
 196             end if
 197
 198             ! Allocate a single level pool reservoir
 199             allocate(levelpool :: this%state%levelpool_ptr)
 200
 201             ! Initialize level pool reservoir
 202             call this%state%levelpool_ptr%init(water_elevation, lake_area, &
 203             weir_elevation, weir_coeffecient, weir_length, dam_length, orifice_elevation, &
 204             orifice_coefficient, orifice_area, lake_max_water_elevation, lake_number)
 205
 206         end if
 207     end subroutine hybrid_init
 208
 209
 210     ! Hybrid Destructor
 211     subroutine hybrid_destroy(this)
 212         implicit none
 213         class(persistence_levelpool_hybrid), intent(inout) :: this ! object being destroyed
 214     end subroutine hybrid_destroy
 215
 216
 217     ! Subroutine for running reservoir for a hybrid reservoir
 218     subroutine run_hybrid_reservoir(this, previous_timestep_inflow, inflow, &
 219         lateral_inflow, water_elevation, outflow, routing_period, dynamic_reservoir_type, &
 220         assimilated_value, assimilated_source_file)
 221         implicit none
 222         class(persistence_levelpool_hybrid), intent(inout) :: this
 223         real, intent(in)    :: previous_timestep_inflow ! cubic meters per second (cms)
 224         real, intent(in)    :: inflow                   ! cubic meters per second (cms)
 225         real, intent(in)    :: lateral_inflow           ! cubic meters per second (cms)
 226         real, intent(inout) :: water_elevation          ! meters
 227         real, intent(out)   :: outflow                  ! cubic meters per second (cms)
 228         real, intent(in)    :: routing_period           ! seconds
 229         integer, intent(out):: dynamic_reservoir_type   ! dynamic reservoir type sent to lake out files
 230         real, intent(out)   :: assimilated_value        ! value assimilated from observation
 231         character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value
 232         real                :: delta_storage            ! timestep change in storage (cubic meters)
 233         real                :: local_water_elevation    ! water elevation passed to levelpool (meters AMSL)
 234         real                :: levelpool_outflow        ! cubic meters per second (cms)
 235         logical             :: max_storage_reached      ! flag for when max storage is reached
 236         integer             :: gage_lookback_seconds    ! seconds between current model
 237                                                         ! time and the time a gage
 238                                                         ! discharge is read
 239
 240         max_storage_reached = .false.
 241         gage_lookback_seconds = 0
 242
 243         ! Update input variables
 244         this%input%inflow = inflow
 245         this%input%lateral_inflow = lateral_inflow
 246
 247         ! The initialization and management of water elevation might diverge between the global state of the model
 248         ! and the internal state. The global state water elevation is passed to/from this subroutine at every
 249         ! timestep. This check is to help guard against an external piece of code adjusting the global state water
 250         ! elevation, which would cause it to no longer be synchronized with internal state water elevation. This
 251         ! check also allows this module to borrow from the existing restart capability by initializing reservoirs
 252         ! to some default water elevation and storage and then waiting for the first timestep to receive the
 253         ! initial condition, so that the reservoir would not need to manage its own restart. If water elevation
 254         ! is refactored to no longer be a global state, then this check is not necessary.
 255         if (this%state%water_elevation .ne. water_elevation) then
 256
 257             ! Update state variables
 258             this%state%water_elevation = water_elevation
 259
 260             this%state%current_storage = (this%state%water_elevation - this%properties%orifice_elevation) * &
 261             this%properties%lake_area
 262
 263         end if
 264
 265         local_water_elevation = this%state%water_elevation
 266
 267         ! If update time to read new timeslice gage discharges. In forecast modes, the update time automatically is set to
 268         ! 1,000,000 to ensure that gage discharges are only read at the first timestep. For analysis runs, a timeslice
 269         ! update time is set as a namelist parameter. A typical update time for an analysis run would be 3600 seconds
 270         ! (1 hour), at which time new observations would be retrieved. The first reservoir for each processor to reach an
 271         ! update time will call the function to read a new timeslice file, and these timeslice discharges will be held in
 272         ! memory for the subsequent reservoirs to use at that timestep.
 273         if (this%state%current_time >= this%state%timeslice_update_time) then
 274
 275             ! Read new timeslice gage discharges
 276             call this%timeslice_data%setup_read_timeslice(this%properties%observation_update_time_interval_seconds, &
 277             this%state%current_time, this%properties%gage_id, gage_lookback_seconds, this%state%gage_discharge, &
 278             this%state%assimilated_source_file)
 279
 280             ! If no good quality gage discharge was found in the given lookback period, the gage discharge is returned
 281             ! as -1.0. The persistence weight index will be set out of bounds larger than the array. This causes the
 282             ! the weight update logic to set the persistence weight to 0.0. Therefore, only levelpool calculations will be used.
 283             if (this%state%gage_discharge < 0.0) then
 284
 285                 ! If at weight update time, then increment persistence weight index.
 286                 ! Otherwise, continue using previous persistence weight index.
 287                 if (this%state%current_time >= this%state%weight_update_time) then
 288                     ! Increment weight index
 289                     this%state%persistence_weight_index = this%state%persistence_weight_index + 1
 290
 291                     ! Set next weight update time
 292                     this%state%weight_update_time = this%state%weight_update_time + this%properties%weight_update_time_interval
 293                 end if
 294
 295                 ! If out of bounds of array, then set persistence weight to 0.0. Otherwise, set persistence weight from array's indexed value.
 296                 if (this%state%persistence_weight_index > SIZE(this%properties%persistence_weighted_coefficients)) then
 297                     this%state%persistence_current_weight = 0.0
 298
 299                 else
 300                     this%state%persistence_current_weight = this%properties%persistence_weighted_coefficients(this%state%persistence_weight_index)
 301                 end if
 302
 303             else
 304
 305                 ! Set persisted outflow to gage discharge
 306                 this%state%persisted_outflow = this%state%gage_discharge
 307
 308                 ! Set assimilated value to gage discharge
 309                 this%state%assimilated_value = this%state%gage_discharge
 310
 311                 ! Start persistence weight index at 1
 312                 this%state%persistence_weight_index = 1
 313
 314                 ! Grab first persistence weight at index 1 from array
 315                 this%state%persistence_current_weight = this%properties%persistence_weighted_coefficients(this%state%persistence_weight_index)
 316
 317                 ! Set weight update time offset by how long back a timeslice discharge was read. For instance, if the weight update interval is
 318                 ! 24 hours, and the gage discharge was read 2 hours back, then the weight update time will be set to 22 hours after the current time.
 319                 this%state%weight_update_time = this%state%current_time + this%properties%weight_update_time_interval - gage_lookback_seconds
 320
 321             end if
 322
 323             ! Calculate levelpool weight
 324             this%state%levelpool_current_weight = 1.0 - this%state%persistence_current_weight
 325
 326             ! Set timeslice update time
 327             this%state%timeslice_update_time = this%state%timeslice_update_time + this%properties%observation_update_time_interval_seconds
 328
 329         ! If update time to change persistence weights
 330         else if (this%state%current_time >= this%state%weight_update_time) then
 331
 332             ! Increment weight index
 333             this%state%persistence_weight_index = this%state%persistence_weight_index + 1
 334
 335             ! Boundary check to not exceed the size of the persistence weights array
 336             if (this%state%persistence_weight_index <= SIZE(this%properties%persistence_weighted_coefficients)) then
 337
 338                 ! Grab indexed persistence weight from array
 339                 this%state%persistence_current_weight = this%properties%persistence_weighted_coefficients(this%state%persistence_weight_index)
 340
 341             else
 342                 ! If boundary of persistence weights array has been exceeded, then set all persistence weights to 0.0
 343                 this%state%persistence_current_weight = 0.0
 344
 345             end if
 346
 347             ! Calculate levelpool weight
 348             this%state%levelpool_current_weight = 1.0 - this%state%persistence_current_weight
 349
 350             ! Set next weight update time
 351             this%state%weight_update_time = this%state%weight_update_time + this%properties%weight_update_time_interval
 352
 353         end if
 354
 355         ! Referencing issue https://github.com/NCAR/wrf_hydro_nwm_public/issues/326 for the uncertainty about
 356         ! previous_timestep_inflow. Further understanding of the exact timing of inflow and previous timestep
 357         ! inflow might affect which values are passed to this calling of levelpool and subsequent mass/balance
 358         ! calculations.
 359         ! Run levelpool reservoir
 360         call this%state%levelpool_ptr%run(previous_timestep_inflow, inflow, &
 361         lateral_inflow, local_water_elevation, levelpool_outflow, routing_period, this%state%levelpool_reservoir_type, &
 362         this%state%levelpool_assimilated_value, this%state%levelpool_assimilated_source_file)
 363
 364         ! If the levelpool weight is within epsilon of 1.0
 365         if (this%state%levelpool_current_weight .ge. 1.0 - epsilon(1.0)) then
 366
 367             ! Set dynamic_reservoir_type to levelpool type
 368             this%state%dynamic_reservoir_type = this%state%levelpool_reservoir_type
 369
 370             ! Set the assimilated_value to sentinel, -9999.0
 371             this%state%assimilated_value = -9999.0
 372
 373             ! Set the assimilated_source_file to empty string
 374             this%state%assimilated_source_file = ""
 375         else
 376
 377             ! Set dynamic_reservoir_type to given USGS or USACE type
 378             this%state%dynamic_reservoir_type = this%properties%reservoir_type
 379         end if
 380
 381         ! Calculate outflow weighted between persistence and levelpool
 382         this%output%outflow = this%state%persistence_current_weight &
 383         * this%state%persisted_outflow + this%state%levelpool_current_weight &
 384         * levelpool_outflow
 385
 386         ! Warn if there is a negative inflow
 387         call warn_negative_inflow(this%input%inflow, this%properties%lake_number, this%state%current_time)
 388
 389         ! Modify if exceeding storage boundary conditions
 390         call modify_for_projected_storage(this%input%inflow, &
 391         this%state%current_storage, this%properties%min_storage, &
 392         this%properties%max_storage, this%properties%lake_number, &
 393         this%state%current_time, int(routing_period), this%output%outflow, max_storage_reached)
 394
 395         ! If max storage has been reached, release greater of persistence or levelpool
 396         if (max_storage_reached .and. this%output%outflow < levelpool_outflow) then
 397             this%output%outflow = levelpool_outflow
 398
 399             ! Set dynamic_reservoir_type to levelpool type
 400             this%state%dynamic_reservoir_type = this%state%levelpool_reservoir_type
 401
 402             ! Set the assimilated_value to sentinel, -9999.0
 403             this%state%assimilated_value = -9999.0
 404
 405             ! Set the assimilated_source_file to empty string
 406             this%state%assimilated_source_file = ""
 407         end if
 408
 409         ! Calculate change in storage
 410         delta_storage = (this%input%inflow - this%output%outflow) * routing_period
 411
 412         ! Update storage from the most recent model states
 413         this%state%current_storage = this%state%current_storage + delta_storage
 414
 415         ! Calculate new water elevation. Delta storage is used as opposed to calculating from current storage in
 416         ! order to minimize floating point error because current storage is a much larger magnitude than water elevation.
 417         this%state%water_elevation = this%state%water_elevation + delta_storage / this%properties%lake_area
 418
 419         ! Update output variable returned from this subroutine
 420         outflow = this%output%outflow
 421
 422         ! Set current inflow to previous_timestep_inflow
 423         this%input%previous_timestep_inflow = inflow
 424
 425         ! Update water_elevation variable returned from this subroutine
 426         water_elevation = this%state%water_elevation
 427
 428         ! Update the dynamic_reservoir_type returned from this subroutine
 429         dynamic_reservoir_type = this%state%dynamic_reservoir_type
 430
 431         ! Update the assimilated_value returned from this subroutine
 432         assimilated_value = this%state%assimilated_value
 433
 434         ! Update the assimilated_source_file returned from this subroutine
 435         assimilated_source_file = this%state%assimilated_source_file
 436
 437         ! Update the current time
 438         this%state%current_time = this%state%current_time + int(routing_period)
 439
 440 #ifdef RESERVOIR_D
 441         ! Log diagnostic data only for development/debugging purposes
 442         call log_hybrid_diagnostic_data(this%properties%lake_number, this%state%current_time, &
 443         this%state%timeslice_update_time, this%state%weight_update_time, gage_lookback_seconds, &
 444         this%properties%gage_id, this%state%persistence_weight_index, &
 445         this%state%persistence_current_weight, this%input%inflow, this%state%current_storage, &
 446         this%state%water_elevation, this%state%gage_discharge, this%state%persisted_outflow, &
 447         levelpool_outflow, this%output%outflow)
 448 #endif
 449
 450     end subroutine run_hybrid_reservoir
 451
 452 end module module_persistence_levelpool_hybrid