Skip to content

Commit

Permalink
Merge pull request #357 from nmizukami/develop
Browse files Browse the repository at this point in the history 
Modify direct-insertion and other minor changes
  • Loading branch information
@nmizukami
nmizukami authored Mar 30, 2023
2 parents eb60504 + e15303c commit 6a96294
Show file tree
Hide file tree
Showing 14 changed files with 232 additions and 246 deletions.
34 changes: 17 additions & 17 deletions docs/source/Input_data.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,23 +32,23 @@ Dimensions required

Minimum variables required

+------------+------------+-----------+-------+-----------------------------------------+
| Variable | Dimension | Unit | Type | Description |
+============+============+===========+=======+=========================================+
| segId | seg | ``-`` | int | unique id of each stream segment |
+------------+------------+-----------+-------+-----------------------------------------+
| HRUid | hru | ``-`` | int | unique hru ID |
+------------+------------+-----------+-------+-----------------------------------------+
| downSegId | seg | ``-`` | int | id of the downstream segment |
+------------+------------+-----------+-------+-----------------------------------------+
| hruSegId | hru | ``-`` | int | id of the stream segment below each HRU |
+------------+------------+-----------+-------+-----------------------------------------+
| area | hru | m2 | real | hru area |
+------------+------------+-----------+-------+-----------------------------------------+
| slope | seg | ``-`` | real | slope of segment |
+------------+------------+-----------+-------+-----------------------------------------+
| length | seg | m | real | length of segment |
+------------+------------+-----------+-------+-----------------------------------------+
+------------+------------+-----------+-------+--------------------------------------------+
| Variable | Dimension | Unit | Type | Description |
+============+============+===========+=======+============================================+
| segId | seg | ``-`` | int | unique id of each stream segment |
+------------+------------+-----------+-------+--------------------------------------------+
| HRUid | hru | ``-`` | int | unique hru ID |
+------------+------------+-----------+-------+--------------------------------------------+
| downSegId | seg | ``-`` | int | id of the downstream segment |
+------------+------------+-----------+-------+--------------------------------------------+
| hruSegId | hru | ``-`` | int | id of the stream segment the HRU flows into|
+------------+------------+-----------+-------+--------------------------------------------+
| area | hru | m2 | real | hru area |
+------------+------------+-----------+-------+--------------------------------------------+
| slope | seg | ``-`` | real | slope of segment |
+------------+------------+-----------+-------+--------------------------------------------+
| length | seg | m | real | length of segment |
+------------+------------+-----------+-------+--------------------------------------------+

Negative or zero (<=0) value for downSegId is reserved for no downstream reach, meaning that such reach or hru does not flow into any reach.
(i.e., basin outlet). For this reason, segID is required to use positive integer value (>0).
Expand Down
1 change: 1 addition & 0 deletions route/build/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,7 @@ IO = \
# CORE
CORE = \
model_finalize.f90 \
data_assimilation.f90 \
accum_runoff.f90 \
basinUH.f90 \
irf_route.f90 \
Expand Down
98 changes: 98 additions & 0 deletions route/build/src/data_assimilation.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
MODULE data_assimilation

! general descriptions
!
!

USE nrtype
! data types
USE dataTypes, ONLY: STRFLX ! fluxes in each reach
USE dataTypes, ONLY: STRSTA ! state in each reach
USE dataTypes, ONLY: RCHTOPO ! Network topology
! global data
USE public_var, ONLY: iulog ! i/o logical unit number
USE public_var, ONLY: qBlendPeriod ! number of time steps for which direct insertion is performed
USE public_var, ONLY: QerrTrend ! temporal discharge error trend: 1->constant,2->linear, 3->logistic

implicit none

private
public::direct_insertion

CONTAINS

! *********************************************************************
! subroutine: perform direct-insertion
! *********************************************************************
SUBROUTINE direct_insertion(iEns, segIndex, & ! input: index of runoff ensemble to be processed
idxRoute, & ! input: reachID to be checked by on-screen pringing
ixDesire, & ! input: reachID to be checked by on-screen pringing
NETOPO_in, & ! input: reach topology data structure
RCHSTA_out, & ! inout: reach state data structure
RCHFLX_out, & ! inout: reach flux data structure
ierr, message) ! output: error control
implicit none
! Argument variables
integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed
integer(i4b), intent(in) :: segIndex ! segment where routing is performed
integer(i4b), intent(in) :: idxRoute ! index of routing method
integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output
type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology
type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data
type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains
integer(i4b), intent(out) :: ierr ! error code
character(*), intent(out) :: message ! error message
! Local variables
logical(lgt) :: verbose ! check details of variables
real(dp) :: Qcorrect ! Discharge correction (when qmodOption=1) [m3/s]
real(dp) :: k ! the logistic decay rate or steepness of the curve.
real(dp) :: y0 !
real(dp) :: x0 !
character(len=strLen) :: cmessage ! error message from subroutine
integer(i4b),parameter :: const=1 ! error reduction type: step function
integer(i4b),parameter :: linear=2 ! error reduction type: linear function
integer(i4b),parameter :: logistic=3 ! error reduction type: logistic function
integer(i4b),parameter :: exponential=4 ! error reduction type: exponential function

ierr=0; message='direct_insertion/'

verbose = .false.
if(NETOPO_in(segIndex)%REACHIX == ixDesire)then
verbose = .true.
end if

if (RCHFLX_out(iens,segIndex)%QOBS>0._dp) then ! there is observation
RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror = RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%REACH_Q - RCHFLX_out(iens,segIndex)%QOBS ! compute error
end if

if (RCHFLX_out(iens,segIndex)%Qelapsed > qBlendPeriod) then
RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror=0._dp
end if

if (RCHFLX_out(iens,segIndex)%Qelapsed <= qBlendPeriod) then
select case(QerrTrend)
case(const)
Qcorrect = RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror
case(linear)
Qcorrect = RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror*(1._dp - real(RCHFLX_out(iens,segIndex)%Qelapsed,dp)/real(qBlendPeriod, dp))
case(logistic)
x0 =0.25; y0 =0.90
k = log(1._dp/y0-1._dp)/(qBlendPeriod/2._dp-qBlendPeriod*x0)
Qcorrect = RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror/(1._dp + exp(-k*(1._dp*RCHFLX_out(iens,segIndex)%Qelapsed-qBlendPeriod/2._dp)))
case(exponential)
if (RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror/=0._dp) then
k = log(0.1_dp/abs(RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror))/(1._dp*qBlendPeriod)
Qcorrect = RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%Qerror*exp(k*RCHFLX_out(iens,segIndex)%Qelapsed)
else
Qcorrect = 0._dp
end if
case default; message=trim(message)//'discharge error trend model must be 1(const),2(liear), or 3(logistic)'; ierr=81; return
end select
else
Qcorrect=0._dp
end if
RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%REACH_Q = max(RCHFLX_out(iens,segIndex)%ROUTE(idxRoute)%REACH_Q-Qcorrect, 0._dp)

END SUBROUTINE direct_insertion

END MODULE data_assimilation
59 changes: 15 additions & 44 deletions route/build/src/dfw_route.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,10 @@ MODULE dfw_route_module
USE public_var, ONLY: realMissing ! missing value for real number
USE public_var, ONLY: integerMissing ! missing value for integer number
USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion)
USE public_var, ONLY: ntsQmodStop ! number of time steps for which direct insertion is performed
USE public_var, ONLY: QerrTrend ! temporal discharge error trend: 1->constant,2->linear, 3->logistic
USE globalData, ONLY: idxDW
! subroutines: general
USE model_finalize, ONLY : handle_err
USE data_assimilation, ONLY: direct_insertion ! qmod option (use 1==direct insertion)
USE model_finalize, ONLY: handle_err

implicit none

Expand Down Expand Up @@ -156,15 +155,7 @@ SUBROUTINE dfw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be pro
integer(i4b) :: iUps ! upstream reach index
integer(i4b) :: iRch_ups ! index of upstream reach in NETOPO
real(dp) :: q_upstream ! total discharge at top of the reach being processed
real(dp) :: Qcorrect ! Discharge correction (when qmodOption=1) [m3/s]
real(dp) :: k ! the logistic decay rate or steepness of the curve.
real(dp) :: y0 !
real(dp) :: x0 !
character(len=strLen) :: cmessage ! error message from subroutine
integer(i4b),parameter :: const=1 ! error reduction type: step function
integer(i4b),parameter :: linear=2 ! error reduction type: linear function
integer(i4b),parameter :: logistic=3 ! error reduction type: logistic function
integer(i4b),parameter :: exponential=4 ! error reduction type: exponential function

ierr=0; message='dfw_rch/'

Expand All @@ -182,39 +173,6 @@ SUBROUTINE dfw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be pro
do iUps = 1,nUps
if (.not. NETOPO_in(segIndex)%goodBas(iUps)) cycle
iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach

if (qmodOption==1) then
if (RCHFLX_out(iens,iRch_ups)%QOBS>0._dp) then
RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%Qerror = RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q - RCHFLX_out(iens,iRch_ups)%QOBS ! compute error
end if
if (RCHFLX_out(iens,iRch_ups)%Qelapsed > ntsQmodStop) then
RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%Qerror=0._dp
end if
if (RCHFLX_out(iens,iRch_ups)%Qelapsed <= ntsQmodStop) then
select case(QerrTrend)
case(const)
Qcorrect = RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror
case(linear)
Qcorrect = RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror*(1._dp - real(RCHFLX_out(iens,iRch_ups)%Qelapsed,dp)/real(ntsQmodStop, dp))
case(logistic)
x0 =0.25; y0 =0.90
k = log(1._dp/y0-1._dp)/(ntsQmodStop/2._dp-ntsQmodStop*x0)
Qcorrect = RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror/(1._dp + exp(-k*(1._dp*RCHFLX_out(iens,iRch_ups)%Qelapsed-ntsQmodStop/2._dp)))
case(exponential)
if (RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror/=0._dp) then
k = log(0.1_dp/abs(RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror))/(1._dp*ntsQmodStop)
Qcorrect = RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror*exp(k*RCHFLX_out(iens,iRch_ups)%Qelapsed)
else
Qcorrect = 0._dp
end if
case default; message=trim(message)//'discharge error trend model must be 1(const),2(liear), or 3(logistic)'; ierr=81; return
end select
else
Qcorrect=0._dp
end if
RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q = max(RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q-Qcorrect, 0._dp)
end if

q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q
end do
endif
Expand Down Expand Up @@ -257,6 +215,19 @@ SUBROUTINE dfw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be pro
write(iulog,'(A,X,G12.5,X,A,X,I9)') ' ---- NEGATIVE VOLUME (Diffusive Wave)= ', RCHFLX_out(iens,segIndex)%ROUTE(idxDW)%REACH_VOL(1), 'at ', NETOPO_in(segIndex)%REACHID
end if

if (qmodOption==1) then
call direct_insertion(iens, segIndex, & ! input: reach index
idxDW, & ! input: routing method id for diffusive wave routing
ixDesire, & ! input: verbose seg index
NETOPO_in, & ! input: reach topology data structure
RCHSTA_out, & ! inout: reach state data structure
RCHFLX_out, & ! inout: reach fluxes datq structure
ierr, cmessage) ! output: error control
if(ierr/=0)then
write(message,'(A,X,I12,X,A)') trim(message)//'/segment=', NETOPO_in(segIndex)%REACHID, '/'//trim(cmessage); return
endif
end if

END SUBROUTINE dfw_rch


Expand Down
13 changes: 6 additions & 7 deletions route/build/src/get_basin_runoff.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,6 @@ SUBROUTINE get_hru_runoff(ierr, message)
select case(qmodOption)
case(no_mod) ! do nothing
case(direct_insert)
RCHFLX(:,:)%QOBS = 0.0_dp
! read gage observation [m3/s] at current time
jx = gage_obs_data%time_ix(simDatetime(1))

Expand Down Expand Up @@ -155,7 +154,7 @@ SUBROUTINE timeMap_sim_forc(tmap_sim_forc, & ! inout: time-steps mapping dat
USE public_var, ONLY: verySmall ! smallest real values
USE public_var, ONLY: secprday ! day to second conversion factor
USE public_var, ONLY: calendar ! calender
USE public_var, ONLY: dt ! simulation time step
USE public_var, ONLY: dt_sim ! simulation time step
USE public_var, ONLY: dt_ro ! input time step

implicit none
Expand Down Expand Up @@ -189,8 +188,8 @@ SUBROUTINE timeMap_sim_forc(tmap_sim_forc, & ! inout: time-steps mapping dat

startRoSec = (juldaySim - juldayRo)*secprday ! day->sec

simLapse(1) = startRoSec+dt*(ixTime-1)
simLapse(2) = simLapse(1) + dt
simLapse(1) = startRoSec+dt_sim*(ixTime-1)
simLapse(2) = simLapse(1) + dt_sim
frcLapse = arth(0._dp, dt_ro, nRo+1)

!-- Find index of runoff time period of which end is within simulation period
Expand Down Expand Up @@ -228,11 +227,11 @@ SUBROUTINE timeMap_sim_forc(tmap_sim_forc, & ! inout: time-steps mapping dat
do iRo=idxFront, idxEnd
tmap_sim_forc%iTime(ctr) = iRo
if (iRo == idxFront)then ! front side of simulation time step
tmap_sim_forc%frac(ctr) = (frcLapse(iRo+1) - simLapse(1))/dt
tmap_sim_forc%frac(ctr) = (frcLapse(iRo+1) - simLapse(1))/dt_sim
else if ( iRo == idxEnd ) then ! end side of simulation time step
tmap_sim_forc%frac(ctr) = (simLapse(2)-frcLapse(iRo))/dt
tmap_sim_forc%frac(ctr) = (simLapse(2)-frcLapse(iRo))/dt_sim
else ! simulation time step is completely with forcing time step
tmap_sim_forc%frac(ctr) = (frcLapse(iRo+1)-frcLapse(iRo))/dt
tmap_sim_forc%frac(ctr) = (frcLapse(iRo+1)-frcLapse(iRo))/dt_sim
endif
ctr = ctr+1
end do
Expand Down
Loading

0 comments on commit 6a96294

Please sign in to comment.