Merge pull request #428 from rcabell/update_reach_lakes

Add support for Lakes in NCAR Reach (v5.1.1)
NCAR · Nov 21, 2019 · 15cd3d9 · 15cd3d9
2 parents 3b7a386 + 0038cce
commit 15cd3d9
Show file tree

Hide file tree

Showing 5 changed files with 72 additions and 36 deletions.
diff --git 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
@@ -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
@@ -984,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_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
@@ -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
@@ -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,38 @@ 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(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))
-
-                   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)
+                        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 +1026,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 +1445,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
 ! ----------------------------------------------------------------