diff --git a/route/build/src/dataTypes.f90 b/route/build/src/dataTypes.f90 index e1512422..59135db4 100644 --- a/route/build/src/dataTypes.f90 +++ b/route/build/src/dataTypes.f90 @@ -351,6 +351,7 @@ MODULE dataTypes ! ---------- reach fluxes -------------------------------------------------------------------- type, public :: hydraulic real(dp) :: REACH_ELE ! water height at current time step [m] + real(dp) :: REACH_INFLOW ! total upstream discharge at current time step [m3/s] real(dp) :: REACH_Q ! discharge at current time step [m3/s] real(dp) :: REACH_VOL(0:1) ! water volume at previous and current time steps [m3] real(dp) :: REACH_WM_FLUX_actual ! water management fluxes to and from each reach [m3/s] diff --git a/route/build/src/dfw_route.f90 b/route/build/src/dfw_route.f90 index 1be0b4af..20ff5ed5 100644 --- a/route/build/src/dfw_route.f90 +++ b/route/build/src/dfw_route.f90 @@ -91,6 +91,8 @@ SUBROUTINE dfw_rch(this, & ! dfw_route_rch object to bound this proced end do endif + RCHFLX_out(iens,segIndex)%ROUTE(idxDW)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches + q_upstream_mod = q_upstream Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1) Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative) diff --git a/route/build/src/histVars_data.f90 b/route/build/src/histVars_data.f90 index a631a603..07eba0af 100644 --- a/route/build/src/histVars_data.f90 +++ b/route/build/src/histVars_data.f90 @@ -12,6 +12,7 @@ MODULE histVars_data USE public_var, ONLY: kinematicWave ! KW routing ID = 3 USE public_var, ONLY: muskingumCunge ! MC routing ID = 4 USE public_var, ONLY: diffusiveWave ! DW routing ID = 5 + USE public_var, ONLY: outputInflow ! logical for outputting upstream inflow in history file USE globalData, ONLY: nRoutes USE globalData, ONLY: routeMethods USE globalData, ONLY: meta_rflx, meta_hflx @@ -54,6 +55,7 @@ MODULE histVars_data real(dp), allocatable :: instRunoff(:) ! instantaneous lateral inflow into each river/lake [m3/s] [nRch] real(dp), allocatable :: dlayRunoff(:) ! lateral inflow into river or lake [m3/s] for each reach [nRch] real(dp), allocatable :: discharge(:,:) ! river/lake discharge [m3/s] for each reach/lake and routing method [nRch,nMethod] + real(dp), allocatable :: inflow(:,:) ! inflow from upstream rivers/lakes [m3/s] for each reach/lake and routing method [nRch,nMethod] real(dp), allocatable :: elev(:,:) ! river/lake surface water elevation [m] for each reach/lake and routing method [nRch,nMethod] real(dp), allocatable :: volume(:,:) ! river/lake volume [m3] for each reach/lake and routing method [nRch,nMethod] @@ -95,31 +97,31 @@ FUNCTION constructor(nHru_local, nRch_local, ierr, message) RESULT(instHistVar) instHistVar%nRch = nRch_local if (meta_hflx(ixHFLX%basRunoff)%varFile) then - allocate(instHistVar%basRunoff(nHRU_local), stat=ierr, errmsg=cmessage) + allocate(instHistVar%basRunoff(nHRU_local), source=0._dp, stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%basRunoff]'; return; endif - instHistVar%basRunoff(1:nHRU_local) = 0._dp end if if (meta_rflx(ixRFLX%instRunoff)%varFile) then - allocate(instHistVar%instRunoff(nRch_local), stat=ierr, errmsg=cmessage) + allocate(instHistVar%instRunoff(nRch_local), source=0._dp, stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif - instHistVar%instRunoff(1:nRch_local) = 0._dp end if if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then - allocate(instHistVar%dlayRunoff(nRch_local), stat=ierr, errmsg=cmessage) + allocate(instHistVar%dlayRunoff(nRch_local), source=0._dp, stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif - instHistVar%dlayRunoff(1:nRch_local) = 0._dp end if if (nRoutes>0) then ! this should be number of methods that ouput - allocate(instHistVar%discharge(nRch_local, nRoutes), stat=ierr, errmsg=cmessage) + allocate(instHistVar%discharge(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%discharge]'; return; endif - instHistVar%discharge(1:nRch_local, 1:nRoutes) = 0._dp - allocate(instHistVar%volume(nRch_local, nRoutes), stat=ierr, errmsg=cmessage) + allocate(instHistVar%volume(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%volume]'; return; endif - instHistVar%volume(1:nRch_local, 1:nRoutes) = 0._dp + + if (outputInflow) then + allocate(instHistVar%inflow(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%inflow]'; return; endif + end if end if END FUNCTION constructor @@ -205,6 +207,9 @@ SUBROUTINE aggregate(this, & ! inout: do ix=1,this%nRch this%discharge(ix,iRoute) = this%discharge(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_Q this%volume(ix,iRoute) = this%volume(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_VOL(1) + if (outputInflow) then + this%inflow(ix,iRoute) = this%inflow(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_INFLOW + end if end do end do @@ -245,6 +250,11 @@ SUBROUTINE finalize(this) this%volume = this%volume/real(this%nt, kind=dp) end if + ! 6. inflow + if (allocated(this%inflow)) then + this%inflow = this%inflow/real(this%nt, kind=dp) + end if + END SUBROUTINE finalize ! --------------------------------------------------------------- @@ -263,6 +273,7 @@ SUBROUTINE refresh(this) if (allocated(this%dlayRunoff)) this%dlayRunoff = 0._dp if (allocated(this%discharge)) this%discharge = 0._dp if (allocated(this%volume)) this%volume = 0._dp + if (allocated(this%inflow)) this%inflow = 0._dp END SUBROUTINE refresh @@ -280,6 +291,7 @@ SUBROUTINE clean(this) if (allocated(this%dlayRunoff)) deallocate(this%dlayRunoff) if (allocated(this%discharge)) deallocate(this%discharge) if (allocated(this%volume)) deallocate(this%volume) + if (allocated(this%inflow)) deallocate(this%inflow) END SUBROUTINE clean @@ -299,6 +311,7 @@ SUBROUTINE read_restart(this, restart_name, ierr, message) real(dp), allocatable :: array_tmp(:) ! temp array integer(i4b) :: ixRoute ! loop index integer(i4b) :: ixFlow, ixVol ! temporal method index + integer(i4b) :: ixInflow ! temporal method index logical(lgt) :: FileStatus ! file open or close type(file_desc_t) :: pioFileDesc ! pio file handle @@ -369,11 +382,16 @@ SUBROUTINE read_restart(this, restart_name, ierr, message) allocate(this%discharge(this%nRch, nRoutes), stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%discharge]'; return; endif allocate(this%volume(this%nRch, nRoutes), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%discharge]'; return; endif + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%volume]'; return; endif this%discharge = 0._dp this%volume = 0._dp + if (outputInflow) then + allocate(this%inflow(this%nRch, nRoutes), source=0.0_dp, stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%volume]'; return; endif + end if + do ixRoute=1,nRoutes select case(routeMethods(ixRoute)) case(accumRunoff) @@ -381,18 +399,23 @@ SUBROUTINE read_restart(this, restart_name, ierr, message) case(impulseResponseFunc) ixFlow=ixRFLX%IRFroutedRunoff ixVol=ixRFLX%IRFvolume + ixInflow=ixRFLX%IRFinflow case(kinematicWaveTracking) ixFlow=ixRFLX%KWTroutedRunoff ixVol=ixRFLX%KWTvolume + ixInflow=ixRFLX%KWTinflow case(kinematicWave) ixFlow=ixRFLX%KWroutedRunoff ixVol=ixRFLX%KWvolume + ixInflow=ixRFLX%KWinflow case(muskingumCunge) ixFlow=ixRFLX%MCroutedRunoff ixVol=ixRFLX%MCvolume + ixInflow=ixRFLX%MCinflow case(diffusiveWave) ixFlow=ixRFLX%DWroutedRunoff ixVol=ixRFLX%DWvolume + ixInflow=ixRFLX%DWinflow case default write(message,'(2A,1X,G0,1X,A)') trim(message), 'routing method index:',routeMethods(ixRoute), 'must be 0-5' ierr=81; return @@ -418,13 +441,25 @@ SUBROUTINE read_restart(this, restart_name, ierr, message) ! need to shift tributary part in main core to after halo reaches (nTribOutlet) if (masterproc) then this%volume(1:nRch_mainstem, ixRoute) = array_tmp(1:nRch_mainstem) - this%volume(nRch_mainstem+nTribOutlet+1:this%nRch,ixRoute) = this%volume(nRch_mainstem+1:nRch_mainstem+nRch_trib,ixRoute) + this%volume(nRch_mainstem+nTribOutlet+1:this%nRch, ixRoute) = array_tmp(nRch_mainstem+1:nRch_mainstem+nRch_trib) else this%volume(:,ixRoute) = array_tmp end if end if - end do - end if + + if (meta_rflx(ixInflow)%varFile) then + call read_dist_array(pioFileDesc, meta_rflx(ixInflow)%varName, array_tmp, ioDesc_hist_rch_double, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + ! need to shift tributary part in main core to after halo reaches (nTribOutlet) + if (masterproc) then + this%inflow(1:nRch_mainstem, ixRoute) = array_tmp(1:nRch_mainstem) + this%inflow(nRch_mainstem+nTribOutlet+1:this%nRch, ixRoute) = array_tmp(nRch_mainstem+1:nRch_mainstem+nRch_trib) + else + this%inflow(:,ixRoute) = array_tmp + end if + end if + end do ! end of ixRoute loop + end if ! end of nRoute>0 if-statement call closeFile(pioFileDesc, FileStatus) diff --git a/route/build/src/historyFile.f90 b/route/build/src/historyFile.f90 index ea03ff8f..a3e1da25 100644 --- a/route/build/src/historyFile.f90 +++ b/route/build/src/historyFile.f90 @@ -462,6 +462,16 @@ SUBROUTINE write_flux_rch(this, hVars_local, ioDescRchFlux, index_write, ierr, m array_temp(1:nRch_write) = hVars_local%volume(index_write, idxMC) case(ixRFLX%DWvolume) array_temp(1:nRch_write) = hVars_local%volume(index_write, idxDW) + case(ixRFLX%KWTinflow) + array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxKWT) + case(ixRFLX%IRFinflow) + array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxIRF) + case(ixRFLX%KWinflow) + array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxKW) + case(ixRFLX%MCinflow) + array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxMC) + case(ixRFLX%DWinflow) + array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxDW) case default write(message,'(2A,1X,G0,1X,1A)') trim(message), 'Invalid RFLX variable index:',iVar,'. Check ixRFLX in var_lookup.f90' ierr=81; return diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index bec98ad7..a19355ae 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -85,15 +85,18 @@ SUBROUTINE irf_rch(this, & ! irf_route_rch object to bound this procedur end if ! get discharge coming from upstream - nUps = size(NETOPO_in(segIndex)%UREACHI) + nUps = count(NETOPO_in(segIndex)%goodBas) ! reminder: goodBas is reach with >0 total contributory area q_upstream = 0.0_dp if (nUps>0) then do iUps = 1,nUps - iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach + if (.not. NETOPO_in(segIndex)%goodBas(iUps)) cycle ! skip upstream reach which does not any flow due to zero total contributory areas + iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%ROUTE(idxIRF)%REACH_Q end do endif + RCHFLX_out(iens,segIndex)%ROUTE(idxIRF)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches + q_upstream_mod = q_upstream Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1) Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative) diff --git a/route/build/src/kwe_route.f90 b/route/build/src/kwe_route.f90 index b431b27a..bca81ea8 100644 --- a/route/build/src/kwe_route.f90 +++ b/route/build/src/kwe_route.f90 @@ -92,6 +92,8 @@ SUBROUTINE kw_rch(this, & ! kwe_route_rch object to bound this procedu end do endif + RCHFLX_out(iens,segIndex)%ROUTE(idxKW)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches + q_upstream_mod = q_upstream Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1) Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative) diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index 1eb31186..851341a3 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -110,9 +110,12 @@ SUBROUTINE kwt_rch(this, & ! kwt_route_rch object to bound this procedur ! local variables ! (1) extract flow from upstream reaches and append to the non-routed flow in JRCH integer(i4b) :: NUPS ! number of upstream reaches + integer(i4b) :: iUps ! upstream loop index + integer(i4b) :: iRch_ups ! upstream reach index real(dp),dimension(:),allocatable :: Q_JRCH ! flow in downstream reach JRCH real(dp),dimension(:),allocatable :: TENTRY ! entry time to JRCH (exit time u/s) integer(i4b) :: NQ1 ! # flow particles + real(dp) :: q_upstream ! total discharge at top of the reach [m3/s] ! (2) route flow within the current [JRCH] river segment integer(i4b) :: ROFFSET ! retrospective offset due to rstep real(dp) :: T_START ! start of time step @@ -162,12 +165,21 @@ SUBROUTINE kwt_rch(this, & ! kwt_route_rch object to bound this procedur ierr=20; message=trim(message)//'negative flow extracted from upstream reach'; return endif + q_upstream = 0._dp + do iUps = 1,nUps + if (.not. NETOPO_in(segIndex)%goodBas(iUps)) cycle ! skip upstream reach which does not any flow due to zero total contributory areas + iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach + q_upstream = q_upstream + RCHFLX_out(iens,iRch_ups)%ROUTE(idxKWT)%REACH_Q + end do + RCHFLX_out(iens,segIndex)%ROUTE(idxKWT)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches + if(segIndex==ixDesire)then write(fmt1,'(A,I5,A)') '(A, 1X',size(Q_JRCH),'(1X,G15.4))' write(iulog,'(a)') ' * Wave discharge from upstream reaches (Q_JRCH) [m2/s]:' write(iulog,fmt1) ' Q_JRCH=',(Q_JRCH(IWV), IWV=0,size(Q_JRCH)-1) endif else + RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_INFLOW = 0._dp ! set flow in headwater reaches to modelled streamflow from time delay histogram RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q = RCHFLX_out(IENS,segIndex)%BASIN_QR(1) if (allocated(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE)) then diff --git a/route/build/src/mc_route.f90 b/route/build/src/mc_route.f90 index 2e4e859d..66c9f188 100644 --- a/route/build/src/mc_route.f90 +++ b/route/build/src/mc_route.f90 @@ -91,6 +91,8 @@ SUBROUTINE mc_rch(this, & ! mc_route_rch object to bound this procedur end do endif + RCHFLX_out(iens,segIndex)%ROUTE(idxMC)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches + q_upstream_mod = q_upstream Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1) Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative) diff --git a/route/build/src/popMetadat.f90 b/route/build/src/popMetadat.f90 index 74366db5..79b72047 100644 --- a/route/build/src/popMetadat.f90 +++ b/route/build/src/popMetadat.f90 @@ -230,20 +230,25 @@ subroutine popMetadat(err,message) meta_PFAF (ixPFAF%code ) = var_info('code' , 'pfafstetter code' ,'-' ,ixDims%seg , .false.) ! ---------- populate history output segment fluxes/state metadata structures ---------------------------------------------------------------------------------------------------- - ! Reach Flux varName varDesc unit, varType, varDim, writeOut - call meta_rflx(ixRFLX%instRunoff )%init('instRunoff' , 'instantaneous runoff in each reach' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%dlayRunoff )%init('dlayRunoff' , 'delayed runoff in each reach' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%sumUpstreamRunoff)%init('sumUpstreamRunoff', 'sum of upstream runoff in each reach' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%IRFroutedRunoff )%init('IRFroutedRunoff' , 'routed runoff in each reach-impulse response function', 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%KWTroutedRunoff )%init('KWTroutedRunoff' , 'routed runoff in each reach-lagrangian kinematic wave', 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%KWroutedRunoff )%init('KWroutedRunoff' , 'routed runoff in each reach-kinematic wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%MCroutedRunoff )%init('MCroutedRunoff' , 'routed runoff in each reach-muskingum-cunge' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%DWroutedRunoff )%init('DWroutedRunoff' , 'routed runoff in each reach-diffusive wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%IRFvolume )%init('IRFvolume' , 'lake and stream volume-impulse response function' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .true.) - call meta_rflx(ixRFLX%KWTvolume )%init('KWTvolume' , 'lake and stream volume-lagrangian kinematic wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%KWvolume )%init('KWvolume' , 'lake and stream volume-kinematic wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%MCvolume )%init('MCvolume' , 'lake and stream volume-muskingum-cunge' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) - call meta_rflx(ixRFLX%DWvolume )%init('DWvolume' , 'lake and stream volume-diffusive wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) + ! Reach Flux varName varDesc unit, varType, varDim, writeOut + call meta_rflx(ixRFLX%instRunoff )%init('instRunoff' , 'instantaneous runoff into stream or lake' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%dlayRunoff )%init('dlayRunoff' , 'delayed runoff into stream or lake' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%sumUpstreamRunoff)%init('sumUpstreamRunoff', 'sum of upstream runoff in each reach' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%IRFroutedRunoff )%init('IRFroutedRunoff' , 'Discharge from lake or stream-impulse response function' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%KWTroutedRunoff )%init('KWTroutedRunoff' , 'Discharge from lake or stream-lagrangian kinematic wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%KWroutedRunoff )%init('KWroutedRunoff' , 'Discharge from lake or stream-kinematic wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%MCroutedRunoff )%init('MCroutedRunoff' , 'Discharge from lake or stream-muskingum-cunge' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%DWroutedRunoff )%init('DWroutedRunoff' , 'Discharge from lake or stream-diffusive wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%IRFvolume )%init('IRFvolume' , 'Volume in lake or stream-impulse response function' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .true.) + call meta_rflx(ixRFLX%KWTvolume )%init('KWTvolume' , 'Volume in lake or stream-lagrangian kinematic wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%KWvolume )%init('KWvolume' , 'Volume in lake or stream-kinematic wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%MCvolume )%init('MCvolume' , 'Volume in lake or stream-muskingum-cunge' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%DWvolume )%init('DWvolume' , 'Volume in lake or stream-diffusive wave' , 'm3' , ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%IRFinflow )%init('IRFinflow' , 'Inflow from upstream lake or streams-impulse response function' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%KWTinflow )%init('KWTinflow' , 'Inflow from upstream lake or streams-lagrangian kinematic wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%KWinflow )%init('KWinflow' , 'Inflow from upstream lake or streams-kinematic wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%MCinflow )%init('MCinflow' , 'Inflow from upstream lake or streams-muskingum-cunge' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) + call meta_rflx(ixRFLX%DWinflow )%init('DWinflow' , 'Inflow from upstream lake or streams-diffusive wave' , 'm3/s', ncd_float, [ixQdims%seg,ixQdims%time], .false.) ! HRU flux varName varDesc unit, varType, varDim, writeOut call meta_hflx(ixHFLX%basRunoff )%init('basRunoff' , 'basin runoff' , 'm/s' , ncd_float, [ixQdims%hru,ixQdims%time], .true.) diff --git a/route/build/src/public_var.f90 b/route/build/src/public_var.f90 index dc0c1573..c4045c70 100644 --- a/route/build/src/public_var.f90 +++ b/route/build/src/public_var.f90 @@ -177,6 +177,8 @@ MODULE public_var character(len=strLen),public :: dname_gageSite = '' ! dimension name for gauge site character(len=strLen),public :: dname_gageTime = '' ! dimension name for time integer(i4b) ,public :: strlen_gageSite = 30 ! maximum character length for site name + ! OUTPUT OPTIONS + logical(lgt) ,public :: outputInflow = .false. ! logical; T-> write upstream inflow in history file output ! USER OPTIONS integer(i4b) ,public :: qmodOption = 0 ! option for streamflow modification integer(i4b) ,public :: hydGeometryOption = compute ! option for hydraulic geometry calculations (0=read from file, 1=compute) diff --git a/route/build/src/read_control.f90 b/route/build/src/read_control.f90 index 9fe5df35..11cea995 100644 --- a/route/build/src/read_control.f90 +++ b/route/build/src/read_control.f90 @@ -230,6 +230,7 @@ SUBROUTINE read_control(ctl_fname, err, message) case(''); read(cData,*,iostat=io_error) meta_rflx(ixRFLX%KWvolume )%varFile ! default: true (turned off if inactive) case(''); read(cData,*,iostat=io_error) meta_rflx(ixRFLX%MCvolume )%varFile ! default: true (turned off if inactive) case(''); read(cData,*,iostat=io_error) meta_rflx(ixRFLX%DWvolume )%varFile ! default: true (turned off if inactive) + case(''); read(cData,*,iostat=io_error) outputInflow ! VARIABLE NAMES for data (overwrite default name in popMeta.f90) ! HRU structure @@ -522,6 +523,14 @@ SUBROUTINE read_control(ctl_fname, err, message) end do ! ---- history Output variables + if (outputInflow) then + meta_rflx(ixRFLX%KWTinflow)%varFile=.true. + meta_rflx(ixRFLX%IRFinflow)%varFile=.true. + meta_rflx(ixRFLX%MCinflow)%varFile=.true. + meta_rflx(ixRFLX%KWinflow)%varFile=.true. + meta_rflx(ixRFLX%DWinflow)%varFile=.true. + end if + ! Make sure turned off if the corresponding routing is not running do iRoute = 0, nRouteMethods-1 select case(iRoute) case(accumRunoff) @@ -532,26 +541,31 @@ SUBROUTINE read_control(ctl_fname, err, message) if (.not. onRoute(iRoute)) then meta_rflx(ixRFLX%KWTroutedRunoff)%varFile=.false. meta_rflx(ixRFLX%KWTvolume)%varFile=.false. + meta_rflx(ixRFLX%KWTinflow)%varFile=.false. end if case(impulseResponseFunc) if (.not. onRoute(iRoute)) then meta_rflx(ixRFLX%IRFroutedRunoff)%varFile=.false. meta_rflx(ixRFLX%IRFvolume)%varFile=.false. + meta_rflx(ixRFLX%IRFinflow)%varFile=.false. end if case(muskingumCunge) if (.not. onRoute(iRoute)) then meta_rflx(ixRFLX%MCroutedRunoff)%varFile=.false. meta_rflx(ixRFLX%MCvolume)%varFile=.false. + meta_rflx(ixRFLX%MCinflow)%varFile=.false. end if case(kinematicWave) if (.not. onRoute(iRoute)) then meta_rflx(ixRFLX%KWroutedRunoff)%varFile=.false. meta_rflx(ixRFLX%KWvolume)%varFile=.false. + meta_rflx(ixRFLX%KWinflow)%varFile=.false. end if case(diffusiveWave) if (.not. onRoute(iRoute)) then meta_rflx(ixRFLX%DWroutedRunoff)%varFile=.false. meta_rflx(ixRFLX%DWvolume)%varFile=.false. + meta_rflx(ixRFLX%DWinflow)%varFile=.false. end if case default; message=trim(message)//'expect digits from 0 and 5'; err=81; return end select diff --git a/route/build/src/var_lookup.f90 b/route/build/src/var_lookup.f90 index 66ae9ef8..79e552e2 100644 --- a/route/build/src/var_lookup.f90 +++ b/route/build/src/var_lookup.f90 @@ -216,6 +216,11 @@ MODULE var_lookup integer(i4b) :: KWvolume = integerMissing ! 11. KW water volume in reach/lake (m3) integer(i4b) :: MCvolume = integerMissing ! 12. MC water volume in reach/lake (m3) integer(i4b) :: DWvolume = integerMissing ! 13. DW water volume in reach/lake (m3) + integer(i4b) :: IRFinflow = integerMissing ! 14. IRF inflow from upstreams into reach/lake (m3/s) + integer(i4b) :: KWTinflow = integerMissing ! 15. KWT inflow flow upstreams into reach/lake (m3/s) + integer(i4b) :: KWinflow = integerMissing ! 16. KW inflow from upstreams into reach/lake (m3/s) + integer(i4b) :: MCinflow = integerMissing ! 17. MC inflow from upstreams into reach/lake (m3/s) + integer(i4b) :: DWinflow = integerMissing ! 18. DW inflow from upstreams into reach/lake (m3/s) endtype iLook_RFLX ! HRU fluxes type, public :: iLook_HFLX @@ -283,7 +288,7 @@ MODULE var_lookup 21,22) type(iLook_PFAF) ,public,parameter :: ixPFAF = iLook_PFAF ( 1) type(iLook_RFLX) ,public,parameter :: ixRFLX = iLook_RFLX ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & - 11,12,13) + 11,12,13,14,15,16,17,18) type(iLook_HFLX) ,public,parameter :: ixHFLX = iLook_HFLX ( 1) type(iLook_basinQ) ,public,parameter :: ixBasinQ = iLook_basinQ ( 1) type(iLook_IRFbas) ,public,parameter :: ixIRFbas = iLook_IRFbas ( 1) diff --git a/route/build/src/write_restart_pio.f90 b/route/build/src/write_restart_pio.f90 index efe6924c..e4b882b9 100644 --- a/route/build/src/write_restart_pio.f90 +++ b/route/build/src/write_restart_pio.f90 @@ -1287,6 +1287,37 @@ SUBROUTINE write_history_state(ierr, message1) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif + if (meta_rflx(ixRFLX%KWTinflow)%varFile) then + array_dp = hVars%inflow(index_write, idxKWT) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call write_pnetcdf(pioFileDesc, 'KWTinflow', array_dp, ioDesc_hist_rch_double, ierr, cmessage) + endif + + if (meta_rflx(ixRFLX%IRFinflow)%varFile) then + array_dp = hVars%inflow(index_write, idxIRF) + call write_pnetcdf(pioFileDesc, 'IRFinflow', array_dp, ioDesc_hist_rch_double, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%KWinflow)%varFile) then + array_dp = hVars%inflow(index_write, idxKW) + call write_pnetcdf(pioFileDesc, 'KWinflow', array_dp, ioDesc_hist_rch_double, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%MCinflow)%varFile) then + array_dp = hVars%inflow(index_write, idxMC) + call write_pnetcdf(pioFileDesc, 'MCinflow', array_dp, ioDesc_hist_rch_double, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%DWinflow)%varFile) then + array_dp = hVars%inflow(index_write, idxDW) + call write_pnetcdf(pioFileDesc, 'DWinflow', array_dp, ioDesc_hist_rch_double, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + END SUBROUTINE write_history_state END SUBROUTINE write_state_nc