Skip to content

Commit

Permalink
fix(rte): Fixed issues pertaining to fhr.
Browse files Browse the repository at this point in the history
One of the changes by MESH-Model/MESH_Code@52c7367 "LongSimFix" was to revise the use of `fhr` from tracking incremental hours in the simulation to always being equal to `1`.

When originally implemented, a temporary fix for observed variables arbitrarily allocated these arrays to `999999` blocks (observation hours). This caused issues for simulations exceeding 999999 hours (approx.. 114 years); for example, climate simulations from 1980 to 2100, which likely would cause an index out-of-bounds error in 2094 or later.

Changes by @mee067 to address this removed the allocation of these variables from `999999` to `(:, 1)`, setting the current record `fhr` always equal to `1`. However, the change only initialized `fhr = 1` when resumed from `seq` format resume files. `fhr` not being properly initialized without this condition seems to have contributed to the corruption of data used by the module, i.e., the record of observations index by gauge ID, causing the indexing issue observed with the streamflow observations printed by `save_basin_output`.

Adding the same condition, `fhr = 1`, during initialization resolves this issue.
```f90
!> Set fhr to 1 (the counting of incremental hours in the routing simulation is not used like it is in standalone WATROUTE).
fhr = 1
```

Setting `fhr = 1` was correctly implemented in ab89db4. However, further modifications have been carried over to allocate the affected variables to `fhr` than to `1`, so the connection to this dimension is clearer. Changes have also been made when resuming from the `seq` format resume files for the case when the file was written by a previous version of the code where `fhr /= 1`.
```f90
real(kind = 4), dimension(:, :), allocatable :: lake_elv_temp
...

if (fms%rsvr%n > 0) then
    allocate(lake_elv_temp(noresv, fhr))
    read(iun) lake_elv_temp(:, fhr)
    lake_elv(:, 1) = lake_elv_temp(:, fhr)
else
    read(iun)
end if
...

!> Set fhr to 1 (the counting of incremental hours in the routing simulation is not used like it is in standalone WATROUTE).
fhr = 1
```
  • Loading branch information
dprincz committed Aug 13, 2024
1 parent ab89db4 commit 5fc2d3b
Showing 1 changed file with 15 additions and 7 deletions.
22 changes: 15 additions & 7 deletions Routing_Model/RPN_watroute/sa_mesh_process/rte_module.f90
Original file line number Diff line number Diff line change
Expand Up @@ -332,12 +332,14 @@ subroutine run_rte_init(fls, shd)
!todo: water class for lakes.
ii_water = ntype

!> Set fhr to 1 (the counting of incremental hours in the routing simulation is not used like it is in standalone WATROUTE).
fhr = 1

!> Reservoirs.
!* noresv: Number of reservoirs/lakes.
!* Nreaches: (rerout.f) Copy of noresv.
noresv = fms%rsvr%n
Nreaches = fms%rsvr%n
fhr = 1
if (fms%rsvr%n > 0) then

!> Book-keeping variables (when fhr > 1).
Expand All @@ -346,7 +348,7 @@ subroutine run_rte_init(fls, shd)
!* lake_stor: Copy of reservoir storage. [m3].
!* lake_outflow: Copy of reservoir outflow. [m3 s-1].
!* del_store: Storage change considering inflow minus outflow. [m3].
allocate(lake_inflow(noresv, 1), lake_stor(noresv, 1), lake_outflow(noresv, 1), del_stor(noresv, 1))
allocate(lake_inflow(noresv, fhr), lake_stor(noresv, fhr), lake_outflow(noresv, fhr), del_stor(noresv, fhr))
lake_inflow = 0.0; lake_stor = 0.0; lake_outflow = 0.0; del_stor = 0.0

!> Reservoir/lake meta-data (from file).
Expand Down Expand Up @@ -382,7 +384,7 @@ subroutine run_rte_init(fls, shd)
!> Measured outflow (from file).
!* qrel: Measured outflow when reservoir releases are replaced with measured values. [m3 s-1].
!* qdwpr: (Local variable in route.f) Used to accumulate flow in reaches that span multiple cells to the reservoir outlet. [m3 s-1].
allocate(qdwpr(noresv, 1), qrel(noresv, 1))
allocate(qdwpr(noresv, fhr), qrel(noresv, fhr))
qdwpr = 0.0
if (count(fms%rsvr%rls%b1 == 0.0) > 0 .and. fms%rsvr%rlsmeas%readmode /= 'n') then
qrel(1:count(fms%rsvr%rls%b1 == 0.0), 1) = real(fms%rsvr%rlsmeas%val(1:count(fms%rsvr%rls%b1 == 0.0)), kind(qrel))
Expand All @@ -394,7 +396,7 @@ subroutine run_rte_init(fls, shd)
!* reach_last: (Used in Great Lakes routing) Stores level of last time-step. [m].
!* lake_area: Used to convert reservoir storage to level (for diagnostic output). [m2].
!* lake_elv: Lake elevation (for diagnostic output). [m].
allocate(reach_last(noresv), lake_area(noresv), lake_elv(noresv, 1))
allocate(reach_last(noresv), lake_area(noresv), lake_elv(noresv, fhr))
lake_area = real(fms%rsvr%rls%area, kind(lake_area))
lake_elv(:, 1) = real(fms%rsvr%rls%zlvl0, kind(lake_elv))
end if
Expand All @@ -404,7 +406,7 @@ subroutine run_rte_init(fls, shd)
if (fms%stmg%n > 0) then
allocate( &
iflowgrid(no), nopt(no), &
qhyd(no, 1))
qhyd(no, fhr))
iflowgrid = fms%stmg%meta%rnk
nopt = -1
qhyd(:, 1) = real(fms%stmg%qomeas%val, kind(qhyd))
Expand Down Expand Up @@ -581,6 +583,7 @@ subroutine run_rte_resume_read(fls, shd)
!> Local variables.
integer(kind = 4) fhr_i4
integer ierr, iun
real(kind = 4), dimension(:, :), allocatable :: lake_elv_temp

!> Return if not the head node or if the process is not active.
if (.not. ISHEADNODE .or. .not. rteflg%PROCESS_ACTIVE) return
Expand All @@ -593,19 +596,24 @@ subroutine run_rte_resume_read(fls, shd)

!> Read inital values from the file.
read(iun) fhr_i4
fhr = 1
fhr = int(fhr_i4)
read(iun) qo2
read(iun) store2
read(iun) qi2
if (fms%rsvr%n > 0) then
read(iun) lake_elv(:, fhr)
allocate(lake_elv_temp(noresv, fhr))
read(iun) lake_elv_temp(:, fhr)
lake_elv(:, 1) = lake_elv_temp(:, fhr)
else
read(iun)
end if

!> Close the file to free the unit.
close(iun)

!> Set fhr to 1 (the counting of incremental hours in the routing simulation is not used like it is in standalone WATROUTE).
fhr = 1

!> Update SA_MESH output variables.
if (associated(out%ts%grid%qi)) out%ts%grid%qi = qi2
if (associated(out%ts%grid%stgch)) out%ts%grid%stgch = store2
Expand Down

0 comments on commit 5fc2d3b

Please sign in to comment.