Merge pull request #462 from nmizukami/cesm-coupling_add_upstreamflow…

…_output

Add upstream inflow in history output

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]

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)

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,30 +382,40 @@ 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)
               ixFlow=ixRFLX%sumUpstreamRunoff
             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)
 

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

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)

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)

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

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)