diff --git a/route/build/src/ascii_util.f90 b/route/build/src/ascii_util.f90 index b3233cf6..17f30278 100644 --- a/route/build/src/ascii_util.f90 +++ b/route/build/src/ascii_util.f90 @@ -97,7 +97,7 @@ subroutine split_line(inline,words,err,message) ! declare local variables character(len=256) :: temp ! temporary line of characters integer(i4b) :: iword ! loop through words - integer(i4b),parameter :: maxWords=100 ! maximum number of words in a line + integer(i4b),parameter :: maxWords=100 ! maximum number of words in a line integer(i4b) :: i1 ! index at the start of a given word character(len=256) :: cword ! the current word integer(i4b) :: nWords ! number of words in the character string @@ -174,7 +174,7 @@ subroutine get_vlines(unt,vlines,err,message) ! start procedure here err=0; message='get_vlines/' ! ***** get the valid lines of data from the file and store in linked lists ***** - icount=0 ! initialize the counter for the valid lines + icount=0 ! initialize the counter for the valid lines do iline=1,maxLines read(unt,'(a)',iostat=iend)temp; if(iend/=0)exit ! read line of data if (temp(1:1)=='!')cycle diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index da290c54..52e80f4a 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -27,7 +27,7 @@ subroutine upstrm_length(nSeg, & ! input ! ---------------------------------------------------------------------------------------- ! Purpose: ! - ! Calculate total length of upstream reach network from each segment + ! Calculate total length of upstream reach network from each segment ! ! ---------------------------------------------------------------------------------------- ! I/O: @@ -42,8 +42,8 @@ subroutine upstrm_length(nSeg, & ! input ! ! ---------------------------------------------------------------------------------------- USE reachparam - - implicit none + + implicit none ! input variables integer(I4B), intent(in) :: nSeg ! number of stream segments ! output variables @@ -51,14 +51,14 @@ subroutine upstrm_length(nSeg, & ! input character(*), intent(out) :: message ! error message ! local variables integer(I4B) :: iSeg ! index for segment loop - integer(I4B) :: iUps ! index for upstream segment loop - integer(I4B) :: jUps ! index for upstream segment loop + integer(I4B) :: iUps ! index for upstream segment loop + integer(I4B) :: jUps ! index for upstream segment loop INTEGER(I4B) :: NASSIGN ! # reaches currently assigned logical(LGT),dimension(:),allocatable :: RCHFLAG ! TRUE if reach is processed integer(I4B) :: nUps ! number of upstream reaches real(DP) :: xLocal ! length of one segement real(DP) :: xTotal ! total length of upstream segment - + ! initialize error control ierr=0; message='strmlength/' ! ---------------------------------------------------------------------------------------- @@ -67,10 +67,10 @@ subroutine upstrm_length(nSeg, & ! input if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHFLAG'; return; endif RCHFLAG(1:nSeg) = .FALSE. ! ---------------------------------------------------------------------------------------- - + seg_loop: do iSeg=1,nSeg !Loop through each segment - - nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment + + nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment allocate(NETOPO(iSeg)%UPSLENG(nUps),stat=ierr) !print *,'--------------------------------------------' !print *,'Seg ID, Num of upstream', iSeg, nUps @@ -78,17 +78,17 @@ subroutine upstrm_length(nSeg, & ! input upstrms_loop: do iUps=1,nUps !Loop through upstream segments of current segment jUps=NETOPO(iSeg)%RCHLIST(iUps) !index of upstreamf segment xTotal = 0.0 !Initialize total length of upstream segments - do + do xLocal=RPARAM(jUps)%RLENGTH ! Get a segment length - xTotal=xTotal+xLocal - if (jUps.eq.NETOPO(iSeg)%REACHIX) exit - jUps = NETOPO(jUps)%DREACHI ! Get index of immediate downstream segment + xTotal=xTotal+xLocal + if (jUps.eq.NETOPO(iSeg)%REACHIX) exit + jUps = NETOPO(jUps)%DREACHI ! Get index of immediate downstream segment enddo NETOPO(iSeg)%UPSLENG(iUps)=xTotal - enddo upstrms_loop + enddo upstrms_loop enddo seg_loop - - end subroutine upstrm_length + + end subroutine upstrm_length ! ********************************************************************* ! subroutine: compute normalized UH from Saint-Venant Eq. at sim. time step @@ -104,7 +104,7 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! ---------------------------------------------------------------------------------------- ! Purpose: ! - ! Calculate UH at given simulation time step from impulse response function + ! Calculate UH at given simulation time step from impulse response function ! for upstream reach network from each segment ! Use Saint-Venant equation -see equations (13)-(15) in Lohmann, et al. (1996) Tellus ! @@ -113,9 +113,9 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! ! INPUTS: ! nSeg: Number of stream segments in the river network (reaches) - ! dt: Simulation time step [sec] - ! velo: velocity in Saint-Venant eq. [m/sec] - ! diff: diffusivity in Ssaint-Venant eq. [m2/sec] + ! dt: Simulation time step [sec] + ! velo: velocity in Saint-Venant eq. [m/sec] + ! diff: diffusivity in Ssaint-Venant eq. [m2/sec] ! ! ---------------------------------------------------------------------------------------- ! Structures modified: @@ -124,14 +124,14 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! ! ---------------------------------------------------------------------------------------- USE reachparam - USE reach_flux + USE reach_flux ! Define variables implicit none ! input variables integer(i4b), intent(in) :: nSeg ! Numer of segment - real(dp), intent(in) :: dt ! Interval of each time step [sec] + real(dp), intent(in) :: dt ! Interval of each time step [sec] real(dp), intent(in) :: velo ! Wave velocity C for each segment [m/sec] real(dp), intent(in) :: diff ! Diffusivity D for each segment [m2/sec] ! output variables @@ -140,19 +140,19 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! local variable INTEGER(I4B) :: NASSIGN ! # reaches currently assigned logical(LGT),dimension(:),allocatable :: RCHFLAG ! TRUE if reach is processed - real(dp),dimension(:),allocatable :: UHM ! - real(dp),dimension(:),allocatable :: UHQ ! - real(dp),dimension(:),allocatable :: fr ! Unit runoff depth evenly distributed over the simulation duration at hourly step - real(dp) :: POT ! Inside of exponential function of IRF + real(dp),dimension(:),allocatable :: UHM ! + real(dp),dimension(:),allocatable :: UHQ ! + real(dp),dimension(:),allocatable :: fr ! Unit runoff depth evenly distributed over the simulation duration at hourly step + real(dp) :: POT ! Inside of exponential function of IRF real(dp) :: H ! IRF, h(x,t) function h function in eq.15 Lohmann et al.1996 - real(dp) :: UHQ0 ! + real(dp) :: UHQ0 ! real(dp) :: INTE ! Accumulation of hrly UH for ordinate normalization real(dp) :: sec ! hr at each time step [hr] - real(dp),parameter :: dTUH=3600.0 ! UH time step [sec] + real(dp),parameter :: dTUH=3600.0 ! UH time step [sec] integer(I4B) :: nUps ! # of upstream reaches integer(i4b) :: iHrStrt ! index of UH time step where rising limb of UH start integer(i4b) :: iHrLast ! index of UH time step where recession limb of UH become zero - integer(i4b) :: nTSub ! # of time step where 1/nTsub [m] of runoff is inserted + integer(i4b) :: nTSub ! # of time step where 1/nTsub [m] of runoff is inserted integer(i4b) :: iSeg ! Loop index integer(i4b) :: iUps ! Loop index integer(i4b) :: iHr ! Loop index @@ -164,7 +164,7 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! Dynamically assigned parameters nTsub=ceiling(dt/dTUH) !nTsub=floor(dt/dTUH) - + ! initialize error control ierr=0; message='make_uh/' ! ---------------------------------------------------------------------------------------- @@ -184,10 +184,10 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input ! Unit Runoff depth [1/nTsub, 1/nTsub, ..., 1/nTsub, 0, 0, ...0] fr=0._dp - fr(1:nTsub) = 1.0 / nTsub + fr(1:nTsub) = 1.0 / nTsub - seg_loop: do iSeg=1,nSeg - nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment + seg_loop: do iSeg=1,nSeg + nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment allocate(NETOPO(iSeg)%UH(nUps),stat=ierr) if(ierr/=0)then; message=trim(message)//'unable to allocate space for NETOPO%UH'; return; endif @@ -197,24 +197,24 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input sec = 0._dp UHM(:) = 0._dp do iHr=1,nHr - sec = sec + dTUH + sec = sec + dTUH if (velo .GT. 0.0) then POT = ((velo*sec-NETOPO(iSeg)%UPSLENG(iUps))**2.0)/(4.0*diff*sec) - if (POT .GT. 69.0) then + if (POT .GT. 69.0) then H = 0.0 - else + else H = 1.0/(2.0*sqrt(PI*diff*sec))*NETOPO(iSeg)%UPSLENG(iUps)*exp(-POT) - endif - else + endif + else H = 0.0 - endif + endif UHM(iHr) = H - INTE = INTE + H + INTE = INTE + H enddo - !Normalize ordinates by sum of IRF - if (INTE > 0.0) then + !Normalize ordinates by sum of IRF + if (INTE > 0.0) then UHM= UHM/INTE - end if + end if !Get hour indices of start of rising part of UHM and end of recession part of UHM INTE = 0._dp do iHr=1,nTMAX @@ -228,27 +228,27 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input iHrStrt=iHr if (INTE > 0.99999) exit enddo - !Convolute with unit runoff depth, fr + !Convolute with unit runoff depth, fr INTE = 0._dp UHQ(:) = 0._dp - do iHr1 = 1,nTMAX - UHQ0=0._dp + do iHr1 = 1,nTMAX + UHQ0=0._dp do iHr = iHrStrt,iHrLast - if ((iHr1-iHr) > 0) then + if ((iHr1-iHr) > 0) then if (iHr1-iHr <= nTsub ) then UHQ0 = UHQ0 + fr(iHr1-iHr)*UHM(iHr) endif - else - exit - endif - enddo + else + exit + endif + enddo UHQ(iHr1) = UHQ0 INTE=INTE+UHQ0 - end do + end do !Normalize ordinates by sum of total flow - if (INTE > 0.0) then + if (INTE > 0.0) then UHQ= UHQ/INTE - end if + end if !Get hour indices of start of rising part of UH and end of recession part of UH INTE = 0._dp do iHr=1,nTMAX @@ -265,24 +265,24 @@ subroutine make_uh(nSeg,dt,velo,diff, & ! input !Aggregate hourly unit hydrograph to simulation time step allocate(NETOPO(iSeg)%UH(iUps)%UH_DATA((iHrLast+nTsub-1)/nTsub),stat=ierr) if(ierr/=0)then; message=trim(message)//'unable to allocate space for UH%UH_DATA'; return; endif - + NETOPO(iSeg)%UH(iUps)%UH_DATA(:)=0._dp do iHr1 = 1,iHrLast iTagg = (iHr1+nTsub-1)/nTsub NETOPO(iSeg)%UH(iUps)%UH_DATA(iTagg) = NETOPO(iSeg)%UH(iUps)%UH_DATA(iTagg)+UHQ(iHr1) - end do + end do enddo upstrms_loop enddo seg_loop deallocate(UHM) deallocate(UHQ) - + end subroutine make_uh ! ******************************************************************************** - ! subroutine: Sum of basin routed runoff from all the immediate upstream basin - ! + ! subroutine: Sum of basin routed runoff from all the immediate upstream basin + ! ! ********************************************************************************* subroutine get_upsbas_qr(nSeg,iEns, & ! input ierr, message) ! output=error control @@ -295,24 +295,24 @@ subroutine get_upsbas_qr(nSeg,iEns, & ! input ! Purpose: ! ! Find sum of routed runoff volume [m3/s] from immediate upstream basins of the segment (at top of - ! the segment). Segment of Headwater basin does not have flow at top of the segment + ! the segment). Segment of Headwater basin does not have flow at top of the segment ! ! ---------------------------------------------------------------------------------------- ! I/O: ! ! INPUTS: ! nSeg: Number of stream segments in the river network (reaches) - ! iEns: index of ensemble flow + ! iEns: index of ensemble flow ! ! ---------------------------------------------------------------------------------------- ! Structures modified: - ! + ! ! Updates structure RCHFLX%UPSBASIN_QR (module reach_flux) - ! + ! ! ---------------------------------------------------------------------------------------- USE reachparam USE reach_flux - + implicit none ! Input INTEGER(I4B), intent(IN) :: iEns ! ensemble member @@ -320,15 +320,15 @@ subroutine get_upsbas_qr(nSeg,iEns, & ! input ! Output integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message - ! Local variables + ! Local variables INTEGER(I4B) :: NASSIGN ! # reaches currently assigned logical(LGT),dimension(:),allocatable :: RCHFLAG ! TRUE if reach is processed INTEGER(I4B) :: iUpsBas ! loop through u/s reaches INTEGER(I4B) :: jUpsBas ! loop through u/s reaches INTEGER(I4B) :: nUpsBas ! number of immediate upstream basins INTEGER(I4B) :: iSeg ! index of reach with the earliest time - real(DP) :: area ! Area of one upstream basin - real(DP) :: qrTotal ! total routed runoff volume of all the upstream basins + real(DP) :: area ! Area of one upstream basin + real(DP) :: qrTotal ! total routed runoff volume of all the upstream basins ! initialize error control ierr=0; message='get_upsbas_qr/' @@ -347,22 +347,22 @@ subroutine get_upsbas_qr(nSeg,iEns, & ! input if (nUpsBas > 0) then ! if segment is not located in headwater basin do iUpsbas=1,nUpsBas jUpsBas = NETOPO(iSeg)%UREACHI(iUpsBas) ! index of the upstream upstream basin - ! Get routed flow from all immediate upstream basins [m3/s] and sum them up + ! Get routed flow from all immediate upstream basins [m3/s] and sum them up area=RPARAM(jUpsBas)%BASAREA if ( area > 0) then - qrTotal = qrTotal + RCHFLX(iEns,jUpsBas)%BASIN_QR(1)!/area + qrTotal = qrTotal + RCHFLX(iEns,jUpsBas)%BASIN_QR(1)!/area endif ! write(*,'(a,1x,i4,1x,f20.2,1x,f20.4)') 'IupBas, area, routed flow= ',jUpsBas, area, RCHFLX(iEns,jUpsBas)%BASIN_QR(1) - enddo !end of upstream basins + enddo !end of upstream basins endif - RCHFLX(iEns,iSeg)%UPSBASIN_QR = qrTotal + RCHFLX(iEns,iSeg)%UPSBASIN_QR = qrTotal ! write(*,'(a,1x,es14.7)') 'Sum of routed flow = ',RCHFLX(iEns,iSeg)%UPSBASIN_QR enddo seg_loop end subroutine get_upsbas_qr ! ********************************************************************* - ! subroutine: Compute delayed runoff from all the upstream segments + ! subroutine: Compute delayed runoff from all the upstream segments ! ********************************************************************* subroutine conv_upsbas_qr(nSeg,iEns, & ! input ierr, message) ! output=error control @@ -374,24 +374,24 @@ subroutine conv_upsbas_qr(nSeg,iEns, & ! input ! ---------------------------------------------------------------------------------------- ! Purpose: ! - ! Convolute routed basisn flow volume at top of each of the upstream segment at one time step and at each segment - ! + ! Convolute routed basisn flow volume at top of each of the upstream segment at one time step and at each segment + ! ! ---------------------------------------------------------------------------------------- ! I/O: ! ! INPUTS: ! nSeg: Number of stream segments in the river network (reaches) - ! iEns: index of ensemble flow + ! iEns: index of ensemble flow ! ! ---------------------------------------------------------------------------------------- ! Structures modified: - ! + ! ! Updates structure RCHFLX%BASIN_QR_IRF (module reach_flux) ! RCHFLX%REACH_Q_IRF ! ---------------------------------------------------------------------------------------- USE reachparam USE reach_flux - + implicit none ! Input INTEGER(I4B), intent(IN) :: iens ! ensemble member @@ -399,14 +399,14 @@ subroutine conv_upsbas_qr(nSeg,iEns, & ! input ! Output integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message - ! Local variables to + ! Local variables to INTEGER(I4B) :: NASSIGN ! # reaches currently assigned logical(LGT),dimension(:),allocatable :: RCHFLAG ! TRUE if reach is processed INTEGER(I4B) :: iUps ! loop through u/s reaches - INTEGER(I4B) :: jUps ! - INTEGER(I4B) :: nUps ! number of all upstream segment + INTEGER(I4B) :: jUps ! + INTEGER(I4B) :: nUps ! number of all upstream segment INTEGER(I4B) :: iSeg ! index of segment - INTEGER(I4B) :: nTDH ! number of UH data + INTEGER(I4B) :: nTDH ! number of UH data INTEGER(I4B) :: iTDH ! index of UH data (i.e.,future time step) ! initialize error control @@ -418,16 +418,16 @@ subroutine conv_upsbas_qr(nSeg,iEns, & ! input RCHFLAG(1:nSeg) = .FALSE. ! ---------------------------------------------------------------------------------------- - ! route one time step routed runoff volue from each up stream basins to segment + ! route one time step routed runoff volue from each up stream basins to segment seg_loop: do iSeg=1,nSeg - nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment + nUps = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment !print *,'--------------------------' !print *,'Working on iSeg, ID, # of ups. basins= ', iSeg, NETOPO(iSeg)%RCHLIST, nUps !print *,'upbasinflow= ',RCHFLX(iens,iSeg)%UPSBASIN_QR upstrms_loop: do iUps=1,nUps ! index of the upstream upstream basin - jUps = NETOPO(iSeg)%RCHLIST(iUps) + jUps = NETOPO(iSeg)%RCHLIST(iUps) ! identify the number of future time steps for a given segment nTDH = size(NETOPO(iSeg)%UH(iUps)%UH_DATA) @@ -436,8 +436,8 @@ subroutine conv_upsbas_qr(nSeg,iEns, & ! input UH_data: do iTDH=1,nTDH RCHFLX(iens,iSeg)%QFUTURE_IRF(iTDH) = RCHFLX(iens,iSeg)%QFUTURE_IRF(iTDH) & + NETOPO(iSeg)%UH(iUps)%UH_DATA(iTDH)*RCHFLX(iEns,jUps)%UPSBASIN_QR - enddo UH_data - enddo upstrms_loop + enddo UH_data + enddo upstrms_loop ! save the routed runoff RCHFLX(iEns,iSeg)%BASIN_QR_IRF(0) = RCHFLX(iEns,iSeg)%BASIN_QR_IRF(1) ! (save the runoff from the previous time step) diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index 70bc6f63..3d44ec0a 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -83,11 +83,11 @@ SUBROUTINE REACHORDER(NRCH, & ! input CYCLE ENDIF ! climb upstream as far as possible - JRCH = IRCH ! the first reach under investigation + JRCH = IRCH ! the first reach under investigation DO ! do until get to a "most upstream" reach that is not assigned NUPS = SIZE(NETOPO(JRCH)%UREACHI) ! get number of upstream reaches IF (NUPS.GE.1) THEN ! (if NUPS = 0, then it is a first-order basin) - KRCH = JRCH ! the reach under investigation + KRCH = JRCH ! the reach under investigation ! loop through upstream reaches DO IUPS=1,NUPS UINDEX = NETOPO(JRCH)%UREACHI(IUPS) ! POSITION of the upstream reach @@ -104,7 +104,7 @@ SUBROUTINE REACHORDER(NRCH, & ! input RCHFLAG(JRCH) = .TRUE. NETOPO(ICOUNT)%RHORDER = JRCH EXIT - ENDIF + ENDIF CYCLE ! if jrch changes, keep looping (move upstream) ELSE ! if the reach is a first-order basin ! assign JRCH @@ -199,7 +199,7 @@ SUBROUTINE REACH_LIST(NRCH,NTOTAL,ierr,message) INTLIST(IRCH)%N_URCH = 0 ! initialize the number of upstream reaches NULLIFY(INTLIST(IRCH)%HPOINT) ! set pointer to a linked list to NULL END DO ! (irch) - + ! build the linked lists for all reaches DO KRCH=1,NRCH ! ensure take streamflow from surrounding basin (a reach is upstream of itself!) @@ -235,14 +235,14 @@ SUBROUTINE REACH_LIST(NRCH,NTOTAL,ierr,message) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif print*, 'jrch, numups, NETOPO(JRCH)%RCHLIST(:) = ', jrch, numups, NETOPO(JRCH)%RCHLIST(:) END DO ! jrch - + ! free up memory DEALLOCATE(INTLIST,STAT=IERR) if(ierr/=0)then; ierr=20; message=trim(message)//'problem deallocating space for INTLIST'; return; endif ! ---------------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------------- CONTAINS - + ! For a down stream reach, add an upstream reach to its list of upstream reaches SUBROUTINE ADD2LIST(D_RCH,U_RCH,ierr,message) INTEGER(I4B),INTENT(IN) :: U_RCH ! upstream reach index @@ -324,7 +324,7 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices ! JRCH: index of stream segment ! T0: start of the time step (seconds) ! T1: end of the time step (seconds) - ! LAKEFLAG: >0 if processing lakes + ! LAKEFLAG: >0 if processing lakes ! RSTEP: retrospective time step offset (optional) ! ! Outputs (in addition to update of data structures): @@ -388,7 +388,7 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices ! ! Most computations were originally performed within calcts in Topnet ver7, with calls ! to subroutines in kinwav_v7.f - ! + ! ! ---------------------------------------------------------------------------------------- ! Modifications to Source (mclark@ucar.edu): ! @@ -401,7 +401,7 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices ! * use of a new data structure (KROUTE) to hold and update the flow particles ! ! * upgrade to F90 (especially structured variables and dynamic memory allocation) - ! + ! ! ---------------------------------------------------------------------------------------- ! Future revisions: ! @@ -475,11 +475,11 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices ! check if(JRCH==ixPrint)then print*, 'JRCH, Q_JRCH = ', JRCH, Q_JRCH - endif + endif ELSE ! set flow in headwater reaches to modelled streamflow from time delay histogram - RCHFLX(IENS,JRCH)%REACH_Q = RCHFLX(IENS,JRCH)%BASIN_QR(1) + RCHFLX(IENS,JRCH)%REACH_Q = RCHFLX(IENS,JRCH)%BASIN_QR(1) RETURN ! no upstream reaches (routing for sub-basins done using time-delay histogram) ENDIF ! ---------------------------------------------------------------------------------------- @@ -515,7 +515,7 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices print*, 'FROUTE = ', FROUTE print*, 'TENTRY = ', TENTRY print*, 'T_EXIT = ', T_EXIT - endif + endif ! ---------------------------------------------------------------------------------------- ! (4) COMPUTE TIME-STEP AVERAGES @@ -523,7 +523,7 @@ SUBROUTINE QROUTE_RCH(IENS,JRCH, & ! input: array indices NR = COUNT(FROUTE)-1 ! -1 because of the zero element (last routed) NN = NQ2-NR ! number of non-routed points TNEW = (/T_START,T_END/) - ! (zero position last routed; use of NR+1 instead of NR keeps next expected routed point) + ! (zero position last routed; use of NR+1 instead of NR keeps next expected routed point) CALL INTERP_RCH(T_EXIT(0:NR+1),Q_JRCH(0:NR+1),TNEW,QNEW,IERR,CMESSAGE) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif if(JRCH == ixPrint) print*, 'QNEW(1) = ', QNEW(1) @@ -695,14 +695,14 @@ SUBROUTINE GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,& ! input ROFFSET = RSTEP END IF IF (LAKEFLAG.EQ.1) THEN ! lakes are enabled - ! get lake outflow and only lake outflow if reach is a lake outlet reach, else do as normal + ! get lake outflow and only lake outflow if reach is a lake outlet reach, else do as normal ILAK = NETOPO(JRCH)%LAKE_IX ! lake index IF (ILAK.GT.0) THEN ! part of reach is in lake IF (NETOPO(JRCH)%REACHIX.EQ.LKTOPO(ILAK)%DREACHI) THEN ! we are in a lake outlet reach ND = 1 ALLOCATE(QD(1),TD(1),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating array for QD and TD'; return; endif - QD(1) = LAKFLX(IENS,ILAK)%LAKE_Q / RPARAM(JRCH)%R_WIDTH ! lake outflow per unit reach width + QD(1) = LAKFLX(IENS,ILAK)%LAKE_Q / RPARAM(JRCH)%R_WIDTH ! lake outflow per unit reach width TD(1) = T1 - DT*ROFFSET ELSE CALL QEXMUL_RCH(IENS,JRCH,T0,T1,ND,QD,TD,ierr,cmessage,RSTEP) ! do as normal for unsubmerged part of inlet reach @@ -827,7 +827,7 @@ SUBROUTINE QEXMUL_RCH(IENS,JRCH,T0,T1,& ! input INTEGER(I4B) :: IR ! index of the upstream reach INTEGER(I4B) :: NS ! size of the wave INTEGER(I4B) :: NR ! # routed particles in u/s reach - INTEGER(I4B) :: NQ ! NR+1, if non-routed particle exists + INTEGER(I4B) :: NQ ! NR+1, if non-routed particle exists TYPE(FPOINT), DIMENSION(:), POINTER, SAVE :: NEW_WAVE ! temporary wave LOGICAL(LGT),SAVE :: INIT=.TRUE. ! used to initialize pointers ! Local variables to merge flow @@ -987,7 +987,7 @@ SUBROUTINE QEXMUL_RCH(IENS,JRCH,T0,T1,& ! input if(ierr/=0)then; message=trim(message)//'problem deallocating array NEW_WAVE'; return; endif NULLIFY(NEW_WAVE) ! save the upstream width - UWIDTH(NUPB+IUPR) = RPARAM(IR)%R_WIDTH ! reach, width = parameter + UWIDTH(NUPB+IUPR) = RPARAM(IR)%R_WIDTH ! reach, width = parameter ! save the time for the first particle in each reach CTIME(NUPB+IUPR) = USFLOW(NUPB+IUPR)%KWAVE(1)%TR ! central time ! keep track of the total number of points that must be routed downstream @@ -1001,7 +1001,7 @@ SUBROUTINE QEXMUL_RCH(IENS,JRCH,T0,T1,& ! input ! *other than* x, we need to estimate (interpolate) flow for the *times* associted with ! each of the flow particles in reach x. Then, at a given time, we can sum the flow ! (routed in reach x plus interpolated flow in all other reaches). This needs to be done - ! for all upstream reaches. + ! for all upstream reaches. ! ---------------------------------------------------------------------------------------- ! We accomplish this as follows. We define a vector of indices (ITIM), where each ! element of ITIM points to a particle in a given upstream reach still to be processed. @@ -1012,7 +1012,7 @@ SUBROUTINE QEXMUL_RCH(IENS,JRCH,T0,T1,& ! input ! reaches by the width of the downstream reach, and sum the flow over all upstream reaches. ! We then move the index forward in ITIM (for the upstream reach just processed), get a ! new vector CTIME, and process the next earliest particle. We continue until all - ! flow particles are processed in all upstream reaches. + ! flow particles are processed in all upstream reaches. ! ---------------------------------------------------------------------------------------- IPRT = 0 ! initialize counter for flow particles in the output array ! allocate space for the merged flow at the downstream reach @@ -1181,7 +1181,7 @@ SUBROUTINE REMOVE_RCH(MAXQPAR,& ! input INTEGER(I4B) :: IPRT ! loop through flow particles REAL(DP), DIMENSION(:), ALLOCATABLE :: Q,T,Z ! copies of Q_JRCH and T_JRCH LOGICAL(LGT), DIMENSION(:), ALLOCATABLE :: PARFLG ! .FALSE. if particle removed - INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: INDEX0 ! indices of original vectors + INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: INDEX0 ! indices of original vectors REAL(DP), DIMENSION(:), ALLOCATABLE :: ABSERR ! absolute error btw interp and orig REAL(DP) :: Q_INTP ! interpolated particle INTEGER(I4B) :: MPRT ! local number of flow particles @@ -1265,8 +1265,8 @@ SUBROUTINE REMOVE_RCH(MAXQPAR,& ! input ! ---------------------------------------------------------------------------------------- CONTAINS FUNCTION INTERP(T0,Q1,Q2,T1,T2) - REAL(DP),INTENT(IN) :: Q1,Q2 ! flow at neighbouring times - REAL(DP),INTENT(IN) :: T1,T2 ! neighbouring times + REAL(DP),INTENT(IN) :: Q1,Q2 ! flow at neighbouring times + REAL(DP),INTENT(IN) :: T1,T2 ! neighbouring times REAL(DP),INTENT(IN) :: T0 ! desired time REAL(DP) :: INTERP ! function name ! dQ/dT dT @@ -1315,7 +1315,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& ! flow either side of a shock -- thus we may have ! fewer elements on output if several particles are ! merged, INTENT(INOUT) - ! TENTRY: array of time elements -- neighbouring times are merged if a shock forms, + ! TENTRY: array of time elements -- neighbouring times are merged if a shock forms, ! then merged times are dis-aggregated, one second is ! added to the time corresponding to the higer merged ! flow (note also fewer elements), INTENT(INOUT) @@ -1351,7 +1351,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& ! ! ---------------------------------------------------------------------------------------- ! Source: - ! + ! ! This routine is based on the subroutine kinwav, located in kinwav_v7.f ! ! ---------------------------------------------------------------------------------------- @@ -1408,7 +1408,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& REAL(DP) :: XXB ! wave break INTEGER(I4B) :: IXB,JXB ! define position of wave break REAL(DP) :: A1,A2 ! stage - different sides of break - REAL(DP) :: CM ! merged celerity + REAL(DP) :: CM ! merged celerity REAL(DP) :: TEXIT ! expected exit time of "current" particle REAL(DP) :: TNEXT ! expected exit time of "next" particle REAL(DP) :: TEXIT2 ! exit time of "bottom" of merged element @@ -1417,7 +1417,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& INTEGER(I4B) :: ICOUNT ! used to account for merged pts character(len=256) :: cmessage ! error message of downwind routine ! ---------------------------------------------------------------------------------------- - ! NOTE: If merged particles DO NOT exit the reach in the current time step, they are + ! NOTE: If merged particles DO NOT exit the reach in the current time step, they are ! disaggregated into the original particles; if the merged particles DO exit the ! reach, then we save only the "slowest" and "fastest" particle. ! ---------------------------------------------------------------------------------------- @@ -1456,7 +1456,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& WC(1:NN) = ALFA*K**(1./ALFA)*Q1(1:NN)**((ALFA-1.)/ALFA) GT_ONE: IF(NN.GT.1) THEN ! no breaking if just one point X = 0. ! altered later to describe "closest" shock - GOTALL: DO ! keep going until all shocks are merged + GOTALL: DO ! keep going until all shocks are merged XB = XMX ! initialized to length of the stream segment ! -------------------------------------------------------------------------------------- ! check for breaking @@ -1496,7 +1496,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& IX(IXB:NN) = IX(IXB+1:NN+1) ! index (minimum index value of each merged particle) T1(IXB:NN) = T1(IXB+1:NN+1) ! entry time WC(IXB:NN) = WC(IXB+1:NN+1) ! wave celerity - Q1(IXB:NN) = Q1(IXB+1:NN+1) ! unmodified flows + Q1(IXB:NN) = Q1(IXB+1:NN+1) ! unmodified flows Q2(IXB:NN) = Q2(IXB+1:NN+1) ! unmodified flows ! update X - already got the "closest shock to start", see if there are any other shocks X = XB @@ -1536,7 +1536,7 @@ SUBROUTINE KINWAV_RCH(JRCH,T_START,T_END,Q_JRCH,TENTRY,FROUTE,T_EXIT,NQ2,& if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ELSE ! merged elements have not exited ! when a merged element does not exit, need to disaggregate into original particles - DO JROUTE=1,NI ! loop thru # original inputs + DO JROUTE=1,NI ! loop thru # original inputs IF(MF(JROUTE).EQ.IROUTE) & CALL RUPDATE(Q0(JROUTE),T0(JROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q0, T0, + run checks if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif diff --git a/route/build/src/lake_param.f90 b/route/build/src/lake_param.f90 index 10ba4e68..314397bc 100644 --- a/route/build/src/lake_param.f90 +++ b/route/build/src/lake_param.f90 @@ -15,7 +15,7 @@ MODULE lake_param REAL(DP) :: DSCHECO ! discharge at "ecological" height REAL(DP) :: DSCHSPL ! discharge at spillway height REAL(DP) :: RATECVA ! discharge rating curve parameter - REAL(DP) :: RATECVB ! discharge rating curve parameter + REAL(DP) :: RATECVB ! discharge rating curve parameter ENDTYPE LAKPRP TYPE(LAKPRP), DIMENSION(:), POINTER :: LPARAM ! Reach Parameters ! Lake topology diff --git a/route/build/src/qtimedelay.f90 b/route/build/src/qtimedelay.f90 index c5ad4ef3..9f51256c 100644 --- a/route/build/src/qtimedelay.f90 +++ b/route/build/src/qtimedelay.f90 @@ -20,7 +20,7 @@ SUBROUTINE QTIMEDELAY(dt, fshape, tscale, IERR, MESSAGE) ! input REAL(DP), INTENT(IN) :: dt ! model time step REAL(SP), INTENT(IN) :: fshape ! shapef parameter in gamma distribution -REAL(DP), INTENT(IN) :: tscale ! time scale parameter +REAL(DP), INTENT(IN) :: tscale ! time scale parameter ! output INTEGER(I4B), INTENT(OUT) :: IERR ! error code CHARACTER(*), INTENT(OUT) :: MESSAGE ! error message diff --git a/route/build/src/read_ntopo.f90 b/route/build/src/read_ntopo.f90 index 6096f973..7113975d 100644 --- a/route/build/src/read_ntopo.f90 +++ b/route/build/src/read_ntopo.f90 @@ -13,11 +13,11 @@ module read_ntopo contains ! ********************************************************************* -! subroutine: get vector dimension from netCDF +! subroutine: get vector dimension from netCDF ! ********************************************************************* subroutine get_vec_dim(fname, & ! input: filename dname, & ! input: variable name - nDim, & ! output: Size of dimension + nDim, & ! output: Size of dimension ierr, message) ! output: error control implicit none ! input variables @@ -81,7 +81,7 @@ subroutine get_vec_ivar(fname, & ! input: filename ierr = nf90_open(trim(fname),nf90_nowrite,ncid) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! allocate space for the output + ! allocate space for the output !allocate(iVec(iCount),stat=ierr) !if(ierr/=0)then; message=trim(message)//'problem allocating space for iVec'; return; endif @@ -222,7 +222,7 @@ subroutine get_vec_dvar(fname, & ! input: filename ierr = nf90_open(trim(fname),nf90_nowrite,ncid) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! allocate space for the output + ! allocate space for the output ! allocate(dVec(iCount),stat=ierr) ! if(ierr/=0)then; message=trim(message)//'problem allocating space for dVec'; return; endif diff --git a/route/build/src/read_simoutput.f90 b/route/build/src/read_simoutput.f90 index 8ab630d2..62b35195 100644 --- a/route/build/src/read_simoutput.f90 +++ b/route/build/src/read_simoutput.f90 @@ -52,7 +52,7 @@ subroutine get_qDims(fname, & ! input: filename ! get the ID of the time variable ierr = nf90_inq_varid(ncid, trim(vname_time), ivarID) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - + ! get the time units ierr = nf90_get_att(ncid, ivarID, 'units', units_time) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif diff --git a/route/build/src/read_streamSeg.f90 b/route/build/src/read_streamSeg.f90 index 0a0b2f2b..a2f25feb 100644 --- a/route/build/src/read_streamSeg.f90 +++ b/route/build/src/read_streamSeg.f90 @@ -19,7 +19,7 @@ subroutine getData(fname, & ! input: file name imap_acil, & ! input-output: ancillary data for mapping hru2basin ntop_acil, & ! input-output: ancillary data for network topology nHRU, & ! output: number of HRUs - nSeg, & ! output: number of stream segments + nSeg, & ! output: number of stream segments ierr, message) ! output: error control USE dataTypes,only:namepvar,nameivar ! provide access to data types implicit none @@ -216,7 +216,7 @@ subroutine hru2basin(nHRU, & ! input: number of HRUs ! compute total area of the HRUs draining to the stream segment totarea = sum(h2b(iSeg)%cHRU(:)%hru_area) - + ! compute the weights h2b(iSeg)%cHRU(:)%wght = h2b(iSeg)%cHRU(:)%hru_area / totarea @@ -245,7 +245,7 @@ subroutine assign_reachparam(nRch, & ! input: number of stream segmen implicit none ! input variables integer(i4b), intent(in) :: nRch ! number of reaches - type(namepvar), intent(in) :: sseg_acil(:) ! ancillary data for stream segments + type(namepvar), intent(in) :: sseg_acil(:) ! ancillary data for stream segments type(nameivar), intent(in) :: ntop_acil(:) ! ancillary data for the network topology ! output variables integer(i4b), intent(out) :: ierr ! error code @@ -271,12 +271,12 @@ subroutine assign_reachparam(nRch, & ! input: number of stream segmen RPARAM(:)%UPSAREA = sseg_acil(ixSEG%upArea )%varData(:) RPARAM(:)%RLENGTH = sseg_acil(ixSEG%length )%varData(:) RPARAM(:)%R_SLOPE = sseg_acil(ixSEG%slope )%varData(:) - + ! compute area draining to each stream segment do iRch=1,nRch RPARAM(iRch)%BASAREA = sum(h2b(iRch)%cHRU(:)%hru_area) end do - + ! compute total area above the bottom of each reach (m2) RPARAM(:)%TOTAREA = RPARAM(:)%UPSAREA + RPARAM(:)%BASAREA @@ -300,7 +300,7 @@ subroutine assign_reachparam(nRch, & ! input: number of stream segmen NETOPO(:)%RCHLAT2 = huge(kind(dp)) ! End latitude NETOPO(:)%RCHLON1 = huge(kind(dp)) ! Start longitude NETOPO(:)%RCHLON2 = huge(kind(dp)) ! End longitude - NETOPO(:)%LAKE_IX = -1 ! Lake index (0,1,2,...,nlak-1) + NETOPO(:)%LAKE_IX = -1 ! Lake index (0,1,2,...,nlak-1) NETOPO(:)%LAKE_ID = -1 ! Lake ID (REC code?) NETOPO(:)%BASULAK = 0._dp ! Area of basin under lake NETOPO(:)%RCHULAK = 0._dp ! Length of reach under lake @@ -335,7 +335,7 @@ subroutine define_topology(nRch, & ! input ! define reach indices NETOPO(:)%REACHIX = arth(1,1,nRch) - NETOPO(:)%DREACHI = -1 + NETOPO(:)%DREACHI = -1 hFlag(:) = .false. ! first identify non-headwater basins with no upstream area diff --git a/route/build/src/route_runoff.f90 b/route/build/src/route_runoff.f90 index 9cdf93f0..820ab0d3 100644 --- a/route/build/src/route_runoff.f90 +++ b/route/build/src/route_runoff.f90 @@ -9,7 +9,7 @@ program route_runoff USE var_lookup,only:ixSEG,nVarsSEG ! index of variables for the stream segments USE var_lookup,only:ixMAP,nVarsMAP ! index of variables for the hru2basin mapping USE var_lookup,only:ixTOP,nVarsTOP ! index of variables for the network topology -USE reachparam ! reach parameters +USE reachparam ! reach parameters USE reachstate ! reach states USE reach_flux ! fluxes in each reach USE nhru2basin ! data structures holding the nhru2basin correspondence @@ -45,7 +45,7 @@ program route_runoff real(dp),parameter :: secprday=86400._dp ! number of seconds in a day logical(lgt),parameter :: doIRFroute=.True. ! switch to turn off/on IRF routing added by NM logical(lgt),parameter :: doKWTroute=.True. -real(dp),parameter :: verySmall=tiny(1.0_dp) ! a very small number +real(dp),parameter :: verySmall=tiny(1.0_dp) ! a very small number ! general guff integer(i4b),parameter :: strLen=256 ! length of character string integer(i4b) :: ierr ! error code @@ -54,7 +54,7 @@ program route_runoff character(len=strLen) :: str ! miscellaneous string ! read control file character(len=strLen) :: cfile_name ! name of the control file -character(len=strLen) :: param_nml ! name of the namelist file +character(len=strLen) :: param_nml ! name of the namelist file integer(i4b) :: iunit ! file unit character(len=strLen),allocatable :: cLines(:) ! vector of character strings integer(i4b) :: ipos ! index of character string @@ -84,20 +84,20 @@ program route_runoff integer(i4b) :: iStart ! start index of the ragged array integer(i4b),dimension(1) :: iDesire ! index of stream segment with maximum upstream area (vector) integer(i4b) :: ixDesire ! index of stream segment with maximum upstream area (scalar) -!integer(i4b) :: iTDH ! index of unit hydrograph data element +!integer(i4b) :: iTDH ! index of unit hydrograph data element ! define stream network information character(len=strLen) :: fname_ntop ! filename containing stream network topology information integer(i4b),allocatable :: REACHIDGV(:) integer(i4b),allocatable :: RCHIXLIST(:) integer(i4b) :: nTotal ! total number of upstream segments for all stream segments -integer(i4b) :: iRchStart -integer(i4b) :: iRchStart1 -integer(i4b),target :: nRchCount -integer(i4b) :: nRchCount1 -integer(i4b) :: iUpRchStart -integer(i4b) :: nUpRchCount -integer(i4b) :: iUpHruStart -integer(i4b) :: nUpHruCount +integer(i4b) :: iRchStart +integer(i4b) :: iRchStart1 +integer(i4b),target :: nRchCount +integer(i4b) :: nRchCount1 +integer(i4b) :: iUpRchStart +integer(i4b) :: nUpRchCount +integer(i4b) :: iUpHruStart +integer(i4b) :: nUpHruCount integer(i4b),allocatable :: upStrmRchList(:) ! define metadata from model output file character(len=strLen) :: fname_qsim ! filename containing simulated runoff @@ -136,7 +136,7 @@ program route_runoff ! route delaied runoff through river network with St.Venant UH real(dp) :: velo ! velocity [m/s] for Saint-Venant equation added by NM real(dp) :: diff ! diffusivity [m2/s] for Saint-Venant equation added by NM -integer(i4b) :: nUH_DATA_MAX ! maximum number of elements in the UH data among all the upstreamfs for a segment +integer(i4b) :: nUH_DATA_MAX ! maximum number of elements in the UH data among all the upstreamfs for a segment ! compute total instantaneous runoff upstream of each reach integer(i4b),allocatable :: iUpstream(:) ! indices for all reaches upstream real(dp),allocatable :: qUpstream(:) ! streamflow for all reaches upstream @@ -193,9 +193,9 @@ program route_runoff case(''); output_dir = trim(cData) ! directory containing output data ! define variables for the network topology case(''); fname_ntop = trim(cData) ! name of file containing stream network topology information - case('' ); read(cData,*,iostat=ierr) iSegOut ! seg_id of outlet streamflow segment + case('' ); read(cData,*,iostat=ierr) iSegOut ! seg_id of outlet streamflow segment if(ierr/=0) call handle_err(70,'problem with internal read of iSegOut info, read from control file') - ! define the file and variable name for runoff netCDF + ! define the file and variable name for runoff netCDF case(''); fname_qsim = trim(cData) ! name of file containing the runoff case(''); vname_qsim = trim(cData) ! name of runoff variable case(''); vname_hruid = trim(cData) ! name of the HRU id @@ -206,7 +206,7 @@ program route_runoff ! define the output filename case(''); fname_output = trim(cData) ! filename for the model output ! define namelist name for routing parameters - case(''); param_nml = trim(cData) ! name of namelist including routing parameter value + case(''); param_nml = trim(cData) ! name of namelist including routing parameter value end select end do ! looping through lines in the control file @@ -250,7 +250,7 @@ program route_runoff ! ***** ! (1.2) Read in the network topology information... ! ********************************************* -! Data type to be populated +! Data type to be populated ! NETOPO(:)%REACHIX index ! NETOPO(:)%REACHID id ! NETOPO(:)%DREACHK id @@ -268,12 +268,12 @@ program route_runoff ! NETOPO(:)%UREACHK(:) id ! NETOPO(:)%UREACHI(:) index ! h2b(:)%cHRU(:)%hru_ix index -! h2b(:)%cHRU(:)%hru_id id -! h2b(:)%cHRU(:)%hru_lon -! h2b(:)%cHRU(:)%hru_lat +! h2b(:)%cHRU(:)%hru_id id +! h2b(:)%cHRU(:)%hru_lon +! h2b(:)%cHRU(:)%hru_lat ! h2b(:)%cHRU(:)%hru_elev ! h2b(:)%cHRU(:)%hru_area -! h2b(:)%cHRU(:)%wght +! h2b(:)%cHRU(:)%wght ! Read global reach id (input = filename, variable name, variable vector, start index; output = error control) allocate(REACHIDGV(nSeg), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for REACHIDGV') @@ -292,19 +292,19 @@ program route_runoff call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachCount', nRchCount, iSegDesire, ierr, cmessage); call handle_err(ierr,cmessage) print*,'iRchStart = ',iRchStart print*,'Number of upstream segment from outlet segment (nRchCount): ',nRchCount - - ! Read reach list of index from global segments (all the upstream reachs for each segment) + + ! Read reach list of index from global segments (all the upstream reachs for each segment) allocate(upStrmRchList(nRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for upStrmRchList') call get_vec_ivar(trim(ancil_dir)//trim(fname_ntop), 'reachList', upStrmRchList, iRchStart,nRchCount, ierr, cmessage); call handle_err(ierr,cmessage) - - ! Reach upstream segment and associated HRU infor from non-ragged vector + + ! Reach upstream segment and associated HRU infor from non-ragged vector allocate(NETOPO(nRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO') allocate(RPARAM(nRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for RPARAM') allocate(h2b(nRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for h2b') - + ! Create REACH index for local segments NETOPO(:)%REACHIX=arth(1,1,nRchCount) - + do iSeg=1,nRchCount ! Reach reach topology and parameters (integer) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachID', NETOPO(iSeg)%REACHID,upStrmRchList(iSeg),ierr,cmessage); call handle_err(ierr,cmessage) @@ -320,39 +320,39 @@ program route_runoff call get_scl_dvar(trim(ancil_dir)//trim(fname_ntop),'reachLon1', NETOPO(iSeg)%RCHLON1,upStrmRchList(iSeg),ierr,cmessage); call handle_err(ierr,cmessage) call get_scl_dvar(trim(ancil_dir)//trim(fname_ntop),'reachLon2', NETOPO(iSeg)%RCHLON2,upStrmRchList(iSeg),ierr,cmessage); call handle_err(ierr,cmessage) enddo - + ! Recompute downstream segment index as local segment list, or NETOPO(:)%REACHID do iSeg=1,nRchCount ! Assign downstream segment ID = 0 at desired outlet segment - if (NETOPO(iSeg)%REACHID == iSegOut) then + if (NETOPO(iSeg)%REACHID == iSegOut) then NETOPO(iSeg)%DREACHK = 0 else ! Identify the index of the desired stream segment from reachID vector (dimension size: nSeg) iSelect = minloc(abs(NETOPO(:)%REACHID - NETOPO(iSeg)%DREACHK)) NETOPO(iSeg)%DREACHI = iSelect(1) ! de-vectorize the desired stream segment - if (NETOPO(NETOPO(iSeg)%DREACHI)%REACHID /= NETOPO(iSeg)%DREACHK) then - print*,'iSeg = ', iSeg + if (NETOPO(NETOPO(iSeg)%DREACHI)%REACHID /= NETOPO(iSeg)%DREACHK) then + print*,'iSeg = ', iSeg print*,'NETOPO(iSeg)%DREACHK = ', NETOPO(iSeg)%DREACHK print*,'NETOPO(NETOPO(iSeg)%DREACHI)%REACHID = ', NETOPO(NETOPO(iSeg)%DREACHI)%REACHID call handle_err(20,'unable to find desired downstream segment') endif endif enddo - - ! Reach upstream segment and associated HRU infor from ragged vector + + ! Reach upstream segment and associated HRU infor from ragged vector nTotal=0 do iSeg=1,nRchCount - ! sAll dimension + ! sAll dimension call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachStart', iRchStart1, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachCount', nRchCount1, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) - + allocate(NETOPO(iSeg)%UPSLENG(nRchCount1), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%UPSLENG') allocate(NETOPO(iSeg)%RCHLIST(nRchCount1), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%RCHLIST') allocate(RCHIXLIST(nRchCount1), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for RCHIXLIST(nRchCount)') - + call get_vec_ivar(trim(ancil_dir)//trim(fname_ntop),'reachList' ,RCHIXLIST,iRchStart1,nRchCount1,ierr,cmessage); call handle_err(ierr,cmessage) call get_vec_dvar(trim(ancil_dir)//trim(fname_ntop),'upReachTotalLength',NETOPO(iSeg)%UPSLENG(:),iRchStart1,nRchCount1,ierr,cmessage); call handle_err(ierr,cmessage) - + ! Recompute all the upstream segment indices as local segment list = NETOPO(:)%REACHID nTotal = nTotal + nRchCount1 do jSeg=1,nRchCount1 @@ -362,43 +362,43 @@ program route_runoff enddo !print*,'NETOPO(iSeg)%RCHLIST(:) = ',NETOPO(iSeg)%RCHLIST(:) deallocate(RCHIXLIST, stat=ierr) - - ! sUps dimension + + ! sUps dimension call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachStart', iUpRchStart, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachCount', nUpRchCount, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) - + allocate(NETOPO(iSeg)%UREACHI(nUpRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%UREACHI') allocate(NETOPO(iSeg)%UREACHK(nUpRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%UREACHK') allocate(NETOPO(iSeg)%goodBas(nUpRchCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%goodBas') if (nUpRchCount > 0) then - + call get_vec_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachID' ,NETOPO(iSeg)%UREACHK(:),iUpRchStart,nUpRchCount,ierr,cmessage); call handle_err(ierr,cmessage) do jSeg=1,nUpRchCount ! Identify the index of the desired stream segment from reachID vector (dimension size: nSeg) iSelect = minloc(abs(NETOPO(:)%REACHID - NETOPO(iSeg)%UREACHK(jSeg))) NETOPO(iSeg)%UREACHI(jSeg) = iSelect(1) ! de-vectorize the desired stream segment ! check that we identify the correct upstream reach - if (NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID /= NETOPO(iSeg)%UREACHK(jSeg)) then - print*,'iSeg = ', iSeg + if (NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID /= NETOPO(iSeg)%UREACHK(jSeg)) then + print*,'iSeg = ', iSeg print*,'NETOPO(iSeg)%UREACHK(jSeg) = ', NETOPO(iSeg)%UREACHK(jSeg) print*,'NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID = ', NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID call handle_err(20,'unable to find desired immediate upstream segment') endif - + ! check that the upstream reach has a basin area > 0 if(RPARAM(NETOPO(iSeg)%UREACHI(jSeg))%TOTAREA > verySmall)then NETOPO(iSeg)%goodBas(jSeg) = .true. else NETOPO(iSeg)%goodBas(jSeg) = .false. endif - + enddo ! looping through the immediate upstream reaches endif ! if not a headwater !print*,'NETOPO(iSeg)%UREACHI(:) = ',NETOPO(iSeg)%UREACHI(:) !print*,'NETOPO(iSeg)%REACHID = ',NETOPO(iSeg)%REACHID !print*,'NETOPO(iSeg)%UREACHK(:) = ',NETOPO(iSeg)%UREACHK(:) - - ! sHrus dimension + + ! sHrus dimension call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upHruStart', iUpHruStart, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upHruCount', nUpHruCount, upStrmRchList(iSeg), ierr, cmessage); call handle_err(ierr,cmessage) allocate(h2b(iSeg)%cHRU(nUpHruCount), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for h2b(iSeg)%cHRU(nUpHruCount)') @@ -416,7 +416,7 @@ program route_runoff else ! if the entire river network routing is selected print*, 'Route all the segments included in network topology' - ! Populate sSeg dimensioned variable + ! Populate sSeg dimensioned variable allocate(NETOPO(nSeg), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO') allocate(RPARAM(nSeg), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for RPARAM') do iSeg=1,nSeg @@ -436,23 +436,23 @@ program route_runoff call get_scl_dvar(trim(ancil_dir)//trim(fname_ntop),'basinArea', RPARAM(iSeg)%BASAREA,iSeg,ierr,cmessage); call handle_err(ierr,cmessage) call get_scl_dvar(trim(ancil_dir)//trim(fname_ntop),'totalArea', RPARAM(iSeg)%TOTAREA,iSeg,ierr,cmessage); call handle_err(ierr,cmessage) enddo - ! Populate sAll dimensioned variable + ! Populate sAll dimensioned variable ! NETOPO%RCHLIST - upstream reach list - ! NETOPO%UPSLENG - total upstream reach length + ! NETOPO%UPSLENG - total upstream reach length nTotal=0 do iSeg=1,nSeg call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachStart', iRchStart1, iSeg, ierr, cmessage); call handle_err(ierr,cmessage) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'reachCount', nRchCount1, iSeg, ierr, cmessage); call handle_err(ierr,cmessage) allocate(NETOPO(iSeg)%UPSLENG(nRchCount1), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%UPSLENG') allocate(NETOPO(iSeg)%RCHLIST(nRchCount1), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for NETOPO%RCHLIST') - + call get_vec_dvar(trim(ancil_dir)//trim(fname_ntop),'upReachTotalLength',NETOPO(iSeg)%UPSLENG(:),iRchStart1,nRchCount1,ierr,cmessage); call handle_err(ierr,cmessage) call get_vec_ivar(trim(ancil_dir)//trim(fname_ntop),'reachList',NETOPO(iSeg)%RCHLIST(:),iRchStart1,nRchCount1,ierr,cmessage); call handle_err(ierr,cmessage) nTotal = nTotal + nRchCount1 - enddo - ! Populate sUps dimensioned variable + enddo + ! Populate sUps dimensioned variable ! NETOPO%UREACHI - Immediate upstream reach index list - ! NETOPO%UREACHK - Immediate upstream reach ID list + ! NETOPO%UREACHK - Immediate upstream reach ID list do iSeg=1,nSeg call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachStart', iUpRchStart, iSeg, ierr, cmessage); call handle_err(ierr,cmessage) call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachCount', nUpRchCount, iSeg, ierr, cmessage); call handle_err(ierr,cmessage) @@ -465,8 +465,8 @@ program route_runoff call get_vec_ivar(trim(ancil_dir)//trim(fname_ntop),'upReachIndex',NETOPO(iSeg)%UREACHI(:),iUpRchStart,nUpRchCount,ierr,cmessage); call handle_err(ierr,cmessage) do jSeg=1,nUpRchCount ! check that we identify the correct upstream reach - if (NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID /= NETOPO(iSeg)%UREACHK(jSeg)) then - print*,'iSeg = ', iSeg + if (NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID /= NETOPO(iSeg)%UREACHK(jSeg)) then + print*,'iSeg = ', iSeg print*,'NETOPO(iSeg)%UREACHK(jSeg) = ', NETOPO(iSeg)%UREACHK(jSeg) print*,'NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID = ', NETOPO(NETOPO(iSeg)%UREACHI(jSeg))%REACHID call handle_err(20,'unable to find desired immediate upstream segment') @@ -479,8 +479,8 @@ program route_runoff endif enddo ! looping through the immediate upstream reaches endif ! if not a headwater - enddo - ! Populate sHrus dimensioned variable + enddo + ! Populate sHrus dimensioned variable allocate(h2b(nSeg), stat=ierr); if(ierr/=0) call handle_err(ierr,'problem allocating space for h2b') do iSeg=1,nSeg call get_scl_ivar(trim(ancil_dir)//trim(fname_ntop),'upHruStart', iUpHruStart, iSeg, ierr, cmessage); call handle_err(ierr,cmessage) @@ -510,7 +510,7 @@ program route_runoff !do iRch=1,nSegRoute ! jRch = NETOPO(iRch)%RHORDER ! write(*,'(a,1x,2(i4,1x),f20.2)') 'iRch, NETOPO(jRch)%DREACHI, RPARAM(jRch)%TOTAREA/1000000._dp = ', iRch, NETOPO(jRch)%DREACHI, RPARAM(jRch)%TOTAREA/1000000._dp - !end do + !end do !pause ' check reachorder' end if @@ -531,10 +531,10 @@ program route_runoff !check ! !do iSeg=1949,1949 ! do iSeg=9,9 - ! nUpstream = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment + ! nUpstream = size(NETOPO(iSeg)%RCHLIST) ! size of upstream segment ! do iUps=1,nUpstream ! jUps=NETOPO(iSeg)%RCHLIST(iUps) - ! nTDH= size(NETOPO(iSeg)%UH(iUps)%UH_DATA) ! size of UH data + ! nTDH= size(NETOPO(iSeg)%UH(iUps)%UH_DATA) ! size of UH data ! write(*,'(a,1x,i4,1x,i4.1,1x,f10.2)') 'upstrm index, upstrmID,length = ', NETOPO(iSeg)%RCHLIST(iUps), NETOPO(jUps)%REACHID, NETOPO(iSeg)%UPSLENG(iUps) ! write(*,'(96f6.3)') (NETOPO(iSeg)%UH(iUps)%UH_DATA(iTDH), iTDH=1,nTDH) ! enddo @@ -572,8 +572,8 @@ program route_runoff call handle_err(ierr, cmessage) print*,'size(qsimHRUid) = ', size(qsimHRUid) -! check all the upstream hrus at the desired outlet exist in runoff file -! Assign hru_ix based on order of hru in runoff file +! check all the upstream hrus at the desired outlet exist in runoff file +! Assign hru_ix based on order of hru in runoff file do iSeg=1,nSegRoute ! (loop through stream segments) nDrain = size(h2b(iSeg)%cHRU) ! (number of HRUs that drain into a given stream segment) if(nDrain > 0)then @@ -582,16 +582,16 @@ program route_runoff if( minval(abs(qsimHRUid - h2b(iSeg)%cHRU(iHRU)%hru_id)) /= 0 )then write (str,'(I10)') h2b(iSeg)%cHRU(iHRU)%hru_id call handle_err(20,'runoff file does not include runoff time series at HRU'//trim(str)) - end if - ! Assign hru index + end if + ! Assign hru index iSelect = minloc(abs(qsimHRUid - h2b(iSeg)%cHRU(iHRU)%hru_id)) h2b(iSeg)%cHRU(iHRU)%hru_ix=iSelect(1) if(h2b(iSeg)%cHRU(iHRU)%hru_id /= qsimHRUid(h2b(iSeg)%cHRU(iHRU)%hru_ix))then call handle_err(20,'mis-match in HRUs') endif - qsimHRUid_mask(iSelect(1)) = .true. + qsimHRUid_mask(iSelect(1)) = .true. qsimHRUarea(iSelect(1)) = h2b(iSeg)%cHRU(iHRU)%hru_area - end do ! (loop through HRUs that drain into a given stream segment) + end do ! (loop through HRUs that drain into a given stream segment) endif ! (if HRUs drain into the stream segment) end do ! (loop through stream segments) @@ -643,7 +643,7 @@ program route_runoff ! initialize the routed elements RCHFLX(:,:)%BASIN_QR(1) = 0._dp -RCHFLX(:,:)%BASIN_QR_IRF(1) = 0._dp +RCHFLX(:,:)%BASIN_QR_IRF(1) = 0._dp ! initialize the time-delay histogram do iens=1,nens @@ -653,14 +653,14 @@ program route_runoff call handle_err(ierr, 'problem allocating space for QFUTURE element') ! initialize to zeroes RCHFLX(iens,ibas)%QFUTURE(:) = 0._dp - + if (doIRFroute) then ! allocate space for the delayed runoff - nUpstream = size(NETOPO(ibas)%RCHLIST) ! size of upstream segment + nUpstream = size(NETOPO(ibas)%RCHLIST) ! size of upstream segment nUH_DATA_MAX=0 upstrms_loop: do iUps=1,nUpstream nUH_DATA_MAX= max(nUH_DATA_MAX ,size(NETOPO(ibas)%UH(iUps)%UH_DATA)) - enddo upstrms_loop + enddo upstrms_loop allocate(RCHFLX(iens,ibas)%QFUTURE_IRF(nUH_DATA_MAX), stat=ierr) call handle_err(ierr, 'problem allocating space for QFUTURE_IRF element') ! initialize to zeroes @@ -726,7 +726,7 @@ program route_runoff write (str,'(I10)') h2b(ibas)%cHRU(iHRU)%hru_id call handle_err(20,'negaive runoff value is invalid') endif - + ! compute the weighted average qsim_basin(ibas) = qsim_basin(ibas) + h2b(ibas)%cHRU(iHRU)%wght*qsim_hru(ix)*time_conv*length_conv ! ensure m/s end do ! (looping through gridpoints associated with each basin) @@ -742,7 +742,7 @@ program route_runoff ! write instantaneous local runoff in each stream segment (m3/s) call write_dVec(trim(output_dir)//trim(fname_output), 'instBasinRunoff', RCHFLX(iens,:)%BASIN_QI, (/1,iTime/), (/nSegRoute,1/), ierr, cmessage) call handle_err(ierr,cmessage) - + ! ***** ! (5c) Delay runoff within local basins... ! **************************************** @@ -760,24 +760,24 @@ program route_runoff ! move array back do JTIM=2,NTDH RCHFLX(iens,ibas)%QFUTURE(JTIM-1) = RCHFLX(iens,ibas)%QFUTURE(JTIM) - end do + end do RCHFLX(iens,ibas)%QFUTURE(NTDH) = 0._dp end do ! (looping through basins) ! write routed local runoff in each stream segment (m3/s) call write_dVec(trim(output_dir)//trim(fname_output), 'dlayBasinRunoff', RCHFLX(iens,:)%BASIN_QR(1), (/1,iTime/), (/nSegRoute,1/), ierr, cmessage) call handle_err(ierr,cmessage) - + ! ***** ! (5d) Route streamflow for each upstream segment through the river network... ! ************************************************** if (doIRFroute) then - ! Get upstreme routed runoff depth at top of each segment for IRF routing + ! Get upstreme routed runoff depth at top of each segment for IRF routing call get_upsbas_qr(nSegRoute,iens,ierr,cmessage) ! write routed runoff (m3/s) call write_dVec(trim(output_dir)//trim(fname_output), 'UpBasRoutedRunoff', RCHFLX(iens,:)%UPSBASIN_QR, (/1,iTime/), (/nSegRoute,1/), ierr, cmessage) call handle_err(ierr,cmessage) - ! Get upstreme routed runoff depth at top of each segment for IRF routing + ! Get upstreme routed runoff depth at top of each segment for IRF routing call conv_upsbas_qr(nSegRoute,iens,ierr,cmessage) ! write routed runoff (m3/s) call write_dVec(trim(output_dir)//trim(fname_output), 'IRFroutedRunoff', RCHFLX(iens,:)%REACH_Q_IRF, (/1,iTime/), (/nSegRoute,1/), ierr, cmessage) @@ -795,10 +795,10 @@ program route_runoff allocate(iUpstream(nUpstream), qUpstream(nUpstream), stat=ierr) if(ierr/=0) call handle_err(ierr,'problem allocating vectors for all upstream basins') ! get indices for all reaches upstream - iUpstream = NETOPO(iSeg)%RCHLIST(1:nUpstream) + iUpstream(1:nUpstream) = NETOPO(iSeg)%RCHLIST(1:nUpstream) ! get streamflow for all reaches upstream - qUpstream = RCHFLX(iens,iUpstream(1:nUpstream))%BASIN_QR(1) - ! get mean streamflow + qUpstream(1:nUpstream) = RCHFLX(iens,iUpstream(1:nUpstream))%BASIN_QR(1) + ! get mean streamflow RCHFLX(IENS,iSeg)%UPSTREAM_QI = sum(qUpstream) ! test if(NETOPO(iSeg)%REACHID == ixPrint)then @@ -839,7 +839,7 @@ program route_runoff !if(iRch==5) pause 'finished stream segment' end do ! (looping through stream segments) !pause 'first time step' - + ! write routed runoff (m3/s) call write_dVec(trim(output_dir)//trim(fname_output), 'KWTroutedRunoff', RCHFLX(iens,:)%REACH_Q, (/1,iTime/), (/nSegRoute,1/), ierr, cmessage) call handle_err(ierr,cmessage) @@ -885,19 +885,19 @@ end subroutine handle_err subroutine findix(scl,vec,iSelect) - ! Find vec index where the value match up with scl + ! Find vec index where the value match up with scl implicit none - - integer(i4b),intent(in) :: scl + + integer(i4b),intent(in) :: scl integer(i4b),intent(in),allocatable :: vec(:) integer(i4b),intent(out) :: iSelect - + integer(i4b) :: i(1) i = minloc(abs(vec - scl)) iSelect = i(1) ! de-vectorize the desired stream segment if(vec(iSelect) /= scl)& call handle_err(20,'unable to find matched value') - end subroutine findix + end subroutine findix end diff --git a/route/build/src/var_lookup.f90 b/route/build/src/var_lookup.f90 index 1267cd61..5b54c21f 100644 --- a/route/build/src/var_lookup.f90 +++ b/route/build/src/var_lookup.f90 @@ -30,7 +30,7 @@ MODULE var_lookup ! (3) define variables for the hru2basin mapping ! *********************************************************************************************************** type, public :: iLook_MAP - integer(i4b) :: HRUid = 1 ! HRU id + integer(i4b) :: HRUid = 1 ! HRU id integer(i4b) :: segHRUid = 2 ! id of the stream segment below each HRU endtype iLook_MAP ! *********************************************************************************************************** @@ -43,10 +43,10 @@ MODULE var_lookup ! *********************************************************************************************************** ! (X) define size of data vectors ! *********************************************************************************************************** - integer(i4b),parameter,public :: nVarsHRU=4 - integer(i4b),parameter,public :: nVarsSEG=7 - integer(i4b),parameter,public :: nVarsMAP=2 - integer(i4b),parameter,public :: nVarsTOP=2 + integer(i4b),parameter,public :: nVarsHRU=4 + integer(i4b),parameter,public :: nVarsSEG=7 + integer(i4b),parameter,public :: nVarsMAP=2 + integer(i4b),parameter,public :: nVarsTOP=2 ! *********************************************************************************************************** ! (X) define data vectors ! ***********************************************************************************************************