hydro/Routing/Reservoirs/module_reservoir_utilities.F

   1 ! This module contains a variety of subroutines/functions
   2 ! used by reservoir objects.
   3
   4 module module_reservoir_utilities
   5     use netcdf
   6     use module_hydro_stop, only: HYDRO_stop
   7     use iso_fortran_env, only: int64
   8     implicit none
   9
  10 #ifndef NCEP_WCOSS
  11     integer, parameter :: log_warning = 6
  12 #else
  13     integer, parameter :: log_warning = 78
  14 #endif
  15
  16 contains
  17
  18     ! Checks and warns for boundary conditions of negative inflow
  19     subroutine warn_negative_inflow(inflow, lake_number, current_time)
  20         real, intent(in)    :: inflow
  21         integer, intent(in) :: lake_number, current_time
  22
  23         if (inflow < 0.0) then
  24             write(log_warning,*) "WARNING: Current inflow ", inflow, &
  25             " cubic meters per second has reached below zero for reservoir ", &
  26             lake_number, " at ", current_time, " seconds after model start time."
  27         end if
  28
  29     end subroutine warn_negative_inflow
  30
  31
  32     ! Checks and warns for boundary conditions of exceeding max or min storage
  33     subroutine warn_storage_boundaries(inflow, end_time_storage, min_storage, max_storage, lake_number, current_time)
  34         real, intent(in)    :: inflow, end_time_storage, min_storage, max_storage
  35         integer, intent(in) :: lake_number, current_time
  36
  37         ! Also have option to set the end_time_storage to the max_storage and change the outflow accordingly.
  38         if (end_time_storage > max_storage) then
  39
  40             write(log_warning,*) "WARNING: Calculated storage for this timestep, ", end_time_storage, &
  41             " cubic meters, will exceed the maximum storage boundary of: ", max_storage, &
  42             " cubic meters, for reservoir ", lake_number," at ", current_time, " seconds after model start time."
  43             write(log_warning,*) "Current inflow is ", inflow, " cubic meters per second."
  44
  45         ! Also have option to set the end_time_storage to the min_storage and change the outflow accordingly.
  46         else if (end_time_storage < min_storage .and. end_time_storage > 0.0) then
  47
  48             write(log_warning,*) "WARNING: Calculated storage for this timestep, ", end_time_storage, &
  49             " cubic meters, will fall below the minimum storage requirement of: ", min_storage, &
  50             " cubic meters, for reservoir ", lake_number, " at ", current_time, " seconds after model start time."
  51             write(log_warning,*) "Current inflow is ", inflow, " cubic meters per second."
  52
  53         end if
  54
  55     end subroutine warn_storage_boundaries
  56
  57
  58     ! Checks and modifies for boundary conditions of negative storage
  59     subroutine modify_negative_storage(inflow, previous_time_storage, routing_period, lake_number, current_time, &
  60         update_time, end_time_storage, outflow)
  61         real, intent(in)    :: inflow, previous_time_storage, routing_period
  62         integer, intent(in) :: lake_number, current_time, update_time
  63         real, intent(inout) :: end_time_storage, outflow
  64
  65         ! Note to user:
  66         ! In order to conserve mass, when a reservoir's storage would become 0.0 or less, the release
  67         ! must be adjusted to avoid releasing water that is not available in the reservoir. To determine
  68         ! conditions when this could happen, this function calculates the projected by the end of the
  69         ! timestep (dt) based on the latest known conditions of the storage, inflow, and release values.
  70
  71         ! If the current conditions would cause the storage to become less than or equal to 0.0 cubic meters,
  72         ! then the outflow/release will be altered to provide the available storage over the timestep
  73         ! (given as R = (S + I*dt)/dt ). This will cause storage at the end of the timestep to be 0.0. From
  74         ! this point, the outflow/release is bounded by the inflow provided until the next update time when
  75         ! the machine learning model will be run again.
  76
  77         ! Therefore, the user's choice of time_interval has a significant impact on this simulated outflow/release
  78         ! behavior in this edge case.  Release cannot increase again until the next update time occurs. This does,
  79         ! however, allow the reservoir to fill back up once inflow increases to exceed current outflow/release
  80         ! since storage would then become > 0.0.
  81
  82         ! If one wants to allow the release to continue to rise as inflow rises, you can set R = I*dt once storage
  83         ! has been depleted. This passes on inflow to outflow until an update time occurs, in which case the ML
  84         ! model and its post conditions are responsible for calculating an appropriate release.
  85
  86         ! The logical 'and' below is not necessary because if otherwise this check updates outflow and
  87         ! end_time_storage before an update time, those variables would be overwritten after a new release
  88         ! is calculated from the machine learning model. The logical 'and' is left here for emphasis on
  89         ! the importance it has on timesteps that are not also update times.
  90         if (end_time_storage <= 0.0 .AND. current_time < update_time) then
  91
  92             write(log_warning,*) "WARNING: End time storage, ", end_time_storage, &
  93             " cubic meters, has reached at or below 0.0 cubic meters for reservoir ", lake_number, &
  94             " at ", current_time, " seconds after model start time."
  95             write(log_warning,*) "Current inflow is ", inflow, " cubic meters per second."
  96
  97             outflow = (previous_time_storage + inflow * routing_period) / routing_period
  98
  99             end_time_storage = 0.0
 100
 101             if (outflow < 0.0) outflow = 0.0
 102
 103         end if
 104
 105     end subroutine modify_negative_storage
 106
 107
 108     ! Checks and modifies for boundary conditions of negative inflow
 109     subroutine modify_negative_inflow(inflow)
 110         real, intent(inout) :: inflow
 111
 112         if (inflow < 0.0) inflow = 0.0
 113
 114     end subroutine modify_negative_inflow
 115
 116
 117     ! Checks and modifies for boundary conditions of negative outflow or exceeding min or max storage
 118     subroutine modify_for_projected_storage(inflow, previous_time_storage, min_storage, max_storage, &
 119         lake_number, current_time, time_interval, outflow, max_storage_reached)
 120         real, intent(in)        :: inflow, previous_time_storage, min_storage, max_storage
 121         integer, intent(in)     :: lake_number, current_time, time_interval
 122         real, intent(inout)     :: outflow
 123         logical, intent(inout)  :: max_storage_reached
 124         real                    :: excess, projected_next_time_interval_storage
 125
 126         if (outflow < 0.0) then
 127             write(log_warning,*) "WARNING: Calculations return a negative outflow for reservoir ", lake_number
 128             write(log_warning,*) "at ", current_time, " seconds after model start time."
 129             outflow = 0.0
 130         end if
 131
 132         ! The projected_next_time_interval_storage is only calculated as a boundary condition to update the release
 133         ! and this storage is calculated assuming a constant flow rate over the time interval. If the time interval
 134         ! is too large relative to the channel routing period, this may not be an appropriate calculation.
 135         projected_next_time_interval_storage = previous_time_storage + (inflow - outflow) * time_interval
 136
 137         if (projected_next_time_interval_storage > max_storage) then
 138
 139             max_storage_reached = .true.
 140
 141             write(log_warning,*) "WARNING: Modified release to prevent maximum storage exceedance", &
 142             " for reservoir ", lake_number
 143             write(log_warning,*) "at ", current_time, " seconds after model start time."
 144
 145         else if (projected_next_time_interval_storage < min_storage &
 146             .AND. projected_next_time_interval_storage > 0.0) then
 147
 148             outflow = (projected_next_time_interval_storage - min_storage) / time_interval + outflow
 149
 150             write(log_warning,*) "WARNING: Modified release to maintain minimum storage for reservoir ", lake_number
 151             write(log_warning,*) "at ", current_time, " seconds after model start time."
 152
 153         else if (projected_next_time_interval_storage <= 0.0) then
 154
 155             outflow = inflow
 156
 157             write(log_warning,*)  "WARNING: Modified release to prevent storage deficit for reservoir ", lake_number
 158             write(log_warning,*) "at ", current_time, " seconds after model start time."
 159
 160         end if
 161
 162         if (outflow < 0.0) outflow = 0.0
 163
 164     end subroutine modify_for_projected_storage
 165
 166     ! Reads the reservoir_type form the reservoir parameter file. This subroutine first checks to ensure that
 167     ! the lake_id array in the reservoir parameter file matches the lake_id array in the lake parameter file and
 168     ! calls hydro_stop if they do not match.
 169     subroutine read_reservoir_type(reservoir_parameter_file, lake_parameter_lake_ids, number_of_lake_parameter_lakes, &
 170         reservoir_types)
 171         character(len=*), intent(in) :: reservoir_parameter_file
 172         integer(kind=int64), dimension(:), intent(in) :: lake_parameter_lake_ids
 173         integer, intent(in) :: number_of_lake_parameter_lakes
 174         integer, dimension(:), intent(inout) :: reservoir_types
 175         integer(kind=int64), allocatable, dimension(:) :: reservoir_parameter_lake_ids
 176
 177         integer, dimension(nf90_max_var_dims) :: dim_ids
 178         integer :: ncid, var_id, lake_index, number_of_reservoirs, number_of_reservoir_types
 179         integer :: status                        ! status of reading NetCDF
 180
 181         ! Open Reservoir Parameter NetCDF file
 182         status = nf90_open(path = reservoir_parameter_file, mode = nf90_nowrite, ncid = ncid)
 183         if (status /= nf90_noerr) call handle_err(status, "Could not open reservoir parameter file " &
 184         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 185
 186         status = nf90_inq_varid(ncid, "lake_id", var_id)
 187         if (status /= nf90_noerr) call handle_err(status, "lake_id read error in reservoir parameter file " &
 188         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 189
 190         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 191         if(status /= nf90_NoErr) call handle_err(status, "lake_id read error in reservoir parameter file " &
 192         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 193
 194         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_reservoirs)
 195         if(status /= nf90_NoErr) call handle_err(status, "lake_id read error in reservoir parameter file " &
 196         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 197
 198         allocate(reservoir_parameter_lake_ids(number_of_reservoirs))
 199
 200         status = nf90_get_var(ncid, var_id, reservoir_parameter_lake_ids)
 201         if (status /= nf90_noerr) call handle_err(status, "lake_id read error in reservoir parameter file " &
 202         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 203
 204         ! Check to ensure that the lake_id array size in the reservoir parameter file matches the lake_id array
 205         ! size in the lake parameter file
 206         if (number_of_lake_parameter_lakes .ne. number_of_reservoirs) then
 207             call hydro_stop("ERROR: Reservoir ID array in the reservoir parameter file does not match the Lake ID array &
 208             in the lake parameter file.")
 209
 210         else
 211             ! Check to ensure that each lake id in the lake_id array in the reservoir parameter file matches
 212             ! each lake id in the lake_id array in the lake parameter file
 213             do lake_index = 1, number_of_reservoirs
 214
 215                 if (lake_parameter_lake_ids(lake_index) .ne. reservoir_parameter_lake_ids(lake_index)) then
 216                     call hydro_stop("Reservoir ID array in the reservoir parameter file does not match the Lake ID array &
 217                     in the lake parameter file.")
 218                 end if
 219             end do
 220         end if
 221
 222         status = nf90_inq_varid(ncid, "reservoir_type", var_id)
 223         if (status /= nf90_noerr) call handle_err(status, "reservoir_type read error in reservoir parameter file " &
 224         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 225
 226         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 227         if(status /= nf90_NoErr) call handle_err(status, "reservoir_type read error in reservoir parameter file " &
 228         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 229
 230         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_reservoir_types)
 231         if(status /= nf90_NoErr) call handle_err(status, "reservoir_type read error in reservoir parameter file " &
 232         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 233
 234         ! Check to ensure that the reservoir_type array size matches the lake_id array size.
 235         if (number_of_reservoir_types .ne. number_of_lake_parameter_lakes) then
 236             call hydro_stop("ERROR: Reservoir Type array size in the reservoir parameter file does not match the Lake ID array size &
 237             in the lake parameter file.")
 238         end if
 239
 240         status = nf90_get_var(ncid, var_id, reservoir_types)
 241         if (status /= nf90_noerr) call handle_err(status, "reservoir_type read error in reservoir parameter file " &
 242         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 243
 244         ! Close Reservoir Parameter NetCDF file
 245         status = nf90_close(ncid)
 246         if (status /= nf90_noerr) call handle_err(status, "Could not close reservoir parameter file " &
 247         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 248
 249         if(allocated(reservoir_parameter_lake_ids)) deallocate(reservoir_parameter_lake_ids)
 250
 251     end subroutine read_reservoir_type
 252
 253
 254     ! Gets the dimensions of the gage arrays from the timeslice file
 255     subroutine get_timeslice_array_dimension(ncid, netcdf_array_name, timeslice_file, number_of_gages)
 256         integer, intent(in) :: ncid
 257         character(len=*), intent(in) :: netcdf_array_name
 258         character(len=256), intent(in) :: timeslice_file
 259         integer, intent(out) :: number_of_gages
 260         integer :: var_id, status
 261         integer, dimension(nf90_max_var_dims) :: dim_ids
 262
 263         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 264         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 265         // trim(ADJUSTL(timeslice_file)) // ".")
 266
 267         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 268         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 269         // trim(ADJUSTL(timeslice_file)) // ".")
 270
 271         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_gages)
 272         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 273         // trim(ADJUSTL(timeslice_file)) // ".")
 274
 275     end subroutine get_timeslice_array_dimension
 276
 277
 278     ! Read gage id array from the timeslice NetCDF
 279     subroutine read_timeslice_gage_ids(ncid, netcdf_array_name, timeslice_file, gage_id_var_array)
 280         integer, intent(in) :: ncid
 281         character(len=*), intent(in) :: netcdf_array_name
 282         character(len=256), intent(in) :: timeslice_file
 283         character(len=*), dimension(:), intent(out) :: gage_id_var_array
 284         integer :: var_id, status
 285
 286         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 287         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 288         // trim(ADJUSTL(timeslice_file)) // ".")
 289
 290         status = nf90_get_var(ncid, var_id, gage_id_var_array)
 291         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 292         // trim(ADJUSTL(timeslice_file)) // ".")
 293
 294     end subroutine read_timeslice_gage_ids
 295
 296
 297     ! Read real one dimension parameter arrays from the timeslice NetCDF
 298     subroutine read_timeslice_netcdf_real_1D_variables(ncid, netcdf_array_name, timeslice_file, real_var_array)
 299         integer, intent(in) :: ncid
 300         character(len=*), intent(in) :: netcdf_array_name
 301         character(len=256), intent(in) :: timeslice_file
 302         real, dimension(:), intent(out) :: real_var_array
 303         integer :: var_id, status
 304
 305         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 306         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 307         // trim(ADJUSTL(timeslice_file)) // ".")
 308
 309         status = nf90_get_var(ncid, var_id, real_var_array)
 310         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 311         // trim(ADJUSTL(timeslice_file)) // ".")
 312
 313     end subroutine read_timeslice_netcdf_real_1D_variables
 314
 315
 316     ! Read integer one dimension parameter arrays from the timeslice NetCDF
 317     subroutine read_timeslice_netcdf_integer_1D_variables(ncid, netcdf_array_name, timeslice_file, integer_var_array)
 318         integer, intent(in) :: ncid
 319         character(len=*), intent(in) :: netcdf_array_name
 320         character(len=256), intent(in) :: timeslice_file
 321         integer, dimension(:), intent(out) :: integer_var_array
 322         integer :: var_id, status
 323
 324         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 325         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 326         // trim(ADJUSTL(timeslice_file)) // ".")
 327
 328         status = nf90_get_var(ncid, var_id, integer_var_array)
 329         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in timeslice file " &
 330         // trim(ADJUSTL(timeslice_file)) // ".")
 331
 332     end subroutine read_timeslice_netcdf_integer_1D_variables
 333
 334
 335     ! Read all gage ids from reservoir parameter file
 336     subroutine read_persistence_netcdf_gage_ids_all(ncid, netcdf_array_name, reservoir_parameter_file, var_id, number_of_gages)
 337         integer, intent(in) :: ncid
 338         character(len=*), intent(in) :: netcdf_array_name
 339         character(len=*), intent(in) :: reservoir_parameter_file
 340         integer, intent(out) ::  var_id, number_of_gages
 341         integer :: status
 342         integer, dimension(nf90_max_var_dims) :: dim_ids
 343
 344         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 345         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 346         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 347
 348         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 349         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 350         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 351
 352         status = nf90_inquire_dimension(ncid, dim_ids(2), len = number_of_gages)
 353         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 354         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 355
 356     end subroutine read_persistence_netcdf_gage_ids_all
 357
 358
 359     ! Read lake id from reservoir parameter file and determine the lake_id_index, which
 360     ! is the index of the particular reservoir's values in all of the parameter arrays.
 361     subroutine read_netcdf_lake_id(ncid, lake_number, netcdf_array_name, reservoir_parameter_file, lake_id_index)
 362         integer, intent(in)          :: ncid
 363         integer(kind=int64), intent(in)  :: lake_number
 364         character(len=*), intent(in) :: netcdf_array_name
 365         character(len=*), intent(in) :: reservoir_parameter_file
 366         integer, intent(out) :: lake_id_index
 367         integer, allocatable, dimension(:) :: temp_lake_id_array
 368         integer :: var_id, status, number_of_lakes, lake_index
 369         integer, dimension(nf90_max_var_dims) :: dim_ids
 370
 371         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 372         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 373         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 374
 375         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 376         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name)
 377
 378         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_lakes)
 379         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 380         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 381
 382         allocate(temp_lake_id_array(number_of_lakes))
 383
 384         status = nf90_get_var(ncid, var_id, temp_lake_id_array)
 385         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 386         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 387
 388         do lake_index = 1, number_of_lakes
 389             if (temp_lake_id_array(lake_index) == lake_number) then
 390                 lake_id_index = lake_index
 391             end if
 392         end do
 393
 394         if(allocated(temp_lake_id_array)) deallocate(temp_lake_id_array)
 395
 396     end subroutine read_netcdf_lake_id
 397
 398
 399     ! Read gage id array from reservoir parameter file
 400     subroutine read_persistence_netcdf_gage_id(ncid, lake_id_index, netcdf_array_name, reservoir_parameter_file, gage_id)
 401         integer, intent(in) :: ncid, lake_id_index
 402         character(len=*), intent(in) :: netcdf_array_name
 403         character(len=*), intent(in) :: reservoir_parameter_file
 404         character(len=*), intent(out) :: gage_id
 405         character(len=15) :: gage_id_string, gage_id_string_trimmed
 406         character(len=15), allocatable, dimension(:) :: temp_gage_id_array
 407         integer :: var_id, status, number_of_lakes
 408         integer, dimension(nf90_max_var_dims) :: dim_ids
 409         integer        :: char_index, number_counter, temp_val
 410         character(len=1)   :: temp_char
 411
 412         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 413         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 414         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 415
 416         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 417         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 418         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 419
 420         status = nf90_inquire_dimension(ncid, dim_ids(2), len = number_of_lakes)
 421         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 422         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 423
 424         allocate(temp_gage_id_array(number_of_lakes))
 425
 426         status = nf90_get_var(ncid, var_id, temp_gage_id_array)
 427         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 428         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 429
 430         gage_id_string = temp_gage_id_array(lake_id_index)
 431
 432         do char_index = 1, 15
 433             temp_val = ichar(gage_id_string(char_index:char_index))
 434
 435             ! Check if null character, then pad with a space
 436             if (temp_val == 0) then
 437                 gage_id_string(char_index:char_index) = ' '
 438             end if
 439         end do
 440
 441         gage_id_string_trimmed = ADJUSTL(trim(gage_id_string))
 442
 443         gage_id = gage_id_string_trimmed
 444
 445         if(allocated(temp_gage_id_array)) deallocate(temp_gage_id_array)
 446
 447     end subroutine read_persistence_netcdf_gage_id
 448
 449
 450     ! Read rfc gage id array from the reservoir parameters file
 451     subroutine read_reservoir_parameters_netcdf_rfc_gage_id(ncid, lake_id_index, &
 452         netcdf_array_name, reservoir_parameter_file, rfc_gage_id)
 453         integer, intent(in) :: ncid, lake_id_index
 454         character(len=*), intent(in) :: netcdf_array_name
 455         character(len=*), intent(in) :: reservoir_parameter_file
 456         character(len=5), intent(out) :: rfc_gage_id
 457         character(len=5) :: gage_id_string, gage_id_string_trimmed
 458         character(len=5), allocatable, dimension(:) :: temp_gage_id_array
 459         integer :: var_id, status, number_of_lakes
 460         integer, dimension(nf90_max_var_dims) :: dim_ids
 461         integer        :: char_index, number_counter
 462         character(len=1)   :: temp_char
 463
 464         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 465         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 466          // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 467
 468         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 469         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name)
 470
 471         status = nf90_inquire_dimension(ncid, dim_ids(2), len = number_of_lakes)
 472         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 473         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 474
 475         allocate(temp_gage_id_array(number_of_lakes))
 476
 477         status = nf90_get_var(ncid, var_id, temp_gage_id_array)
 478         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 479         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 480
 481         gage_id_string = temp_gage_id_array(lake_id_index)
 482
 483         gage_id_string_trimmed = ADJUSTL(trim(gage_id_string))
 484
 485         rfc_gage_id = gage_id_string_trimmed
 486
 487         if(allocated(temp_gage_id_array)) deallocate(temp_gage_id_array)
 488
 489     end subroutine read_reservoir_parameters_netcdf_rfc_gage_id
 490
 491
 492     ! Read an RFC forecast integer variable for a corresponding lake id from the reservoir parameter file
 493     subroutine read_reservoir_parameters_netcdf_rfc_integer(ncid, lake_id_index, &
 494         netcdf_array_name, reservoir_parameter_file, integer_var)
 495         integer, intent(in) :: ncid, lake_id_index
 496         character(len=*), intent(in) :: netcdf_array_name
 497         character(len=*), intent(in) :: reservoir_parameter_file
 498         integer, intent(out) :: integer_var
 499         integer, allocatable, dimension(:) :: temp_integer_array
 500         integer :: var_id, status, number_of_lakes
 501         integer, dimension(nf90_max_var_dims) :: dim_ids
 502
 503         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 504         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 505         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 506
 507         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 508         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 509         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 510
 511         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_lakes)
 512         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 513         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 514
 515         allocate(temp_integer_array(number_of_lakes))
 516
 517         status = nf90_get_var(ncid, var_id, temp_integer_array)
 518         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 519         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 520
 521         integer_var = temp_integer_array(lake_id_index)
 522
 523         if(allocated(temp_integer_array)) deallocate(temp_integer_array)
 524
 525     end subroutine read_reservoir_parameters_netcdf_rfc_integer
 526
 527
 528     ! Read real two dimension parameter arrays from the reservoir parameter file
 529     subroutine read_persistence_netcdf_real_2D_parameters(ncid, lake_id_index, &
 530         netcdf_array_name, reservoir_parameter_file, var_id, number_of_weights, number_of_lakes)
 531         integer, intent(in) :: ncid, lake_id_index
 532         character(len=*), intent(in) :: netcdf_array_name
 533         character(len=*), intent(in) :: reservoir_parameter_file
 534         integer, intent(out) :: var_id, number_of_weights, number_of_lakes
 535         integer :: status
 536         integer, dimension(nf90_max_var_dims) :: dim_ids
 537
 538         status = nf90_inq_varid(ncid, netcdf_array_name, var_id)
 539         if (status /= nf90_noerr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 540         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 541
 542         status = nf90_inquire_variable(ncid, var_id, dimids = dim_ids)
 543         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 544         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 545
 546         status = nf90_inquire_dimension(ncid, dim_ids(1), len = number_of_weights)
 547         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 548         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 549
 550         status = nf90_inquire_dimension(ncid, dim_ids(2), len = number_of_lakes)
 551         if(status /= nf90_NoErr) call handle_err(status, netcdf_array_name // " read error in reservoir parameter file " &
 552         // trim(ADJUSTL(reservoir_parameter_file)) // ".")
 553
 554     end subroutine read_persistence_netcdf_real_2D_parameters
 555
 556
 557     ! Read integer parameters from the RFC Time Series NetCDF
 558     subroutine read_timeslice_netcdf_integer_variables(ncid, netcdf_variable_name, timeslice_file, integer_variable)
 559         integer, intent(in) :: ncid
 560         character(len=*), intent(in) :: netcdf_variable_name
 561         character(len=256), intent(in) :: timeslice_file
 562         integer, intent(out) :: integer_variable
 563         integer :: var_id, status
 564
 565         status = nf90_inq_varid(ncid, netcdf_variable_name, var_id)
 566         if (status /= nf90_noerr) call handle_err(status, netcdf_variable_name // " read error in timeslice file " &
 567         // trim(ADJUSTL(timeslice_file)) // ".")
 568
 569         status = nf90_get_var(ncid, var_id, integer_variable)
 570         if (status /= nf90_noerr) call handle_err(status, netcdf_variable_name // " read error in timeslice file " &
 571         // trim(ADJUSTL(timeslice_file)) // ".")
 572
 573     end subroutine read_timeslice_netcdf_integer_variables
 574
 575
 576     ! Handle NetCDF error
 577     subroutine handle_err(status, variable_error)
 578         implicit none
 579         integer, intent (in) :: status
 580         character (len=*), intent(in) :: variable_error
 581         if(status /= nf90_noerr) then
 582             call hydro_stop("ERROR: " // trim(nf90_strerror(status)) // ': ' // variable_error)
 583         end if
 584
 585     end subroutine handle_err
 586
 587
 588     ! Create Persistence Levelpool Hybrid Reservoir Function Diagnostic Log CSV File
 589     subroutine create_hybrid_diagnostic_log_file(lake_number)
 590         integer, intent(in) :: lake_number
 591         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 592         character(len=50)   :: filename_string
 593
 594         write(lake_number_string, "(I15)") lake_number
 595         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 596         filename_string = "hybrid_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 597
 598         open (113, file=filename_string, status="unknown")
 599
 600         write (113, "(A256)") "Current Time, Timeslice Update Time, Weight Update Time, &
 601         Gage Lookback Seconds, Gage ID, Current Persistence Weight Index, &
 602         Current Persistence Weight, Inflow, Storage, Water Elevation, Gage Discharge, &
 603         Persisted Outflow, Levelpool Outflow, Weighted Outflow"
 604
 605         close(113)
 606
 607     end subroutine create_hybrid_diagnostic_log_file
 608
 609
 610     ! Log data to Persistence Levelpool Hybrid Reservoir Function Diagnostic Log CSV File
 611     subroutine log_hybrid_diagnostic_data(lake_number, current_time, timeslice_update_time, &
 612         weight_update_time, gage_lookback_seconds, gage_id, persistence_weight_index, &
 613         persistence_current_weight, inflow, current_storage, water_elevation, gage_discharge, &
 614         persisted_outflow, levelpool_outflow, outflow)
 615         integer, intent(in) :: lake_number, current_time, timeslice_update_time
 616         integer, intent(in) :: weight_update_time, gage_lookback_seconds
 617         integer, intent(in) :: persistence_weight_index
 618         real,    intent(in) :: persistence_current_weight, inflow, current_storage
 619         real,    intent(in) :: water_elevation, gage_discharge, persisted_outflow
 620         real,    intent(in) :: levelpool_outflow, outflow
 621         character(len=*), intent(in) :: gage_id
 622         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 623         character(len=50)   :: filename_string
 624
 625         write(lake_number_string, "(I15)") lake_number
 626         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 627         filename_string = "hybrid_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 628
 629         open (113, file=filename_string, status="unknown", position = 'append')
 630
 631         write (113, "(I20, A1, I20, A1, I20, A1, I20, A1, A15, A1, I20, A1, F15.5, A1, F15.5, &
 632         A1, F15.15, A1, F15.5, A1, F15.5, A1, F15.5, A1, F15.5, A1, F15.5)") &
 633         current_time, ',', timeslice_update_time, ',', weight_update_time, ',', &
 634         gage_lookback_seconds, ',', gage_id, ',', persistence_weight_index, ',', &
 635         persistence_current_weight, ',', inflow, ',', current_storage, ',', &
 636         water_elevation, ',', gage_discharge, ',', persisted_outflow, ',', levelpool_outflow, &
 637         ',', outflow
 638
 639         close (113)
 640
 641     end subroutine log_hybrid_diagnostic_data
 642
 643
 644     ! Create RFC Forecasts Reservoir Function Diagnostic Log CSV File
 645     subroutine create_rfc_forecasts_diagnostic_log_file(lake_number)
 646         integer, intent(in) :: lake_number
 647         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 648         character(len=50)   :: filename_string
 649
 650         write(lake_number_string, "(I15)") lake_number
 651         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 652         filename_string = "rfc_forecasts_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 653
 654         open (113, file=filename_string, status="unknown")
 655
 656         write (113, "(A256)") "Current_Time, Time_Step, Time_Series_Update_Time, &
 657         Lookback_Seconds, Time_Series_Index, Gage_ID, Inflow, Water_Elevation, &
 658         Levelpool_Outflow, Outflow"
 659
 660         close(113)
 661
 662     end subroutine create_rfc_forecasts_diagnostic_log_file
 663
 664
 665     ! Log data to RFC Forecasts Reservoir Function Diagnostic Log CSV File
 666     subroutine log_rfc_forecasts_diagnostic_data(lake_number, current_time, time_step, &
 667         time_series_update_time, lookback_seconds, time_series_index, gage_id,  &
 668         inflow, water_elevation, levelpool_outflow, outflow)
 669         integer, intent(in) :: lake_number, current_time, time_step, time_series_update_time
 670         integer, intent(in) :: lookback_seconds, time_series_index
 671         character(len=*), intent(in) :: gage_id
 672         real,    intent(in) :: inflow, water_elevation, levelpool_outflow, outflow
 673         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 674         character(len=50)   :: filename_string
 675
 676         write(lake_number_string, "(I15)") lake_number
 677         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 678         filename_string = "rfc_forecasts_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 679
 680         open (113, file=filename_string, status="unknown", position = 'append')
 681
 682         write (113, "(I20, A1, I20, A1, I20, A1, I20, A1, I20, A1, A5, A1, F15.5, A1, &
 683         F15.5, A1, F15.5, A1, F15.5)") &
 684         current_time, ',', time_step, ',', time_series_update_time, ',', &
 685         lookback_seconds, ',', time_series_index, ',', gage_id, ',', &
 686         inflow, ',', water_elevation, ',', levelpool_outflow, ',', outflow
 687
 688         close (113)
 689
 690     end subroutine log_rfc_forecasts_diagnostic_data
 691
 692
 693     ! Create Levelpool Reservoir Function Diagnostic Log CSV File
 694     subroutine create_levelpool_diagnostic_log_file(lake_number)
 695         integer, intent(in) :: lake_number
 696         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 697         character(len=50)   :: filename_string
 698
 699         write(lake_number_string, "(I15)") lake_number
 700         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 701
 702         filename_string = "levelpool_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 703
 704         open (113, file=filename_string, status="unknown")
 705
 706         write (113, "(A256)") "Inflow, Water Elevation, Outflow"
 707
 708         close(113)
 709
 710     end subroutine create_levelpool_diagnostic_log_file
 711
 712
 713     ! Log data to Levelpool Reservoir Function Diagnostic Log CSV File
 714     subroutine log_levelpool_diagnostic_data(lake_number, inflow, water_elevation, outflow)
 715         integer, intent(in) :: lake_number
 716         real,    intent(in) :: inflow
 717         real,    intent(in) :: water_elevation
 718         real,    intent(in) :: outflow
 719         character(len=15)   :: lake_number_string, lake_number_string_trimmed
 720         character(len=50)   :: filename_string
 721
 722         write(lake_number_string, "(I15)") lake_number
 723         lake_number_string_trimmed = ADJUSTL(trim(lake_number_string))
 724         filename_string = "levelpool_logs_"//ADJUSTL(trim(lake_number_string_trimmed))//".csv"
 725
 726         open (113, file=filename_string, status="unknown", position = 'append')
 727
 728         write (113, "(F15.8, A1, F15.8, A1, F15.8)") &
 729         inflow, ',', water_elevation, ',', outflow
 730
 731         close (113)
 732
 733     end subroutine log_levelpool_diagnostic_data
 734
 735
 736     ! Get a new date given an existing date and time interval
 737     subroutine geth_newdate (ndate, odate, idt)
 738     implicit none
 739
 740     !  From old date ("YYYY-MM-DD HH:MM:SS.ffff" or "YYYYMMDDHHMMSSffff") and
 741     !  delta-time, compute the new date.
 742
 743     !  on entry     -  odate  -  the old hdate.
 744     !                  idt    -  the change in time
 745
 746     !  on exit      -  ndate  -  the new hdate.
 747
 748     integer, intent(in)           :: idt
 749     character (len=*), intent(out) :: ndate
 750     character (len=*), intent(in)  :: odate
 751
 752     !  Local Variables
 753
 754     !  yrold    -  indicates the year associated with "odate"
 755     !  moold    -  indicates the month associated with "odate"
 756     !  dyold    -  indicates the day associated with "odate"
 757     !  hrold    -  indicates the hour associated with "odate"
 758     !  miold    -  indicates the minute associated with "odate"
 759     !  scold    -  indicates the second associated with "odate"
 760
 761     !  yrnew    -  indicates the year associated with "ndate"
 762     !  monew    -  indicates the month associated with "ndate"
 763     !  dynew    -  indicates the day associated with "ndate"
 764     !  hrnew    -  indicates the hour associated with "ndate"
 765     !  minew    -  indicates the minute associated with "ndate"
 766     !  scnew    -  indicates the second associated with "ndate"
 767
 768     !  mday     -  a list assigning the number of days in each month
 769
 770     !  i        -  loop counter
 771     !  nday     -  the integer number of days represented by "idt"
 772     !  nhour    -  the integer number of hours in "idt" after taking out
 773     !              all the whole days
 774     !  nmin     -  the integer number of minutes in "idt" after taking out
 775     !              all the whole days and whole hours.
 776     !  nsec     -  the integer number of minutes in "idt" after taking out
 777     !              all the whole days, whole hours, and whole minutes.
 778
 779     integer :: newlen, oldlen
 780     integer :: yrnew, monew, dynew, hrnew, minew, scnew, frnew
 781     integer :: yrold, moold, dyold, hrold, miold, scold, frold
 782     integer :: nday, nhour, nmin, nsec, nfrac, i, ifrc
 783     logical :: opass
 784     character (len=10) :: hfrc
 785     character (len=1) :: sp
 786     logical :: punct
 787     integer :: yrstart, yrend, mostart, moend, dystart, dyend
 788     integer :: hrstart, hrend, mistart, miend, scstart, scend, frstart
 789     integer :: units
 790     integer, dimension(12) :: mday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
 791
 792     ! Determine if odate is "YYYY-MM-DD_HH ... " or "YYYYMMDDHH...."
 793     if (odate(5:5) == "-") then
 794        punct = .TRUE.
 795     else
 796        punct = .FALSE.
 797     endif
 798
 799     !  Break down old hdate into parts
 800
 801     hrold = 0
 802     miold = 0
 803     scold = 0
 804     frold = 0
 805     oldlen = LEN(odate)
 806     if (punct) then
 807        yrstart = 1
 808        yrend = 4
 809        mostart = 6
 810        moend = 7
 811        dystart = 9
 812        dyend = 10
 813        hrstart = 12
 814        hrend = 13
 815        mistart = 15
 816        miend = 16
 817        scstart = 18
 818        scend = 19
 819        frstart = 21
 820        select case (oldlen)
 821        case (10)
 822           ! Days
 823           units = 1
 824        case (13)
 825           ! Hours
 826           units = 2
 827        case (16)
 828           ! Minutes
 829           units = 3
 830        case (19)
 831           ! Seconds
 832           units = 4
 833        case (21)
 834           ! Tenths
 835           units = 5
 836        case (22)
 837           ! Hundredths
 838           units = 6
 839        case (23)
 840           ! Thousandths
 841           units = 7
 842        case (24)
 843           ! Ten thousandths
 844           units = 8
 845        case default
 846
 847           write(*,*) 'FATAL ERROR: geth_newdate:  odd length: #'//trim(odate)//'#'
 848           call hydro_stop("geth_newdate")
 849
 850        end select
 851
 852        if (oldlen.ge.11) then
 853           sp = odate(11:11)
 854        else
 855           sp = ' '
 856        end if
 857
 858     else
 859
 860        yrstart = 1
 861        yrend = 4
 862        mostart = 5
 863        moend = 6
 864        dystart = 7
 865        dyend = 8
 866        hrstart = 9
 867        hrend = 10
 868        mistart = 11
 869        miend = 12
 870        scstart = 13
 871        scend = 14
 872        frstart = 15
 873
 874        select case (oldlen)
 875        case (8)
 876           ! Days
 877           units = 1
 878        case (10)
 879           ! Hours
 880           units = 2
 881        case (12)
 882           ! Minutes
 883           units = 3
 884        case (14)
 885           ! Seconds
 886           units = 4
 887        case (15)
 888           ! Tenths
 889           units = 5
 890        case (16)
 891           ! Hundredths
 892           units = 6
 893        case (17)
 894           ! Thousandths
 895           units = 7
 896        case (18)
 897           ! Ten thousandths
 898           units = 8
 899        case default
 900
 901           write(*,*) 'FATAL ERROR: geth_newdate:  odd length: #'//trim(odate)//'#'
 902           call hydro_stop("geth_newdate")
 903
 904        end select
 905     endif
 906
 907     !  Use internal READ statements to convert the CHARACTER string
 908     !  date into INTEGER components.
 909
 910     read(odate(yrstart:yrend),  '(i4)') yrold
 911     read(odate(mostart:moend),  '(i2)') moold
 912     read(odate(dystart:dyend), '(i2)') dyold
 913     if (units.ge.2) then
 914        read(odate(hrstart:hrend),'(i2)') hrold
 915        if (units.ge.3) then
 916           read(odate(mistart:miend),'(i2)') miold
 917           if (units.ge.4) then
 918              read(odate(scstart:scend),'(i2)') scold
 919              if (units.ge.5) then
 920                 read(odate(frstart:oldlen),*) frold
 921              end if
 922           end if
 923        end if
 924     end if
 925
 926     !  Set the number of days in February for that year.
 927
 928     mday(2) = nfeb(yrold)
 929
 930     !  Check that ODATE makes sense.
 931
 932     opass = .TRUE.
 933
 934     !  Check that the month of ODATE makes sense.
 935
 936     if ((moold.gt.12).or.(moold.lt.1)) then
 937        write(*,*) 'GETH_NEWDATE:  Month of ODATE = ', moold
 938        opass = .FALSE.
 939     end if
 940
 941     !  Check that the day of ODATE makes sense.
 942
 943     if ((dyold.gt.mday(moold)).or.(dyold.lt.1)) then
 944        write(*,*) 'GETH_NEWDATE:  Day of ODATE = ', dyold
 945        opass = .FALSE.
 946     end if
 947
 948     !  Check that the hour of ODATE makes sense.
 949
 950     if ((hrold.gt.23).or.(hrold.lt.0)) then
 951        write(*,*) 'GETH_NEWDATE:  Hour of ODATE = ', hrold
 952        opass = .FALSE.
 953     end if
 954
 955     !  Check that the minute of ODATE makes sense.
 956
 957     if ((miold.gt.59).or.(miold.lt.0)) then
 958        write(*,*) 'GETH_NEWDATE:  Minute of ODATE = ', miold
 959        opass = .FALSE.
 960     end if
 961
 962     !  Check that the second of ODATE makes sense.
 963
 964     if ((scold.gt.59).or.(scold.lt.0)) then
 965        write(*,*) 'GETH_NEWDATE:  Second of ODATE = ', scold
 966        opass = .FALSE.
 967     end if
 968
 969     !  Check that the fractional part  of ODATE makes sense.
 970
 971
 972     if (.not.opass) then
 973
 974        write(*,*) 'FATAL ERROR: Crazy ODATE: ', odate(1:oldlen), oldlen
 975        stop
 976
 977     end if
 978
 979     !  Date Checks are completed.  Continue.
 980
 981
 982     !  Compute the number of days, hours, minutes, and seconds in idt
 983
 984     if (units.ge.5) then !idt should be in fractions of seconds
 985        ifrc = oldlen-(frstart)+1
 986        ifrc = 10**ifrc
 987        nday   = abs(idt)/(86400*ifrc)
 988        nhour  = mod(abs(idt),86400*ifrc)/(3600*ifrc)
 989        nmin   = mod(abs(idt),3600*ifrc)/(60*ifrc)
 990        nsec   = mod(abs(idt),60*ifrc)/(ifrc)
 991        nfrac = mod(abs(idt), ifrc)
 992     else if (units.eq.4) then  !idt should be in seconds
 993        ifrc = 1
 994        nday   = abs(idt)/86400 ! integer number of days in delta-time
 995        nhour  = mod(abs(idt),86400)/3600
 996        nmin   = mod(abs(idt),3600)/60
 997        nsec   = mod(abs(idt),60)
 998        nfrac  = 0
 999     else if (units.eq.3) then !idt should be in minutes
1000        ifrc = 1
1001        nday   = abs(idt)/1440 ! integer number of days in delta-time
1002        nhour  = mod(abs(idt),1440)/60
1003        nmin   = mod(abs(idt),60)
1004        nsec   = 0
1005        nfrac  = 0
1006     else if (units.eq.2) then !idt should be in hours
1007        ifrc = 1
1008        nday   = abs(idt)/24 ! integer number of days in delta-time
1009        nhour  = mod(abs(idt),24)
1010        nmin   = 0
1011        nsec   = 0
1012        nfrac  = 0
1013     else if (units.eq.1) then !idt should be in days
1014        ifrc = 1
1015        nday   = abs(idt)    ! integer number of days in delta-time
1016        nhour  = 0
1017        nmin   = 0
1018        nsec   = 0
1019        nfrac  = 0
1020     else
1021
1022        write(*,'(''GETH_NEWDATE: Strange length for ODATE: '', i3)') &
1023             oldlen
1024        write(*,*) '#'//odate(1:oldlen)//'#'
1025        call hydro_stop("geth_newdate")
1026
1027     end if
1028
1029     if (idt.ge.0) then
1030
1031        frnew = frold + nfrac
1032        if (frnew.ge.ifrc) then
1033           frnew = frnew - ifrc
1034           nsec = nsec + 1
1035        end if
1036
1037        scnew = scold + nsec
1038        if (scnew .ge. 60) then
1039           scnew = scnew - 60
1040           nmin  = nmin + 1
1041        end if
1042
1043        minew = miold + nmin
1044        if (minew .ge. 60) then
1045           minew = minew - 60
1046           nhour  = nhour + 1
1047        end if
1048
1049        hrnew = hrold + nhour
1050        if (hrnew .ge. 24) then
1051           hrnew = hrnew - 24
1052           nday  = nday + 1
1053        end if
1054
1055        dynew = dyold
1056        monew = moold
1057        yrnew = yrold
1058        do i = 1, nday
1059           dynew = dynew + 1
1060           if (dynew.gt.mday(monew)) then
1061              dynew = dynew - mday(monew)
1062              monew = monew + 1
1063              if (monew .gt. 12) then
1064                 monew = 1
1065                 yrnew = yrnew + 1
1066                 ! If the year changes, recompute the number of days in February
1067                 mday(2) = nfeb(yrnew)
1068              end if
1069           end if
1070        end do
1071
1072     else if (idt.lt.0) then
1073
1074        frnew = frold - nfrac
1075        if (frnew .lt. 0) then
1076           frnew = frnew + ifrc
1077           nsec = nsec + 1
1078        end if
1079
1080        scnew = scold - nsec
1081        if (scnew .lt. 00) then
1082           scnew = scnew + 60
1083           nmin  = nmin + 1
1084        end if
1085
1086        minew = miold - nmin
1087        if (minew .lt. 00) then
1088           minew = minew + 60
1089           nhour  = nhour + 1
1090        end if
1091
1092        hrnew = hrold - nhour
1093        if (hrnew .lt. 00) then
1094           hrnew = hrnew + 24
1095           nday  = nday + 1
1096        end if
1097
1098        dynew = dyold
1099        monew = moold
1100        yrnew = yrold
1101        do i = 1, nday
1102           dynew = dynew - 1
1103           if (dynew.eq.0) then
1104              monew = monew - 1
1105              if (monew.eq.0) then
1106                 monew = 12
1107                 yrnew = yrnew - 1
1108                 ! If the year changes, recompute the number of days in February
1109                 mday(2) = nfeb(yrnew)
1110              end if
1111              dynew = mday(monew)
1112           end if
1113        end do
1114     end if
1115
1116     !  Now construct the new mdate
1117
1118     newlen = LEN(ndate)
1119
1120     if (punct) then
1121
1122        if (newlen.gt.frstart) then
1123           write(ndate(1:scend),19) yrnew, monew, dynew, hrnew, minew, scnew
1124           write(hfrc,'(i10)') frnew+1000000000
1125           ndate = ndate(1:scend)//'.'//hfrc(31-newlen:10)
1126
1127        else if (newlen.eq.scend) then
1128           write(ndate(1:scend),19) yrnew, monew, dynew, hrnew, minew, scnew
1129     19    format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2,':',i2.2)
1130
1131        else if (newlen.eq.miend) then
1132           write(ndate,16) yrnew, monew, dynew, hrnew, minew
1133     16    format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2)
1134
1135        else if (newlen.eq.hrend) then
1136           write(ndate,13) yrnew, monew, dynew, hrnew
1137     13    format(i4,'-',i2.2,'-',i2.2,'_',i2.2)
1138
1139        else if (newlen.eq.dyend) then
1140           write(ndate,10) yrnew, monew, dynew
1141     10    format(i4,'-',i2.2,'-',i2.2)
1142
1143        end if
1144
1145     else
1146
1147        if (newlen.gt.frstart) then
1148           write(ndate(1:scend),119) yrnew, monew, dynew, hrnew, minew, scnew
1149           write(hfrc,'(i10)') frnew+1000000000
1150           ndate = ndate(1:scend)//'.'//hfrc(31-newlen:10)
1151
1152        else if (newlen.eq.scend) then
1153           write(ndate(1:scend),119) yrnew, monew, dynew, hrnew, minew, scnew
1154     119   format(i4,i2.2,i2.2,i2.2,i2.2,i2.2)
1155
1156        else if (newlen.eq.miend) then
1157           write(ndate,116) yrnew, monew, dynew, hrnew, minew
1158     116   format(i4,i2.2,i2.2,i2.2,i2.2)
1159
1160        else if (newlen.eq.hrend) then
1161           write(ndate,113) yrnew, monew, dynew, hrnew
1162     113   format(i4,i2.2,i2.2,i2.2)
1163
1164        else if (newlen.eq.dyend) then
1165           write(ndate,110) yrnew, monew, dynew
1166     110   format(i4,i2.2,i2.2)
1167
1168        end if
1169
1170     endif
1171
1172     if (punct .and. (oldlen.ge.11) .and. (newlen.ge.11)) ndate(11:11) = sp
1173
1174     end subroutine geth_newdate
1175
1176     integer function nfeb(year)
1177     !
1178     ! Compute the number of days in February for the given year.
1179     !
1180     implicit none
1181     integer, intent(in) :: year ! Four-digit year
1182
1183     nfeb = 28 ! By default, February has 28 days ...
1184     if (mod(year,4).eq.0) then
1185        nfeb = 29  ! But every four years, it has 29 days ...
1186        if (mod(year,100).eq.0) then
1187           nfeb = 28  ! Except every 100 years, when it has 28 days ...
1188           if (mod(year,400).eq.0) then
1189              nfeb = 29  ! Except every 400 years, when it has 29 days ...
1190              if (mod(year,3600).eq.0) then
1191                 nfeb = 28  ! Except every 3600 years, when it has 28 days.
1192              endif
1193           endif
1194        endif
1195     endif
1196     end function nfeb
1197
1198 end module module_reservoir_utilities