From 8466a7e4e12f8befca7729cdfa05938a33e29b44 Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Tue, 19 Nov 2019 13:24:29 -0700 Subject: [PATCH 1/3] Add Yates' updates to Reach-Lakes --- trunk/NDHMS/Data_Rec/module_namelist.F | 14 +++++++------- trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F | 1 + trunk/NDHMS/Routing/module_HYDRO_io.F | 2 +- trunk/NDHMS/Routing/module_RT.F | 6 +++--- trunk/NDHMS/Routing/module_channel_routing.F | 6 +++--- 5 files changed, 15 insertions(+), 14 deletions(-) diff --git a/trunk/NDHMS/Data_Rec/module_namelist.F b/trunk/NDHMS/Data_Rec/module_namelist.F index 310ed99c2..37c0b4c5d 100644 --- a/trunk/NDHMS/Data_Rec/module_namelist.F +++ b/trunk/NDHMS/Data_Rec/module_namelist.F @@ -713,13 +713,13 @@ subroutine rt_nlst_check(nlst) if (.not. fileExists) call hydro_stop('hydro.namelist ERROR: route_lake_f not found.') endif ! Only allow lakes to be ran with gridded routing or NWM routing - if(len(trim(nlst%route_lake_f)) .ne. 0) then - if(nlst%channel_option .ne. 3) then - if(nlst%UDMP_OPT .ne. 1) then - call hydro_stop('hydro.namelist ERROR: Currently lakes only work with gridded channel routing or UDMP=1. Please change your namelist settings.') - endif - endif - endif +! if(len(trim(nlst%route_lake_f)) .ne. 0) then +! if(nlst%channel_option .ne. 3) then +! if(nlst%UDMP_OPT .ne. 1) then +! call hydro_stop('hydro.namelist ERROR: Currently lakes only work with gridded channel routing or UDMP=1. Please change your namelist settings.') +! endif +! endif +! endif if((nlst%channel_option .eq. 3) .and. (nlst%compound_channel)) then call hydro_stop("Compound channel option not available for diffusive wave routing. ") diff --git a/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F b/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F index ebb4eea73..c5f4c76b0 100644 --- a/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F +++ b/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F @@ -449,6 +449,7 @@ subroutine HYDRO_out(did, rstflag) call output_chrtout_grd_NWM(did,nlst_rt(did)%igrid) endif endif + if (RT_DOMAIN(did)%NLAKES.gt.0) then if(nlst_rt(did)%io_form_outputs .ne. 0) then ! Output lakes in NWM format diff --git a/trunk/NDHMS/Routing/module_HYDRO_io.F b/trunk/NDHMS/Routing/module_HYDRO_io.F index 7eb06853d..5fe28bb90 100644 --- a/trunk/NDHMS/Routing/module_HYDRO_io.F +++ b/trunk/NDHMS/Routing/module_HYDRO_io.F @@ -1004,7 +1004,7 @@ SUBROUTINE READ_ROUTEDIM(IXRT,JXRT,route_chan_f,route_link_f, & !!-- no longer find the lakes from the 2-d hi res grid !DJG inv do j=jxrt,1,-1 ! follwoing is modified by Wei Yu 03/24/2017 - if(UDMP_OPT .eq. 0) then + if(UDMP_OPT .eq. 0 .and. channel_option .eq. 3) then NLAKES = 0 do j=1,jxrt do i = 1,ixrt diff --git a/trunk/NDHMS/Routing/module_RT.F b/trunk/NDHMS/Routing/module_RT.F index 5b0ad441f..86f9e2018 100644 --- a/trunk/NDHMS/Routing/module_RT.F +++ b/trunk/NDHMS/Routing/module_RT.F @@ -532,7 +532,7 @@ subroutine getChanDim(did) rt_domain(did)%GNLINKSL = 1 rt_domain(did)%NLINKSL = 1 endif - if(nlst_rt(did)%UDMP_OPT .eq. 1) & + if(nlst_rt(did)%UDMP_OPT .eq. 1 .or. nlst_rt(did)%channel_option .eq. 1 .or. nlst_rt(did)%channel_option .eq. 2) & call read_NSIMLAKES(rt_domain(did)%NLAKES,nlst_rt(did)%route_lake_f) call rt_allocate(did,rt_domain(did)%ix,rt_domain(did)%jx,& @@ -588,7 +588,7 @@ subroutine getChanDim(did) endif -if(nlst_rt(did)%UDMP_OPT .eq. 1) then +if(nlst_rt(did)%UDMP_OPT .eq. 1 .or. nlst_rt(did)%channel_option .eq. 1 .or. nlst_rt(did)%channel_option .eq. 2) then call read_NSIMLAKES(rt_domain(did)%NLAKES,nlst_rt(did)%route_lake_f) endif @@ -1460,7 +1460,7 @@ subroutine LandRT_ini(did) endif !End if channel option eq 3; else; -! Initialize Lake Elevations for Gridded and NWM routing. +! Initialize Lake Elevations for Gridded and NWM routing. AND REACH Community? do j=1,rt_domain(did)%NLAKES rt_domain(did)%RESHT(j) = rt_domain(did)%ORIFICEE(j) + & ((rt_domain(did)%LAKEMAXH(j) - rt_domain(did)%ORIFICEE(j) )* rt_domain(did)%ELEVLAKE(j)) diff --git a/trunk/NDHMS/Routing/module_channel_routing.F b/trunk/NDHMS/Routing/module_channel_routing.F index 404aab25f..e868d4934 100644 --- a/trunk/NDHMS/Routing/module_channel_routing.F +++ b/trunk/NDHMS/Routing/module_channel_routing.F @@ -964,9 +964,9 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & if(TYPEL(k) .eq. 1) then !--link is a reservoir - ! CALL LEVELPOOL(1,QLINK(k,1), Qup, QLINK(k,1), QLINK(k,2), & - ! QLateral(k), DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), & - ! WEIRC(k), WEIRL(k),ORIFICEE(k), ORIFICEC(k), ORIFICEA(k)) + CALL LEVELPOOL(linkid(k),Quc, Qup, QLINK(k,1), & + QLateral(k), DT, RESHT(k), HRZAREA(k), WEIRH(k), LAKEMAXH(k), & + WEIRC(k), WEIRL(k), ORIFICEE(k), ORIFICEC(k), ORIFICEA(k)) elseif (channel_option .eq. 1) then !muskingum routing Km = MUSK(k) From 966df72350b9531b6cd2891202b325a7394e5620 Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Wed, 20 Nov 2019 16:25:05 -0700 Subject: [PATCH 2/3] Update NCAR Reach routing for lakes * Add call to LEVELPOOL for NCAR Reach routing in the DRIVE_CHANNEL subroutine in module_channel_routing.F * Add update to DRIVE_CHANNEL interface to include RT_DOMAIN(did)%LAKEIDX for mapping links indices to lake indices --- trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F | 14 ++-- trunk/NDHMS/Routing/module_channel_routing.F | 74 +++++++++++++++----- 2 files changed, 63 insertions(+), 25 deletions(-) diff --git a/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F b/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F index c5f4c76b0..035387e45 100644 --- a/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F +++ b/trunk/NDHMS/HYDRO_drv/module_HYDRO_drv.F @@ -985,16 +985,16 @@ subroutine driveChannelRouting(did) RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, & RT_DOMAIN(did)%CHANYJ, nlst_rt(did)%channel_option, & RT_DOMAIN(did)%RETDEP_CHAN, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, & - RT_DOMAIN(did)%node_area & + RT_DOMAIN(did)%node_area, RT_DOMAIN(did)%LAKEIDX, & #ifdef MPP_LAND - ,RT_DOMAIN(did)%lake_index,RT_DOMAIN(did)%link_location,& + RT_DOMAIN(did)%lake_index,RT_DOMAIN(did)%link_location,& RT_DOMAIN(did)%mpp_nlinks,RT_DOMAIN(did)%nlinks_index, & - RT_DOMAIN(did)%yw_mpp_nlinks & - , RT_DOMAIN(did)%LNLINKSL,RT_DOMAIN(did)%LLINKID & - , rt_domain(did)%gtoNode,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd & + RT_DOMAIN(did)%yw_mpp_nlinks, & + RT_DOMAIN(did)%LNLINKSL,RT_DOMAIN(did)%LLINKID, & + rt_domain(did)%gtoNode,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd, & #endif - , rt_domain(did)%CH_LNKRT_SL & - ,nlst_rt(did)%GwBaseSwCRT, gw2d(did)%ho, gw2d(did)%qgw_chanrt, & + rt_domain(did)%CH_LNKRT_SL, & + nlst_rt(did)%GwBaseSwCRT, gw2d(did)%ho, gw2d(did)%qgw_chanrt, & nlst_rt(did)%gwChanCondSw, nlst_rt(did)%gwChanCondConstIn, & nlst_rt(did)%gwChanCondConstOut,rt_domain(did)%velocity & ) diff --git a/trunk/NDHMS/Routing/module_channel_routing.F b/trunk/NDHMS/Routing/module_channel_routing.F index e868d4934..ee7c69019 100644 --- a/trunk/NDHMS/Routing/module_channel_routing.F +++ b/trunk/NDHMS/Routing/module_channel_routing.F @@ -699,7 +699,7 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & RESHT, HRZAREA, LAKEMAXH, WEIRH, WEIRC, WEIRL, ORIFICEC, ORIFICEA, & ORIFICEE, ZELEV, CVOL, NLAKES, QLAKEI, QLAKEO, LAKENODE, & dist, QINFLOWBASE, CHANXI, CHANYJ, channel_option, RETDEP_CHAN, & - NLINKSL, LINKID, node_area & + NLINKSL, LINKID, node_area, lake_lookup & #ifdef MPP_LAND , lake_index,link_location,mpp_nlinks,nlinks_index,yw_mpp_nlinks & , LNLINKSL, LLINKID & @@ -792,8 +792,10 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & REAL, DIMENSION(NLAKES) :: QLLAKE !-- lateral inflow to lake in diffusion scheme REAL*8, DIMENSION(NLAKES) :: QLLAKE8 !-- lateral inflow to lake in diffusion scheme + integer, intent(in), dimension(:) :: lake_lookup !-- inverse lake index for k->lake mapping + !-- Local Variables - INTEGER :: i,j,k,t,m,jj,kk,KRT,node + INTEGER :: i,j,k,t,m,jj,kk,KRT,node,l_idx INTEGER :: DT_STEPS !-- number of timestep in routing REAL :: Qup,Quc !--Q upstream Previous, Q Upstream Current, downstream Previous REAL :: bo !--critical depth, bnd outflow just for testing @@ -812,19 +814,28 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & real ywtmp(ixrt,jxrt) integer LNLINKSL integer, dimension(LNLINKSL) :: LLINKID - real*8, dimension(LNLINKSL) :: LQLateral -! real*4, dimension(LNLINKSL) :: LQLateral + real(kind=8), dimension(LNLINKSL) :: LQLateral integer, dimension(:) :: toNodeInd integer, dimension(:,:) :: gtoNode integer :: nToNodeInd real, dimension(nToNodeInd,2) :: gQLINK + real, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI, tmpRESHT #else - real*8, dimension(NLINKS) :: LQLateral !--lateral flow + real(kind=8), dimension(NLINKS) :: LQLateral !--lateral flow #endif integer flag - integer :: n, kk2, nt, nsteps ! tmp + integer :: n, kk2, nt, nsteps ! tmp +#ifdef MPP_LAND + if(my_id == io_id) then +#endif + allocate(tmpQLAKEO(NLAKES)) + allocate(tmpQLAKEI(NLAKES)) + allocate(tmpRESHT(NLAKES)) +#ifdef MPP_LAND + endif +#endif QLAKEIP = 0 QLAKEI8 = 0 HLINKTMP = 0 @@ -932,6 +943,15 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & !---------- route other reaches, with upstream inflow tmpQlink = 0.0 +#ifdef MPP_LAND + if(my_id .eq. io_id) then +#endif + tmpQLAKEO = QLAKEO + tmpQLAKEI = QLAKEI + tmpRESHT = RESHT +#ifdef MPP_LAND + endif +#endif do k = 1,NLINKSL ! if (ORDER(k) .gt. 1 ) then !-- exclude first order stream Quc = 0.0 @@ -962,29 +982,41 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & end do ! do m #endif - if(TYPEL(k) .eq. 1) then !--link is a reservoir - - CALL LEVELPOOL(linkid(k),Quc, Qup, QLINK(k,1), & - QLateral(k), DT, RESHT(k), HRZAREA(k), WEIRH(k), LAKEMAXH(k), & - WEIRC(k), WEIRL(k), ORIFICEE(k), ORIFICEC(k), ORIFICEA(k)) - - elseif (channel_option .eq. 1) then !muskingum routing + if(TYPEL(k) == 1) then !--link is a reservoir + l_idx = lake_lookup(k) + if (l_idx >= 0) then !-- -999 if not a reservoir in the lookup table (belt-and-suspenders check) + print *, "CALL LEVELPOOL", l_idx,Quc, Qup, QLINK(k,2), QLateral(k), DT, & + RESHT(l_idx), HRZAREA(l_idx), WEIRH(l_idx), LAKEMAXH(l_idx), & + WEIRC(l_idx), WEIRL(l_idx), ORIFICEE(l_idx), ORIFICEC(l_idx), ORIFICEA(l_idx) + CALL LEVELPOOL(l_idx,Quc, Qup, QLINK(k,2), QLateral(k), DT, & + RESHT(l_idx), HRZAREA(l_idx), WEIRH(l_idx), LAKEMAXH(l_idx), & + WEIRC(l_idx), WEIRL(l_idx), ORIFICEE(l_idx), ORIFICEC(l_idx), ORIFICEA(l_idx)) + QLAKEO(l_idx) = QLINK(k,2) !save outflow to lake + QLAKEI(l_idx) = Quc !save inflow to lake + end if + elseif (channel_option .eq. 1) then !muskingum routing Km = MUSK(k) X = MUSX(k) tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow - elseif (channel_option .eq. 2) then ! muskingum cunge + elseif (channel_option .eq. 2) then ! muskingum cunge call SUBMUSKINGCUNGE(tmpQLINK(k,2), velocity(k), LINKID(k), & Qup,Quc, QLINK(k,1), QLateral(k), DTRT_CH, So(k), & CHANLEN(k), MannN(k), ChSSlp(k), Bw(k), Tw(k),Tw_CC(k), n_CC(k), HLINK(k) ) - else + else print *, "FATAL ERROR: no channel option selected" call hydro_stop("In drive_CHANNEL() - no channel option selected") - endif + endif ! endif !!! order(1) .ne. 1 end do !--k links + #ifdef MPP_LAND + call updateLake_seq(RESHT,nlakes,tmpRESHT) + call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO) + call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI) + #endif + !yw check ! gQLINK = 0.0 ! call ReachLS_write_io(tmpQLINK(:,2), gQLINK(:,2)) @@ -997,7 +1029,7 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & ! endif do k = 1, NLINKSL - if(TYPEL(k) .ne. 1) then + if(TYPEL(k) .ne. 2) then QLINK(k,2) = tmpQLINK(k,2) endif QLINK(k,1) = QLINK(k,2) !assing link flow of current to be previous for next time step @@ -1416,8 +1448,14 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & #endif if (KT .eq. 1) KT = KT + 1 - +#ifdef MPP_LAND + if (my_id == io_id) then + if(allocated(tmpRESHT)) deallocate(tmpRESHT) + if(allocated(tmpQLAKEO)) deallocate(tmpQLAKEO) + if(allocated(tmpQLAKEI)) deallocate(tmpQLAKEI) + endif +#endif end subroutine drive_CHANNEL ! ---------------------------------------------------------------- From 0038ccee394142c045519591c97a0278e73810cc Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Thu, 21 Nov 2019 09:19:03 -0700 Subject: [PATCH 3/3] Remove debug statement and left-align #ifdef * Remove extraneous debug print statement * Left-align an #ifdef to make gfortran happy --- trunk/NDHMS/Routing/module_channel_routing.F | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/trunk/NDHMS/Routing/module_channel_routing.F b/trunk/NDHMS/Routing/module_channel_routing.F index ee7c69019..e1fdcc18f 100644 --- a/trunk/NDHMS/Routing/module_channel_routing.F +++ b/trunk/NDHMS/Routing/module_channel_routing.F @@ -985,9 +985,6 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & if(TYPEL(k) == 1) then !--link is a reservoir l_idx = lake_lookup(k) if (l_idx >= 0) then !-- -999 if not a reservoir in the lookup table (belt-and-suspenders check) - print *, "CALL LEVELPOOL", l_idx,Quc, Qup, QLINK(k,2), QLateral(k), DT, & - RESHT(l_idx), HRZAREA(l_idx), WEIRH(l_idx), LAKEMAXH(l_idx), & - WEIRC(l_idx), WEIRL(l_idx), ORIFICEE(l_idx), ORIFICEC(l_idx), ORIFICEA(l_idx) CALL LEVELPOOL(l_idx,Quc, Qup, QLINK(k,2), QLateral(k), DT, & RESHT(l_idx), HRZAREA(l_idx), WEIRH(l_idx), LAKEMAXH(l_idx), & WEIRC(l_idx), WEIRL(l_idx), ORIFICEE(l_idx), ORIFICEC(l_idx), ORIFICEA(l_idx)) @@ -1011,11 +1008,11 @@ Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & ! endif !!! order(1) .ne. 1 end do !--k links - #ifdef MPP_LAND +#ifdef MPP_LAND call updateLake_seq(RESHT,nlakes,tmpRESHT) call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO) call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI) - #endif +#endif !yw check ! gQLINK = 0.0