diff --git a/docs/source/Control_file.rst b/docs/source/Control_file.rst index 12b01780..75b05633 100644 --- a/docs/source/Control_file.rst +++ b/docs/source/Control_file.rst @@ -65,7 +65,7 @@ The following variables (not pre-defined in the code) need to be defined in cont +--------+------------------------+--------------------------------------------------------------------------------------------------+ | 1,2,3 | | units of input runoff. e.g., mm/s | +--------+------------------------+--------------------------------------------------------------------------------------------------+ -| 1,2,3 | | time interval of input runoff in second. e.g., 86400 sec for daily step | +| 1,2,3 | | time interval of simulation time step in second. e.g., 86400 sec for daily step | +--------+------------------------+--------------------------------------------------------------------------------------------------+ | 1,2,3 | | Logical to indicate runoff needs to be remapped to RN_HRU. T or F | +--------+------------------------+--------------------------------------------------------------------------------------------------+ @@ -208,6 +208,48 @@ The output file name convension: .h.yyyy-mm-dd-sssss.nc 3. routed runoff corresponding to the scheme is not ouput if users deactivate a particular routing scheme with 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 | ++=====================+=========================================================================================================+ +| | gauge meta data (two column csv format): gauge_id (non-numeric ID is accepted), seg_id | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | gauge discharge data | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | variable name for discharge [m3/s] | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | variable name for gauge site name (character array) | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | variable name for time | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | dimension name for site | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | imension name for time | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | maximum gauge name string length | ++---------------------+---------------------------------------------------------------------------------------------------------+ + + +Data assimilation is the direct insertion that is performed at a list of reaches in the metadata. Two parameters- and are needed. + tells how bias computed at observation time at each reach evolves in the subsequent future time steps. +To activate data assimilation of observed discharge into simulated discharge, the following control variables need to be specified. + ++---------------------+---------------------------------------------------------------------------------------------------------+ +| tag | Description | ++=====================+=========================================================================================================+ +| | activation of direct insertion. 0 -> do nothing, 1=> discharge direct insertion | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | temporal discharge error trend. 1->constant, 2->linear, 3->logistic, 4->exponential | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | the number of time steps when flow correction stops | ++---------------------+---------------------------------------------------------------------------------------------------------+ + + Control file examples --------------------- diff --git a/docs/source/Input_data.rst b/docs/source/Input_data.rst index d2de7d5f..5993779b 100644 --- a/docs/source/Input_data.rst +++ b/docs/source/Input_data.rst @@ -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). diff --git a/docs/source/conf.py b/docs/source/conf.py index a5561f18..17ffcc08 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -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 diff --git a/route/build/Makefile b/route/build/Makefile index 2d93e95e..132a8525 100644 --- a/route/build/Makefile +++ b/route/build/Makefile @@ -164,6 +164,7 @@ IO = \ # CORE CORE = \ model_finalize.f90 \ + data_assimilation.f90 \ accum_runoff.f90 \ basinUH.f90 \ irf_route.f90 \ diff --git a/route/build/src/data_assimilation.f90 b/route/build/src/data_assimilation.f90 new file mode 100644 index 00000000..561c0ad4 --- /dev/null +++ b/route/build/src/data_assimilation.f90 @@ -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 diff --git a/route/build/src/dfw_route.f90 b/route/build/src/dfw_route.f90 index 79280a24..846ebf43 100644 --- a/route/build/src/dfw_route.f90 +++ b/route/build/src/dfw_route.f90 @@ -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 @@ -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 @@ -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 diff --git a/route/build/src/get_basin_runoff.f90 b/route/build/src/get_basin_runoff.f90 index 949ff64b..c3e70d6b 100644 --- a/route/build/src/get_basin_runoff.f90 +++ b/route/build/src/get_basin_runoff.f90 @@ -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)) @@ -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 @@ -176,7 +175,6 @@ SUBROUTINE timeMap_sim_forc(tmap_sim_forc, & ! inout: time-steps mapping dat real(dp) :: juldaySim ! starting julian day in simulation time step [day] integer(i4b) :: ctr ! counter integer(i4b) :: nRoSub ! - integer(i4b) :: iFile ! loop index of input file integer(i4b) :: iRo ! loop index of runoff time step integer(i4b) :: idxFront ! index of r of which top is within ith model layer (i=1..nLyr) integer(i4b) :: idxEnd ! index of the lowest soil layer of which bottom is within ith model layer (i=1..nLyr) @@ -190,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 @@ -217,7 +215,7 @@ SUBROUTINE timeMap_sim_forc(tmap_sim_forc, & ! inout: time-steps mapping dat nRoSub = idxEnd-idxFront + 1 allocate(tmap_sim_forc%iTime(nRoSub), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//'tmap_sim_forc%iTime or iFile'; return; endif + if(ierr/=0)then; message=trim(message)//trim(cmessage)//'tmap_sim_forc%iTime'; return; endif if (idxFront == idxEnd)then ! if simulation period is completely within runoff period - one runoff time period per simulation period tmap_sim_forc%iTime(ctr) = idxFront @@ -229,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 diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index bc96f836..c8b93fe6 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -3,21 +3,22 @@ MODULE irf_route_module USE nrtype ! data type USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach USE dataTypes, ONLY: RCHTOPO ! Network topology USE dataTypes, ONLY: RCHPRP ! Reach parameter USE dataTypes, ONLY: irfRCH ! irf specific state data structure - USE dataTypes,ONLY: subbasin_omp ! mainstem+tributary data structures +USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures ! global parameters USE public_var, ONLY: iulog ! i/o logical unit number USE public_var, ONLY: realMissing ! missing value for real number USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: dt ! routing time step duration [sec] +USE public_var, ONLY: dt=>dt_sim ! routing time step duration [sec] 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: nThreads ! number of threads used for openMP USE globalData, ONLY: idxIRF ! index of IRF method ! 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 @@ -34,6 +35,7 @@ SUBROUTINE irf_route(iEns, & ! input: index of runoff ensemble to be p ixDesire, & ! input: reachID to be checked by on-screen pringing NETOPO_in, & ! input: reach topology data structure RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure RCHFLX_out, & ! inout: reach flux data structure ierr, message, & ! output: error control ixSubRch) ! optional input: subset of reach indices to be processed @@ -45,7 +47,8 @@ SUBROUTINE irf_route(iEns, & ! input: index of runoff ensemble to be p integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output ! Output type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - TYPE(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + 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 integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed @@ -103,6 +106,7 @@ SUBROUTINE irf_route(iEns, & ! input: index of runoff ensemble to be p !$OMP shared(doRoute) & ! data array shared !$OMP shared(NETOPO_in) & ! data structure shared !$OMP shared(RPARAM_in) & ! data structure shared +!$OMP shared(RCHSTA_out) & ! data structure shared !$OMP shared(RCHFLX_out) & ! data structure shared !$OMP shared(ix, iEns, ixDesire) & ! indices shared !$OMP firstprivate(nTrib) @@ -116,7 +120,7 @@ SUBROUTINE irf_route(iEns, & ! input: index of runoff ensemble to be p seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) if (.not. doRoute(jSeg)) cycle - call irf_rch(iEns, jSeg, ixDesire, NETOPO_IN, RPARAM_in, RCHFLX_out, ierr, cmessage) + call irf_rch(iEns, jSeg, ixDesire, NETOPO_IN, RPARAM_in, RCHSTA_out, RCHFLX_out, ierr, cmessage) if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) end do seg ! call system_clock(openMPend(iTrib)) @@ -143,19 +147,21 @@ SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be proce ixDesire, & ! input: reachID to be checked by on-screen pringing NETOPO_in, & ! input: reach topology data structure RPARAM_in, & ! input: reach parameter 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) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - 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 + integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed + integer(i4b), intent(in) :: segIndex ! segment where routing is performed + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter + 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 real(dp) :: q_upstream ! total discharge at top of the reach being processed real(dp) :: WB_error ! water balance error [m3/s] @@ -171,14 +177,10 @@ SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be proce ! initialize future discharge array at first time if (.not.allocated(RCHFLX_out(iens,segIndex)%QFUTURE_IRF))then - - ntdh = size(NETOPO_in(segIndex)%UH) - - allocate(RCHFLX_out(iens,segIndex)%QFUTURE_IRF(ntdh), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': RCHFLX_out(iens,segIndex)%QFUTURE_IRF'; return; endif - - RCHFLX_out(iens,segIndex)%QFUTURE_IRF(:) = 0._dp - + ntdh = size(NETOPO_in(segIndex)%UH) + allocate(RCHFLX_out(iens,segIndex)%QFUTURE_IRF(ntdh), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//': RCHFLX_out(iens,segIndex)%QFUTURE_IRF'; return; endif + RCHFLX_out(iens,segIndex)%QFUTURE_IRF(:) = 0._dp end if ! get discharge coming from upstream @@ -187,17 +189,6 @@ SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be proce if (nUps>0) then do iUps = 1,nUps 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 ! there is observation - RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%Qerror = RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%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(idxIRF)%Qerror=0._dp - end if - RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%REACH_Q = max(RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%REACH_Q-RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%Qerror, 0.0001) - end if - q_upstream = q_upstream + RCHFLX_out(iens,iRch_ups)%ROUTE(idxIRF)%REACH_Q end do endif @@ -237,6 +228,20 @@ SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be proce write(iulog,'(A,X,G12.5,X,A,X,I9)') ' ---- NEGATIVE VOLUME [m3]= ', RCHFLX_out(iens,segIndex)%ROUTE(idxIRF)%REACH_VOL(1), 'at ', NETOPO_in(segIndex)%REACHID ! RCHFLX_out(iens,segIndex)%ROUTE(idxIRF)%REACH_VOL(1) = 0._dp end if + + if (qmodOption==1) then + call direct_insertion(iens, segIndex, & ! input: reach index + idxIRF, & ! input: routing method id for Euler kinematic wave + 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 irf_rch diff --git a/route/build/src/kw_route.f90 b/route/build/src/kw_route.f90 index fe704228..ae0ef060 100644 --- a/route/build/src/kw_route.f90 +++ b/route/build/src/kw_route.f90 @@ -15,10 +15,10 @@ MODULE kw_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: idxKW ! index of kw routing ! 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 @@ -177,16 +177,6 @@ SUBROUTINE kw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc 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 ! there is observation - RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%Qerror = RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%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(idxKW)%Qerror=0._dp - end if - RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%REACH_Q = max(RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%REACH_Q-RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%Qerror, 0.0001) - end if q_upstream = q_upstream + RCHFLX_out(iens,iRch_ups)%ROUTE(idxKW)%REACH_Q end do endif @@ -228,6 +218,18 @@ SUBROUTINE kw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc write(iulog,'(A,X,G12.5,X,A,X,I9)') ' ---- NEGATIVE VOLUME (Kinematic Wave)= ', RCHFLX_out(iens,segIndex)%ROUTE(idxKW)%REACH_VOL(1), 'at ', NETOPO_in(segIndex)%REACHID end if + if (qmodOption==1) then + call direct_insertion(iens, segIndex, & ! input: reach index + idxKW, & ! input: routing method id for Euler kinematic wave + 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 kw_rch diff --git a/route/build/src/main_route.f90 b/route/build/src/main_route.f90 index 583680cf..b3893729 100644 --- a/route/build/src/main_route.f90 +++ b/route/build/src/main_route.f90 @@ -145,6 +145,7 @@ SUBROUTINE main_route(iens, & ! input: ensemble index ixPrint, & ! input: index of the desired reach NETOPO, & ! input: reach topology data structure RPARAM, & ! input: reach parameter data structure + RCHSTA, & ! inout: reach state data structure RCHFLX, & ! inout: reach flux data structure ierr,cmessage) ! output: error control if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif diff --git a/route/build/src/mc_route.f90 b/route/build/src/mc_route.f90 index 3b93dba2..51086a2c 100644 --- a/route/build/src/mc_route.f90 +++ b/route/build/src/mc_route.f90 @@ -15,10 +15,10 @@ MODULE mc_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: idxMC ! index of IRF method ! 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 @@ -175,16 +175,6 @@ SUBROUTINE mc_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc 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 ! there is observation - RCHFLX_out(iens,iRch_ups)%ROUTE(idxMC)%Qerror = RCHFLX_out(iens,iRch_ups)%ROUTE(idxMC)%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(idxMC)%Qerror=0._dp - end if - RCHFLX_out(iens,iRch_ups)%ROUTE(idxMC)%REACH_Q = max(RCHFLX_out(iens,iRch_ups)%ROUTE(idxMC)%REACH_Q-RCHFLX_out(iens,iRch_ups)%ROUTE(idxMC)%Qerror, 0.0001) - end if q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%ROUTE(idxMC)%REACH_Q end do endif @@ -226,6 +216,19 @@ SUBROUTINE mc_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc write(iulog,'(A,X,G12.5,X,A,X,I9)') ' ---- NEGATIVE VOLUME (Muskingum-Cunge)= ', RCHFLX_out(iens,segIndex)%ROUTE(idxMC)%REACH_VOL(1), 'at ', NETOPO_in(segIndex)%REACHID end if + if (qmodOption==1) then + call direct_insertion(iens, segIndex, & ! input: reach index + idxMC, & ! input: routing method id for Muskingum-cunge + 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 mc_rch diff --git a/route/build/src/model_setup.f90 b/route/build/src/model_setup.f90 index 24284a06..e5a87000 100644 --- a/route/build/src/model_setup.f90 +++ b/route/build/src/model_setup.f90 @@ -190,7 +190,7 @@ END SUBROUTINE init_data ! ********************************************************************* SUBROUTINE update_time(finished, ierr, message) - USE public_var, ONLY : dt + USE public_var, ONLY : dt_sim USE public_var, ONLY : calendar USE public_var, ONLY : time_units ! netcdf time units - t_unit since yyyy-mm-dd hh:mm:ss USE globalData, ONLY : iTime ! current simulation time step index @@ -222,23 +222,23 @@ SUBROUTINE update_time(finished, ierr, message) endif ! update model time step bound - TSEC(0) = TSEC(0) + dt - TSEC(1) = TSEC(0) + dt + TSEC(0) = TSEC(0) + dt_sim + TSEC(1) = TSEC(0) + dt_sim ! update model time index iTime=iTime+1 ! increment simulation datetime simDatetime(0) = simDatetime(1) - simDatetime(1) = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + simDatetime(1) = simDatetime(1)%add_sec(dt_sim, calendar, ierr, cmessage) ! model time stamp variable for output t_unit = trim( time_units(1:index(time_units,' ')) ) select case( trim(t_unit) ) - case('seconds','second','sec','s'); timeVar = timeVar+ dt - case('minutes','minute','min'); timeVar = timeVar+ dt/60._dp - case('hours','hour','hr','h'); timeVar = timeVar+ dt/3600._dp - case('days','day','d'); timeVar = timeVar+ dt/86400._dp + case('seconds','second','sec','s'); timeVar = timeVar+ dt_sim + case('minutes','minute','min'); timeVar = timeVar+ dt_sim/60._dp + case('hours','hour','hr','h'); timeVar = timeVar+ dt_sim/3600._dp + case('days','day','d'); timeVar = timeVar+ dt_sim/86400._dp case default ierr=20; message=trim(message)//'= '//trim(t_unit)//': must be seconds, minutes, hours or days.'; return end select @@ -253,7 +253,7 @@ SUBROUTINE init_state(ierr, message) USE ascii_util_module, ONLY : lower ! convert string to lower case USE read_restart, ONLY : read_state_nc ! read netcdf state output file - USE public_var, ONLY : dt ! simulation time step (seconds) + USE public_var, ONLY : dt_sim ! simulation time step (seconds) USE public_var, ONLY : impulseResponseFunc ! IRF routing ID = 1 USE public_var, ONLY : kinematicWaveTracking ! KWT routing ID = 2 USE public_var, ONLY : kinematicWave ! KW routing ID = 3 @@ -290,6 +290,7 @@ SUBROUTINE init_state(ierr, message) RCHFLX(:,:)%BASIN_QR(0) = 0._dp RCHFLX(:,:)%BASIN_QR(1) = 0._dp RCHFLX(:,:)%Qelapsed = 0 + RCHFLX(:,:)%QOBS = 0._dp do iRoute = 1, nRoutes if (routeMethods(iRoute)==impulseResponseFunc) then @@ -333,7 +334,7 @@ SUBROUTINE init_state(ierr, message) end do ! initialize time - TSEC(0)=0._dp; TSEC(1)=dt + TSEC(0)=0._dp; TSEC(1)=dt_sim else call read_state_nc(trim(restart_dir)//trim(fname_state_in), T0, T1, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -361,13 +362,15 @@ SUBROUTINE init_time(nRoTime, & ! input: number of time steps USE public_var, ONLY: simStart ! date string defining the start of the simulation USE public_var, ONLY: simEnd ! date string defining the end of the simulation USE public_var, ONLY: calendar ! calendar name - USE public_var, ONLY: dt + USE public_var, ONLY: dt_sim ! simulation time step [sec] + USE public_var, ONLY: dt_ro ! runoff input time step [sec] USE public_var, ONLY: secprday USE public_var, ONLY: restart_write ! restart write option USE public_var, ONLY: restart_date ! restart date USE public_var, ONLY: restart_month ! USE public_var, ONLY: restart_day ! USE public_var, ONLY: restart_hour ! + USE public_var, ONLY: maxTimeDiff ! time difference tolerance for input checks USE globalData, ONLY: timeVar ! model time variables in time unit since reference time USE globalData, ONLY: iTime ! time index at simulation time step USE globalData, ONLY: simDatetime ! current model time data (yyyy:mm:dd:hh:mm:sec) @@ -394,6 +397,8 @@ SUBROUTINE init_time(nRoTime, & ! input: number of time steps real(dp) :: convTime2sec real(dp) :: roTimeVar(nRoTime) real(dp) :: sec(nRoTime) + real(dp) :: sec1(nRoTime) + real(dp) :: dt_ro_array(nRoTime-1) character(len=7) :: t_unit character(len=strLen) :: cmessage ! error message of downwind routine character(len=50) :: fmt1='(a,I4,a,I2.2,a,I2.2,x,I2.2,a,I2.2,a,F5.2)' @@ -427,8 +432,23 @@ SUBROUTINE init_time(nRoTime, & ! input: number of time steps call endDatetime%str2datetime(simEnd, ierr, cmessage) if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [endDatetime]'; return; endif - ! calendar in runoff data sec(:) = roTimeVar(:)*convTime2sec + + ! Get input (runoff) time step [sec] from input data + if (abs(dt_ro-realMissing)<=epsilon(dt_ro)) then + dt_ro = sec(2)-sec(1) + write(iulog,'(2a,x,F9.2)') new_line('a'),'INFO: input time step [s]:',dt_ro + + ! check runoff time interval is consistent + sec1 = cshift(sec, 1) + dt_ro_array = sec1(1:nRoTime-1) - sec(1:nRoTime-1) + if (any(abs(dt_ro_array-dt_ro)>maxTimeDiff)) then + write(iulog,'(2a)') new_line('a'),'WARNING: time step [s] in runoff input is not consistent' + write(iulog,'(2a)') new_line('a'),' expect the time step to be equal in the runoff time series' + end if + end if + + ! runoff data datetime [yyyy-mm-dd hh:mm:ss] do ix = 1, nRoTime roCal(ix) = refCal%add_sec(sec(ix), calendar, ierr, cmessage) end do @@ -515,12 +535,12 @@ SUBROUTINE init_time(nRoTime, & ! input: number of time steps end if call restDatetime%str2datetime(restart_date, ierr, cmessage) if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restart_date]'; return; endif - dropDatetime = restDatetime%add_sec(-dt, calendar, ierr, cmessage) + dropDatetime = restDatetime%add_sec(-dt_sim, calendar, ierr, cmessage) if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restDatetime->dropDatetime]'; return; endif restart_month = dropDatetime%month(); restart_day = dropDatetime%day(); restart_hour = dropDatetime%hour() case('yearly','monthly','daily') restDatetime = datetime(2000, restart_month, restart_day, restart_hour, 0, 0._dp) - dropDatetime = restDatetime%add_sec(-dt, calendar, ierr, cmessage) + dropDatetime = restDatetime%add_sec(-dt_sim, calendar, ierr, cmessage) if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [ dropDatetime for periodical restart]'; return; endif restart_month = dropDatetime%month(); restart_day = dropDatetime%day(); restart_hour = dropDatetime%hour() case('never') @@ -708,7 +728,6 @@ SUBROUTINE init_runoff(remap_flag, & ! input: logical whether or not runno USE public_var, ONLY : fname_remap ! name of runoff mapping netCDF name USE public_var, ONLY : calendar ! name of calendar USE public_var, ONLY : time_units ! time units - USE public_var, ONLY : realMissing ! real missing value (-9999._dp) USE dataTypes, ONLY : remap ! remapping data type USE dataTypes, ONLY : runoff ! runoff data type USE read_runoff, ONLY : read_runoff_metadata ! read meta data from runoff data @@ -850,6 +869,7 @@ SUBROUTINE init_qmod(ierr, message) character(*), intent(out) :: message ! error message ! local variables character(len=strLen) :: cmessage ! error message from subroutine + logical(lgt) :: fileExist ! file exists or not integer(i4b), parameter :: no_mod=0 integer(i4b), parameter :: direct_insert=1 @@ -863,15 +883,20 @@ SUBROUTINE init_qmod(ierr, message) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! initialize gage obs data - gage_obs_data = gageObs(trim(ancil_dir)//trim(fname_gageObs), & - vname_gageFlow, & - vname_gageTime, vname_gageSite, & - dname_gageTime, dname_gageSite, & - ierr, cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! compute link between gage ID and reach ID (river network domain) - index of reachID for each gage ID - call gage_obs_data%comp_link(reachID, gage_meta_data) + inquire(file=trim(ancil_dir)//trim(fname_gageObs), exist=fileExist) + if (fileExist) then + gage_obs_data = gageObs(trim(ancil_dir)//trim(fname_gageObs), & + vname_gageFlow, & + vname_gageTime, vname_gageSite, & + dname_gageTime, dname_gageSite, & + ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ! compute link between gage ID and reach ID (river network domain) - index of reachID for each gage ID + call gage_obs_data%comp_link(reachID, gage_meta_data) + else + qmodOption=no_mod + end if case default ierr=1; message=trim(message)//"Error: qmodOption invalid"; return end select diff --git a/route/build/src/process_ntopo.f90 b/route/build/src/process_ntopo.f90 index 0b4c3da1..c4af4a2a 100644 --- a/route/build/src/process_ntopo.f90 +++ b/route/build/src/process_ntopo.f90 @@ -65,7 +65,7 @@ SUBROUTINE augment_ntopo(nHRU, & ! input: number of HRUs USE routing_param, ONLY : make_uh ! construct reach unit hydrograph USE globalData, ONLY : mann_n, wscale ! KWT routing parameters (Transfer function parameters) USE globalData, ONLY : velo, diff ! IRF routing parameters (Transfer function parameters) - USE public_var, ONLY : dt ! simulation time step [sec] + USE public_var, ONLY : dt_sim ! simulation time step [sec] implicit none ! Argument variables @@ -243,7 +243,7 @@ SUBROUTINE augment_ntopo(nHRU, & ! input: number of HRUs end do ! compute lag times in the channel unit hydrograph - call make_uh(seg_length, dt, velo, diff, temp_dat, ierr, cmessage) + call make_uh(seg_length, dt_sim, velo, diff, temp_dat, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! put the lag times in the data structures @@ -371,7 +371,7 @@ SUBROUTINE put_data_struct(nSeg, structSEG, structNTOPO, & USE dataTypes, ONLY : RCHTOPO ! Network topology USE globalData, ONLY : fshape, tscale ! basin IRF routing parameters (Transfer function parameters) USE public_var, ONLY : min_slope ! minimum slope - USE public_var, ONLY : dt ! simulation time step [sec] + USE public_var, ONLY : dt_sim ! simulation time step [sec] USE routing_param, ONLY : basinUH ! construct basin unit hydrograph implicit none @@ -392,7 +392,7 @@ SUBROUTINE put_data_struct(nSeg, structSEG, structNTOPO, & ierr=0; message='put_data_struct/' ! get lag times in the basin unit hydrograph (not sure this is right place...) - call basinUH(dt, fshape, tscale, ierr, cmessage) + call basinUH(dt_sim, fshape, tscale, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! allocate space diff --git a/route/build/src/public_var.f90 b/route/build/src/public_var.f90 index 5443ee5a..36362694 100644 --- a/route/build/src/public_var.f90 +++ b/route/build/src/public_var.f90 @@ -26,6 +26,7 @@ MODULE public_var real(dp), parameter,public :: min_slope=1.e-6_dp ! minimum slope real(dp), parameter,public :: runoffMin=1.e-15_dp ! minimum runoff from each basin real(dp), parameter,public :: negRunoffTol=-1.e-3_dp ! nagative runoff tolerance + real(dp), parameter,public :: maxTimeDiff=1/secprday ! time difference tolerance for input checks ! routing related constants integer(i4b),parameter,public :: MAXQPAR=20 ! maximum number of particles @@ -113,7 +114,7 @@ MODULE public_var character(len=10) ,public :: routOpt = '0' ! routing scheme options 0: accum runoff, 1:IRF, 2:KWT, 3:KW, 4:MC, 5:DW integer(i4b) ,public :: doesBasinRoute = 1 ! basin routing options 0-> no, 1->IRF, otherwise error character(len=strLen),public :: newFileFrequency = 'yearly' ! frequency for new output files (daily, monthly, yearly, single) - real(dp) ,public :: dt = realMissing ! simulation time step (seconds) + real(dp) ,public :: dt_sim = realMissing ! simulation time step (seconds) ! STATES character(len=strLen),public :: restart_write = 'never' ! restart write option: never-> never write, last -> write at last time step, specified, yearly, monthly, daily character(len=strLen),public :: restart_date = charMissing ! specifed restart date @@ -142,7 +143,8 @@ MODULE public_var character(len=strLen),public :: dname_wtTime = '' ! dimension name for time ! USER OPTIONS integer(i4b) ,public :: qmodOption = 0 ! options for streamflow modification (DA): 0-> no DA, 1->direct insertion - integer(i4b) ,public :: ntsQmodStop = 10 ! number of time steps for which streamflow modification is performed + integer(i4b) ,public :: QerrTrend = 1 ! temporal discharge error decreasing trend: 1->constant, 2->linear, 3->logistic, 4->exponential + integer(i4b) ,public :: qBlendPeriod = 10 ! number of time steps for which streamflow modification is performed through blending observation logical(lgt) ,public :: takeWater = .false. ! switch for water abstraction and injection integer(i4b) ,public :: hydGeometryOption = compute ! option for hydraulic geometry calculations (0=read from file, 1=compute) integer(i4b) ,public :: topoNetworkOption = compute ! option for network topology calculations (0=read from file, 1=compute) diff --git a/route/build/src/read_control.f90 b/route/build/src/read_control.f90 index ef5fcca3..025b28b6 100644 --- a/route/build/src/read_control.f90 +++ b/route/build/src/read_control.f90 @@ -139,7 +139,7 @@ SUBROUTINE read_control(ctl_fname, err, message) case(''); routOpt = trim(cData) ! routing scheme options 0-> accumRunoff, 1->IRF, 2->KWT, 3-> KW, 4->MC, 5->DW case(''); read(cData,*,iostat=io_error) doesBasinRoute ! basin routing options 0-> no, 1->IRF, otherwise error case(''); newFileFrequency = trim(cData) ! frequency for new output options (case-insensitive): daily, monthly, yearly, or single - case(''); read(cData,*,iostat=io_error) dt ! time interval of the simulation (To-do: change dt to dt_sim) + case(''); read(cData,*,iostat=io_error) dt_sim ! time interval of the simulation ! RESTART case(''); restart_write = trim(cData) ! restart write option (case-insensitive): never, last, specified, yearly, monthly, or daily case(''); restart_date = trim(cData) ! specified restart date, yyyy-mm-dd (hh:mm:ss) for Specified option @@ -151,7 +151,8 @@ SUBROUTINE read_control(ctl_fname, err, message) case(''); param_nml = trim(cData) ! name of namelist including routing parameter value ! USER OPTIONS: Define options to include/skip calculations case(''); read(cData,*,iostat=io_error) qmodOption ! option for streamflow modification (DA): 0->no DA, 1->direct insertion - case(''); read(cData,*,iostat=io_error) ntsQmodStop ! number of time steps for which streamflow modification is performed + case(''); read(cData,*,iostat=io_error) qBlendPeriod ! number of time steps for which streamflow modification is performed through blending observation + case(''); read(cData,*,iostat=io_error) QerrTrend ! temporal dischargge error decreasing trend. 1->constant, 2->linear, 3->logistic, 4->exponential case(''); read(cData,*,iostat=io_error) takeWater ! switch for abstraction/injection case(''); read(cData,*,iostat=io_error) hydGeometryOption ! option for hydraulic geometry calculations (0=read from file, 1=compute) case(''); read(cData,*,iostat=io_error) topoNetworkOption ! option for network topology calculations (0=read from file, 1=compute) diff --git a/route/build/src/write_restart.f90 b/route/build/src/write_restart.f90 index a336f476..acda81b6 100644 --- a/route/build/src/write_restart.f90 +++ b/route/build/src/write_restart.f90 @@ -17,7 +17,7 @@ MODULE write_restart USE public_var, ONLY: iulog ! i/o logical unit number USE public_var, ONLY: integerMissing USE public_var, ONLY: realMissing -USE public_var, ONLY: dt +USE public_var, ONLY: dt_sim USE public_var, ONLY: doesBasinRoute USE public_var, ONLY: impulseResponseFunc USE public_var, ONLY: kinematicWaveTracking @@ -147,8 +147,8 @@ SUBROUTINE restart_output(ierr, message) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! update model time step bound - TSEC1 = TSEC(0) + dt - TSEC2 = TSEC1 + dt + TSEC1 = TSEC(0) + dt_sim + TSEC2 = TSEC1 + dt_sim call write_state_nc(fnameRestart, & ! Input: state netcdf name TSEC1, TSEC2, & ! Input: time, time step, start and end time [sec] @@ -189,7 +189,7 @@ SUBROUTINE restart_fname(fnameRestart, timeStamp, ierr, message) select case(timeStamp) case(currTimeStep); timeStampCal = simDatetime(1) - case(nextTimeStep); timeStampCal = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + case(nextTimeStep); timeStampCal = simDatetime(1)%add_sec(dt_sim, calendar, ierr, cmessage) case default; ierr=20; message=trim(message)//'time stamp option in restart filename: invalid -> 1: current time Step or 2: next time step'; return end select @@ -270,7 +270,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get which time is for restarting - timeStampCal = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + timeStampCal = simDatetime(1)%add_sec(dt_sim, calendar, ierr, cmessage) write(globalDesc, fmtYMDHMS) timeStampCal%year(),'-',timeStampCal%month(),'-',timeStampCal%day(),timeStampCal%hour(),':',timeStampCal%minute(),':',nint(timeStampCal%sec()) call put_global_attr(ncid, 'Restart time', trim(globalDesc), ierr, cmessage)