Skip to content

Commit

Permalink
Merge pull request #358 from ESCOMP/develop
Browse files Browse the repository at this point in the history 
v1.2.3- Take 3
  • Loading branch information
@nmizukami
nmizukami authored Mar 30, 2023
2 parents 8cc5bc0 + 6a96294 commit 104ef34
Show file tree
Hide file tree
Showing 16 changed files with 316 additions and 134 deletions.
44 changes: 43 additions & 1 deletion docs/source/Control_file.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ The following variables (not pre-defined in the code) need to be defined in cont
+--------+------------------------+--------------------------------------------------------------------------------------------------+
| 1,2,3 | <units_qsim> | units of input runoff. e.g., mm/s |
+--------+------------------------+--------------------------------------------------------------------------------------------------+
| 1,2,3 | <dt_qsim> | time interval of input runoff in second. e.g., 86400 sec for daily step |
| 1,2,3 | <dt_qsim> | time interval of simulation time step in second. e.g., 86400 sec for daily step |
+--------+------------------------+--------------------------------------------------------------------------------------------------+
| 1,2,3 | <is_remap> | Logical to indicate runoff needs to be remapped to RN_HRU. T or F |
+--------+------------------------+--------------------------------------------------------------------------------------------------+
Expand Down Expand Up @@ -208,6 +208,48 @@ The output file name convension: <case_name>.h.yyyy-mm-dd-sssss.nc
3. routed runoff corresponding to the scheme is not ouput if users deactivate a particular routing scheme with <route_opt> tag.


Data assimilation options
---------------------

mizuRoute can read gauge observed discharge data (in netCDF) along with gauge meta ascii data. To read gauge observation and gauge metadata, the following control variables need to be specified.


+---------------------+---------------------------------------------------------------------------------------------------------+
| tag | Description |
+=====================+=========================================================================================================+
| <gageMetaFile> | gauge meta data (two column csv format): gauge_id (non-numeric ID is accepted), seg_id |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <fname_gageObs> | gauge discharge data |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <vname_gageFlow> | variable name for discharge [m3/s] |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <vname_gageSite> | variable name for gauge site name (character array) |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <vname_gageTime> | variable name for time |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <dname_gageSite> | dimension name for site |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <dname_gageTime> | imension name for time |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <strlen_gageSite> | maximum gauge name string length |
+---------------------+---------------------------------------------------------------------------------------------------------+


Data assimilation is the direct insertion that is performed at a list of reaches in the metadata. Two parameters-<QerrTrend> and <ntsQmodStop> are needed.
<QerrTrend> tells how bias computed at observation time at each reach evolves in the subsequent future <ntsQmodStop> time steps.
To activate data assimilation of observed discharge into simulated discharge, the following control variables need to be specified.

+---------------------+---------------------------------------------------------------------------------------------------------+
| tag | Description |
+=====================+=========================================================================================================+
| <qmodOption> | activation of direct insertion. 0 -> do nothing, 1=> discharge direct insertion |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <QerrTrend> | temporal discharge error trend. 1->constant, 2->linear, 3->logistic, 4->exponential |
+---------------------+---------------------------------------------------------------------------------------------------------+
| <ntsQmodStop> | the number of time steps when flow correction stops |
+---------------------+---------------------------------------------------------------------------------------------------------+


Control file examples
---------------------

Expand Down
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
3 changes: 2 additions & 1 deletion docs/source/conf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,8 @@
# a list of builtin themes.
#
#html_theme = 'alabaster'
html_theme = 'sphinxdoc'
#html_theme = 'sphinxdoc'
html_theme = 'sphinx_rtd_theme'

# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
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
27 changes: 15 additions & 12 deletions route/build/src/dfw_route.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +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 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 @@ -173,16 +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
RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q = max(RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q-RCHFLX_out(iens,iRch_ups)%ROUTE(idxDW)%Qerror, 0.0001)
end if
q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%ROUTE(idxDW)%REACH_Q
end do
endif
Expand Down Expand Up @@ -225,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
Loading

0 comments on commit 104ef34

Please sign in to comment.