Skip to content

Commit 2ff305a

Browse files
nmizukaminmizukami
authored
Merge pull request #462 from nmizukami/cesm-coupling_add_upstreamflow_output
Add upstream inflow in history output
2 parents 3098538 + 121918a commit 2ff305a

13 files changed

+155
-31
lines changed

‎route/build/src/dataTypes.f90

+1
Original file line numberDiff line numberDiff line change
@@ -351,6 +351,7 @@ MODULE dataTypes
351351
! ---------- reach fluxes --------------------------------------------------------------------
352352
type, public :: hydraulic
353353
real(dp) :: REACH_ELE ! water height at current time step [m]
354+
real(dp) :: REACH_INFLOW ! total upstream discharge at current time step [m3/s]
354355
real(dp) :: REACH_Q ! discharge at current time step [m3/s]
355356
real(dp) :: REACH_VOL(0:1) ! water volume at previous and current time steps [m3]
356357
real(dp) :: REACH_WM_FLUX_actual ! water management fluxes to and from each reach [m3/s]

‎route/build/src/dfw_route.f90

+2
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,8 @@ SUBROUTINE dfw_rch(this, & ! dfw_route_rch object to bound this proced
9191
end do
9292
endif
9393

94+
RCHFLX_out(iens,segIndex)%ROUTE(idxDW)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches
95+
9496
q_upstream_mod = q_upstream
9597
Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1)
9698
Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative)

‎route/build/src/histVars_data.f90

+49-14
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ MODULE histVars_data
1212
USE public_var, ONLY: kinematicWave ! KW routing ID = 3
1313
USE public_var, ONLY: muskingumCunge ! MC routing ID = 4
1414
USE public_var, ONLY: diffusiveWave ! DW routing ID = 5
15+
USE public_var, ONLY: outputInflow ! logical for outputting upstream inflow in history file
1516
USE globalData, ONLY: nRoutes
1617
USE globalData, ONLY: routeMethods
1718
USE globalData, ONLY: meta_rflx, meta_hflx
@@ -54,6 +55,7 @@ MODULE histVars_data
5455
real(dp), allocatable :: instRunoff(:) ! instantaneous lateral inflow into each river/lake [m3/s] [nRch]
5556
real(dp), allocatable :: dlayRunoff(:) ! lateral inflow into river or lake [m3/s] for each reach [nRch]
5657
real(dp), allocatable :: discharge(:,:) ! river/lake discharge [m3/s] for each reach/lake and routing method [nRch,nMethod]
58+
real(dp), allocatable :: inflow(:,:) ! inflow from upstream rivers/lakes [m3/s] for each reach/lake and routing method [nRch,nMethod]
5759
real(dp), allocatable :: elev(:,:) ! river/lake surface water elevation [m] for each reach/lake and routing method [nRch,nMethod]
5860
real(dp), allocatable :: volume(:,:) ! river/lake volume [m3] for each reach/lake and routing method [nRch,nMethod]
5961

@@ -95,31 +97,31 @@ FUNCTION constructor(nHru_local, nRch_local, ierr, message) RESULT(instHistVar)
9597
instHistVar%nRch = nRch_local
9698

9799
if (meta_hflx(ixHFLX%basRunoff)%varFile) then
98-
allocate(instHistVar%basRunoff(nHRU_local), stat=ierr, errmsg=cmessage)
100+
allocate(instHistVar%basRunoff(nHRU_local), source=0._dp, stat=ierr, errmsg=cmessage)
99101
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%basRunoff]'; return; endif
100-
instHistVar%basRunoff(1:nHRU_local) = 0._dp
101102
end if
102103

103104
if (meta_rflx(ixRFLX%instRunoff)%varFile) then
104-
allocate(instHistVar%instRunoff(nRch_local), stat=ierr, errmsg=cmessage)
105+
allocate(instHistVar%instRunoff(nRch_local), source=0._dp, stat=ierr, errmsg=cmessage)
105106
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif
106-
instHistVar%instRunoff(1:nRch_local) = 0._dp
107107
end if
108108

109109
if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then
110-
allocate(instHistVar%dlayRunoff(nRch_local), stat=ierr, errmsg=cmessage)
110+
allocate(instHistVar%dlayRunoff(nRch_local), source=0._dp, stat=ierr, errmsg=cmessage)
111111
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif
112-
instHistVar%dlayRunoff(1:nRch_local) = 0._dp
113112
end if
114113

115114
if (nRoutes>0) then ! this should be number of methods that ouput
116-
allocate(instHistVar%discharge(nRch_local, nRoutes), stat=ierr, errmsg=cmessage)
115+
allocate(instHistVar%discharge(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage)
117116
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%discharge]'; return; endif
118-
instHistVar%discharge(1:nRch_local, 1:nRoutes) = 0._dp
119117

120-
allocate(instHistVar%volume(nRch_local, nRoutes), stat=ierr, errmsg=cmessage)
118+
allocate(instHistVar%volume(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage)
121119
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%volume]'; return; endif
122-
instHistVar%volume(1:nRch_local, 1:nRoutes) = 0._dp
120+
121+
if (outputInflow) then
122+
allocate(instHistVar%inflow(nRch_local, nRoutes), source=0._dp, stat=ierr, errmsg=cmessage)
123+
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%inflow]'; return; endif
124+
end if
123125
end if
124126

125127
END FUNCTION constructor
@@ -205,6 +207,9 @@ SUBROUTINE aggregate(this, & ! inout:
205207
do ix=1,this%nRch
206208
this%discharge(ix,iRoute) = this%discharge(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_Q
207209
this%volume(ix,iRoute) = this%volume(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_VOL(1)
210+
if (outputInflow) then
211+
this%inflow(ix,iRoute) = this%inflow(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_INFLOW
212+
end if
208213
end do
209214
end do
210215

@@ -245,6 +250,11 @@ SUBROUTINE finalize(this)
245250
this%volume = this%volume/real(this%nt, kind=dp)
246251
end if
247252

253+
! 6. inflow
254+
if (allocated(this%inflow)) then
255+
this%inflow = this%inflow/real(this%nt, kind=dp)
256+
end if
257+
248258
END SUBROUTINE finalize
249259

250260
! ---------------------------------------------------------------
@@ -263,6 +273,7 @@ SUBROUTINE refresh(this)
263273
if (allocated(this%dlayRunoff)) this%dlayRunoff = 0._dp
264274
if (allocated(this%discharge)) this%discharge = 0._dp
265275
if (allocated(this%volume)) this%volume = 0._dp
276+
if (allocated(this%inflow)) this%inflow = 0._dp
266277

267278
END SUBROUTINE refresh
268279

@@ -280,6 +291,7 @@ SUBROUTINE clean(this)
280291
if (allocated(this%dlayRunoff)) deallocate(this%dlayRunoff)
281292
if (allocated(this%discharge)) deallocate(this%discharge)
282293
if (allocated(this%volume)) deallocate(this%volume)
294+
if (allocated(this%inflow)) deallocate(this%inflow)
283295

284296
END SUBROUTINE clean
285297

@@ -299,6 +311,7 @@ SUBROUTINE read_restart(this, restart_name, ierr, message)
299311
real(dp), allocatable :: array_tmp(:) ! temp array
300312
integer(i4b) :: ixRoute ! loop index
301313
integer(i4b) :: ixFlow, ixVol ! temporal method index
314+
integer(i4b) :: ixInflow ! temporal method index
302315
logical(lgt) :: FileStatus ! file open or close
303316
type(file_desc_t) :: pioFileDesc ! pio file handle
304317

@@ -369,30 +382,40 @@ SUBROUTINE read_restart(this, restart_name, ierr, message)
369382
allocate(this%discharge(this%nRch, nRoutes), stat=ierr, errmsg=cmessage)
370383
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%discharge]'; return; endif
371384
allocate(this%volume(this%nRch, nRoutes), stat=ierr, errmsg=cmessage)
372-
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%discharge]'; return; endif
385+
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%volume]'; return; endif
373386

374387
this%discharge = 0._dp
375388
this%volume = 0._dp
376389

390+
if (outputInflow) then
391+
allocate(this%inflow(this%nRch, nRoutes), source=0.0_dp, stat=ierr, errmsg=cmessage)
392+
if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [hVars%volume]'; return; endif
393+
end if
394+
377395
do ixRoute=1,nRoutes
378396
select case(routeMethods(ixRoute))
379397
case(accumRunoff)
380398
ixFlow=ixRFLX%sumUpstreamRunoff
381399
case(impulseResponseFunc)
382400
ixFlow=ixRFLX%IRFroutedRunoff
383401
ixVol=ixRFLX%IRFvolume
402+
ixInflow=ixRFLX%IRFinflow
384403
case(kinematicWaveTracking)
385404
ixFlow=ixRFLX%KWTroutedRunoff
386405
ixVol=ixRFLX%KWTvolume
406+
ixInflow=ixRFLX%KWTinflow
387407
case(kinematicWave)
388408
ixFlow=ixRFLX%KWroutedRunoff
389409
ixVol=ixRFLX%KWvolume
410+
ixInflow=ixRFLX%KWinflow
390411
case(muskingumCunge)
391412
ixFlow=ixRFLX%MCroutedRunoff
392413
ixVol=ixRFLX%MCvolume
414+
ixInflow=ixRFLX%MCinflow
393415
case(diffusiveWave)
394416
ixFlow=ixRFLX%DWroutedRunoff
395417
ixVol=ixRFLX%DWvolume
418+
ixInflow=ixRFLX%DWinflow
396419
case default
397420
write(message,'(2A,1X,G0,1X,A)') trim(message), 'routing method index:',routeMethods(ixRoute), 'must be 0-5'
398421
ierr=81; return
@@ -418,13 +441,25 @@ SUBROUTINE read_restart(this, restart_name, ierr, message)
418441
! need to shift tributary part in main core to after halo reaches (nTribOutlet)
419442
if (masterproc) then
420443
this%volume(1:nRch_mainstem, ixRoute) = array_tmp(1:nRch_mainstem)
421-
this%volume(nRch_mainstem+nTribOutlet+1:this%nRch,ixRoute) = this%volume(nRch_mainstem+1:nRch_mainstem+nRch_trib,ixRoute)
444+
this%volume(nRch_mainstem+nTribOutlet+1:this%nRch, ixRoute) = array_tmp(nRch_mainstem+1:nRch_mainstem+nRch_trib)
422445
else
423446
this%volume(:,ixRoute) = array_tmp
424447
end if
425448
end if
426-
end do
427-
end if
449+
450+
if (meta_rflx(ixInflow)%varFile) then
451+
call read_dist_array(pioFileDesc, meta_rflx(ixInflow)%varName, array_tmp, ioDesc_hist_rch_double, ierr, cmessage)
452+
if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
453+
! need to shift tributary part in main core to after halo reaches (nTribOutlet)
454+
if (masterproc) then
455+
this%inflow(1:nRch_mainstem, ixRoute) = array_tmp(1:nRch_mainstem)
456+
this%inflow(nRch_mainstem+nTribOutlet+1:this%nRch, ixRoute) = array_tmp(nRch_mainstem+1:nRch_mainstem+nRch_trib)
457+
else
458+
this%inflow(:,ixRoute) = array_tmp
459+
end if
460+
end if
461+
end do ! end of ixRoute loop
462+
end if ! end of nRoute>0 if-statement
428463

429464
call closeFile(pioFileDesc, FileStatus)
430465

‎route/build/src/historyFile.f90

+10
Original file line numberDiff line numberDiff line change
@@ -462,6 +462,16 @@ SUBROUTINE write_flux_rch(this, hVars_local, ioDescRchFlux, index_write, ierr, m
462462
array_temp(1:nRch_write) = hVars_local%volume(index_write, idxMC)
463463
case(ixRFLX%DWvolume)
464464
array_temp(1:nRch_write) = hVars_local%volume(index_write, idxDW)
465+
case(ixRFLX%KWTinflow)
466+
array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxKWT)
467+
case(ixRFLX%IRFinflow)
468+
array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxIRF)
469+
case(ixRFLX%KWinflow)
470+
array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxKW)
471+
case(ixRFLX%MCinflow)
472+
array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxMC)
473+
case(ixRFLX%DWinflow)
474+
array_temp(1:nRch_write) = hVars_local%inflow(index_write, idxDW)
465475
case default
466476
write(message,'(2A,1X,G0,1X,1A)') trim(message), 'Invalid RFLX variable index:',iVar,'. Check ixRFLX in var_lookup.f90'
467477
ierr=81; return

‎route/build/src/irf_route.f90

+5-2
Original file line numberDiff line numberDiff line change
@@ -85,15 +85,18 @@ SUBROUTINE irf_rch(this, & ! irf_route_rch object to bound this procedur
8585
end if
8686

8787
! get discharge coming from upstream
88-
nUps = size(NETOPO_in(segIndex)%UREACHI)
88+
nUps = count(NETOPO_in(segIndex)%goodBas) ! reminder: goodBas is reach with >0 total contributory area
8989
q_upstream = 0.0_dp
9090
if (nUps>0) then
9191
do iUps = 1,nUps
92-
iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach
92+
if (.not. NETOPO_in(segIndex)%goodBas(iUps)) cycle ! skip upstream reach which does not any flow due to zero total contributory areas
93+
iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach
9394
q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%ROUTE(idxIRF)%REACH_Q
9495
end do
9596
endif
9697

98+
RCHFLX_out(iens,segIndex)%ROUTE(idxIRF)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches
99+
97100
q_upstream_mod = q_upstream
98101
Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1)
99102
Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative)

‎route/build/src/kwe_route.f90

+2
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,8 @@ SUBROUTINE kw_rch(this, & ! kwe_route_rch object to bound this procedu
9292
end do
9393
endif
9494

95+
RCHFLX_out(iens,segIndex)%ROUTE(idxKW)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches
96+
9597
q_upstream_mod = q_upstream
9698
Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1)
9799
Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative)

‎route/build/src/kwt_route.f90

+12
Original file line numberDiff line numberDiff line change
@@ -110,9 +110,12 @@ SUBROUTINE kwt_rch(this, & ! kwt_route_rch object to bound this procedur
110110
! local variables
111111
! (1) extract flow from upstream reaches and append to the non-routed flow in JRCH
112112
integer(i4b) :: NUPS ! number of upstream reaches
113+
integer(i4b) :: iUps ! upstream loop index
114+
integer(i4b) :: iRch_ups ! upstream reach index
113115
real(dp),dimension(:),allocatable :: Q_JRCH ! flow in downstream reach JRCH
114116
real(dp),dimension(:),allocatable :: TENTRY ! entry time to JRCH (exit time u/s)
115117
integer(i4b) :: NQ1 ! # flow particles
118+
real(dp) :: q_upstream ! total discharge at top of the reach [m3/s]
116119
! (2) route flow within the current [JRCH] river segment
117120
integer(i4b) :: ROFFSET ! retrospective offset due to rstep
118121
real(dp) :: T_START ! start of time step
@@ -162,12 +165,21 @@ SUBROUTINE kwt_rch(this, & ! kwt_route_rch object to bound this procedur
162165
ierr=20; message=trim(message)//'negative flow extracted from upstream reach'; return
163166
endif
164167

168+
q_upstream = 0._dp
169+
do iUps = 1,nUps
170+
if (.not. NETOPO_in(segIndex)%goodBas(iUps)) cycle ! skip upstream reach which does not any flow due to zero total contributory areas
171+
iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach
172+
q_upstream = q_upstream + RCHFLX_out(iens,iRch_ups)%ROUTE(idxKWT)%REACH_Q
173+
end do
174+
RCHFLX_out(iens,segIndex)%ROUTE(idxKWT)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches
175+
165176
if(segIndex==ixDesire)then
166177
write(fmt1,'(A,I5,A)') '(A, 1X',size(Q_JRCH),'(1X,G15.4))'
167178
write(iulog,'(a)') ' * Wave discharge from upstream reaches (Q_JRCH) [m2/s]:'
168179
write(iulog,fmt1) ' Q_JRCH=',(Q_JRCH(IWV), IWV=0,size(Q_JRCH)-1)
169180
endif
170181
else
182+
RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_INFLOW = 0._dp
171183
! set flow in headwater reaches to modelled streamflow from time delay histogram
172184
RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q = RCHFLX_out(IENS,segIndex)%BASIN_QR(1)
173185
if (allocated(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE)) then

‎route/build/src/mc_route.f90

+2
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,8 @@ SUBROUTINE mc_rch(this, & ! mc_route_rch object to bound this procedur
9191
end do
9292
endif
9393

94+
RCHFLX_out(iens,segIndex)%ROUTE(idxMC)%REACH_INFLOW = q_upstream ! total inflow from the upstream reaches
95+
9496
q_upstream_mod = q_upstream
9597
Qlat = RCHFLX_out(iens,segIndex)%BASIN_QR(1)
9698
Qabs = RCHFLX_out(iens,segIndex)%REACH_WM_FLUX ! initial water abstraction (positive) or injection (negative)

0 commit comments

Comments
 (0)