!  Program Name:
!  Author(s)/Contact(s):
!  Abstract:
!  History Log:
! 
!  Usage:
!  Parameters: <Specify typical arguments passed>
!  Input Files:
!        <list file names and briefly describe the data they include>
!  Output Files:
!        <list file names and briefly describe the information they include>
! 
!  Condition codes:
!        <list exit condition or error codes returned >
!        If appropriate, descriptive troubleshooting instructions or
!        likely causes for failures could be mentioned here with the
!        appropriate error code
! 
!  User controllable options: <if applicable>

MODULE module_channel_routing
#ifdef MPP_LAND
  USE module_mpp_land
  use MODULE_mpp_ReachLS, only : updatelinkv,                   & 
                                 ReachLS_write_io, gbcastvalue, &
                                 gbcastreal2

#endif
  IMPLICIT NONE

  contains

! ------------------------------------------------
!   FUNCTION MUSKING
! ------------------------------------------------
	REAL FUNCTION MUSKING(idx,qup,quc,qdp,dt,Km,X)

	IMPLICIT NONE

!--local variables
        REAL    :: C1, C2, C3
        REAL    :: Km         !K travel time in hrs in reach
        REAL    :: X          !weighting factors 0<=X<=0.5 
        REAL    :: dt         !routing period in hrs
        REAL    :: avgbf      !average base flow for initial condition
        REAL    :: qup        !inflow from previous timestep
        REAL    :: quc        !inflow  of current timestep
        REAL    :: qdp        !outflow of previous timestep
        REAL    :: dth        !timestep in hours
        INTEGER :: idx       ! index

        dth = dt/3600    !hours in timestep
        C1 = (dth - 2*Km*X)/(2*Km*(1-X)+dth)
        C2 = (dth+2*Km*X)/(2*Km*(1-X)+dth)
        C3 = (2*Km*(1-X)-dth)/(2*Km*(1-X)+dth)
        MUSKING = (C1*quc)+(C2*qup)+(C3*qdp)

! ----------------------------------------------------------------
  END FUNCTION MUSKING
! ----------------------------------------------------------------

! ------------------------------------------------
!   SUBROUTINE LEVELPOOL
! ------------------------------------------------

SUBROUTINE LEVELPOOL(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,oe,oc,oa)

    !! ----------------------------  argument variables
    !! All elevations should be relative to a common base (often belev(k))

    real, intent(INOUT) :: H       ! water elevation height (m)
    real, intent(IN)    :: dt      ! routing period [s]
    real, intent(IN)    :: qi0     ! inflow at previous timestep (cms)
    real, intent(IN)    :: qi1     ! inflow at current timestep (cms)
    real, intent(OUT)   :: qo1     ! outflow at current timestep
    real, intent(IN)    :: ql      ! lateral inflow
    real, intent(IN)    :: ar      ! area of reservoir (km^2)
    real, intent(IN)    :: we      ! bottom of weir elevation
    real, intent(IN)    :: wc      ! weir coeff.
    real, intent(IN)    :: wl      ! weir length (m)
    real, intent(IN)    :: oe      ! orifice elevation
    real, intent(IN)    :: oc      ! orifice coeff.
    real, intent(IN)    :: oa      ! orifice area (m^2)
    real, intent(IN)    :: maxh    ! max depth of reservoir before overtop (m)                     
    integer, intent(IN) :: ln      ! lake number

    !!DJG Add lake option switch here...move up to namelist in future versions...
    integer :: LAKE_OPT            ! Lake model option (move to namelist later)
    real    :: Htmp                ! Temporary assign of incoming lake el. (m)

    !! ----------------------------  local variables
    real :: sap                    ! local surface area values
    real :: discharge              ! storage discharge m^3/s
    real :: tmp1, tmp2
    real :: dh, dh1, dh2, dh3      ! height function and 3 order RK
    real :: It, Itdt_3, Itdt_2_3
    real :: maxWeirDepth           !maximum capacity of weir
    !! ----------------------------  subroutine body: from chow, mad mays. pg. 252
    !! -- determine from inflow hydrograph


    !!DJG Set hardwire for LAKE_OPT...move specification of this to namelist in
    !future versions...
    LAKE_OPT = 2
    Htmp = H   !temporary set of incoming lake water elevation...
    
    
    !!DJG IF-block for lake model option  1 - outflow=inflow, 2 - Chow et al level
    !pool, .....
    if (LAKE_OPT.eq.1) then     ! If-block for simple pass through scheme....
       
       qo1 = qi1                 ! Set outflow equal to inflow at current time      
       H = Htmp                  ! Set new lake water elevation to incoming lake el.
       
    else if (LAKE_OPT.eq.2) then   ! If-block for Chow et al level pool scheme
 
    It = qi0
    Itdt_3   = (qi0 + (qi1 + ql))/3
    Itdt_2_3 = (qi0 + (qi1 + ql))/3 + Itdt_3
    maxWeirDepth =  maxh - we   

    !-- determine Q(dh) from elevation-discharge relationship
    !-- and dh1
    dh = H - we
    if (dh .gt. maxWeirDepth) then 
       dh = maxWeirDepth 
    endif

    if (dh .gt. 0.0 ) then              !! orifice and overtop discharge
        tmp1 = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        tmp2 = wc * wl * (dh ** 2./3.)
        discharge = tmp1 + tmp2

        if (H .gt. 0.0) then
          sap = (ar * 1.0E6 ) * (1 + (H - we) / H)
        else
          sap  = 0.0
        endif

    else if ( H .gt. oe ) then     !! only orifice flow,not full
        discharge = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        sap = ar * 1.0E6
    else
        discharge = 0.0
        sap = ar * 1.0E6
    endif

    if (sap .gt. 0) then 
      dh1 = ((It - discharge)/sap)*dt
    else
      dh1 = 0.0
    endif

    !-- determine Q(H + dh1/3) from elevation-discharge relationship
    !-- dh2
    dh = (H+dh1/3) - we
    if (dh .gt. maxWeirDepth) then 
       dh = maxWeirDepth 
    endif

    if (dh .gt. 0.0 ) then              !! orifice and overtop discharge
        tmp1 = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        tmp2 = wc * wl * (dh ** 2./3.) 
        discharge = tmp1 + tmp2

        if (H .gt. 0.0) then 
         sap = (ar * 1.0E6 ) * (1 + (H - we) / H)
        else
         sap  = 0.0
        endif

    else if ( H .gt. oe ) then     !! only orifice flow,not full
        discharge = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        sap = ar * 1.0E6
    else
        discharge = 0.0
        sap = ar * 1.0E6
    endif

    if (sap .gt. 0.0) then 
     dh2 = ((Itdt_3 - discharge)/sap)*dt
    else
     dh2 = 0.0
    endif

    !-- determine Q(H + 2/3 dh2) from elevation-discharge relationship
    !-- dh3
    dh = (H + (0.667*dh2)) - we
    if (dh .gt. maxWeirDepth) then 
       dh = maxWeirDepth 
    endif

    if (dh .gt. 0.0 ) then              !! orifice and overtop discharge
        tmp1 = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        tmp2 = wc * wl * (dh ** 2./3.) 
        discharge = tmp1 + tmp2

        if (H .gt. 0.0) then
         sap = (ar * 1.0E6 ) * (1 + (H - we) / H)
        else
         sap = 0.0
        endif 

    else if ( H .gt. oe ) then     !! only orifice flow,not full
        discharge = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        sap = ar * 1.0E6
    else
        discharge = 0.0
        sap = ar * 1.0E6
    endif
    
    if (sap .gt. 0.0) then 
      dh3 = ((Itdt_2_3 - discharge)/sap)*dt
    else
      dh3 = 0.0
    endif

    !-- determine dh and H
    dh = (dh1/4.) + (0.75*dh3)
    H = H + dh

    !-- compute final discharge
    dh = H - we
    if (dh .gt. maxWeirDepth) then 
       dh = maxWeirDepth 
    endif
    if (dh .gt. 0.0 ) then              !! orifice and overtop discharge
        tmp1 = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        tmp2 = wc * wl * (dh ** 2./3.)
        discharge = tmp1 + tmp2

        if (H .gt. 0.0) then 
         sap = (ar * 1.0E6 ) * (1 + (H - we) / H)
        else
         sap = 0.0
        endif

    else if ( H .gt. oe ) then     !! only orifice flow,not full
        discharge = oc * oa * sqrt(2 * 9.81 * ( H - oe ) )
        sap = ar * 1.0E6
    else
        discharge = 0.0
        sap = ar * 1.0E6
    endif

    if(H .ge. maxh) then  ! overtop condition
     discharge = qi1
     H = maxh
    endif

    qo1  = discharge  ! return the flow rate from reservoir

23 format('botof H dh orf wr Q',f8.4,2x,f8.4,2x,f8.3,2x,f8.3,2x,f8.2)
24 format('ofonl H dh sap Q ',f8.4,2x,f8.4,2x,f8.0,2x,f8.2)

 
     ELSE   ! ELSE for LAKE_OPT....
     ENDIF  ! ENDIF for LAKE_OPT....
 
  RETURN

! ----------------------------------------------------------------
  END SUBROUTINE LEVELPOOL
! ----------------------------------------------------------------


! ------------------------------------------------
!   FUNCTION Diffusive wave
! ------------------------------------------------
        REAL FUNCTION DIFFUSION(nod,z1,z20,h1,h2,dx,n, &
                                Bw, Cs)
        IMPLICIT NONE
!-- channel geometry and characteristics
        REAL    :: Bw         !-bottom width (meters)
        REAL    :: Cs         !-Channel side slope slope
        REAL    :: dx         !-channel lngth (m)
        REAL,intent(in)    :: n          !-mannings coefficient
        REAL    :: R          !-Hydraulic radius
        REAL    :: AREA       !- wetted area
        REAL    :: h1,h2      !-tmp height variables
        REAL    :: z1,z2      !-z1 is 'from', z2 is 'to' elevations
        REAL    :: z          !-channel side distance
        REAL    :: w          !-upstream weight
        REAL    :: Ku,Kd      !-upstream and downstream conveyance
        REAL    :: Kf         !-final face conveyance
        REAL    :: Sf         !-friction slope
        REAL    :: sgn        !-0 or 1 
        INTEGER :: nod         !- node
        REAL ::  z20, dzx

! added by Wei Yu for bad data.

        dzx = (z1 - z20)/dx
        if(dzx .lt. 0.002) then
           z2 = z1 - dx*0.002  
        else
           z2 = z20
        endif
!end 

        if (n.le.0.0.or.Cs.le.0.or.Bw.le.0) then
         print *, "Error in Diffusion function ->channel coefficients"
         print *, "nod, n, Cs, Bw", nod, n, Cs, Bw 
         call hydro_stop("In DIFFUSION() - Error channel coefficients.")
        endif

!        Sf = ((z1+h1)-(z2+h2))/dx  !-- compute the friction slope
       !if(z1 .eq. z2) then
       ! Sf = ((z1-(z2-0.01))+(h1-h2))/dx  !-- compute the friction slope
       !else
!         Sf = ((z1-z2)+(h1-h2))/dx  !-- compute the friction slope
       !endif

!modifieed by Wei Yu for false geography data
         if(abs(z1-z2) .gt. 1.0E5) then
#ifdef HYDRO_D
             print*, "WARNING: huge slope rest to 0 for channel grid.", z1,z2
#endif
             Sf = ((h1-h2))/dx  !-- compute the friction slope
         else
             Sf = ((z1-z2)+(h1-h2))/dx  !-- compute the friction slope
         endif
!end  modfication

        sgn = SGNf(Sf)             !-- establish sign

        w = 0.5*(sgn + 1.)         !-- compute upstream or downstream weighting
        
        z = 1/Cs                   !--channel side distance (m)
        R = ((Bw+z*h1)*h1)/(Bw+2*h1*sqrt(1+z*z)) !-- Hyd Radius
        AREA = (Bw+z*h1)*h1        !-- Flow area
        Ku = (1/n)*(R**(2./3.))*AREA     !-- convenyance

        R = ((Bw+z*h2)*h2)/(Bw+2*h2*sqrt(1+z*z)) !-- Hyd Radius
        AREA = (Bw+z*h2)*h2        !-- Flow area
        Kd = (1/n)*(R**(2./3.))*AREA     !-- convenyance

        Kf =  (1-w)*Kd + w*Ku      !-- conveyance 
        DIFFUSION = Kf * sqrt(abs(Sf))*sgn


100     format('z1,z2,h1,h2,kf,Dif, Sf, sgn  ',f8.3,2x,f8.3,2x,f8.4,2x,f8.4,2x,f8.3,2x,f8.3,2x,f8.3,2x,f8.0)

  END FUNCTION DIFFUSION
! ----------------------------------------------------------------

! ------------------------------------------------
!   FUNCTION MUSKINGUM CUNGE
! ------------------------------------------------
        REAL FUNCTION MUSKINGCUNGE(idx,qup, quc, qdp, ql,&
                                   dt,So,dx,n,Cs,Bw)
        IMPLICIT NONE

!--local variables
        REAL    :: C1, C2, C3, C4
        REAL    :: Km          !K travel time in hrs in reach
        REAL    :: X          !weighting factors 0<=X<=0.5
        REAL    :: dt         !routing period in  seconds
        REAL    :: qup        !flow upstream previous timestep
        REAL    :: quc        !flow upstream current timestep
        REAL    :: qdp        !flow downstream previous timestep
!       REAL    :: qdc        !flow downstream current timestep
        REAL    :: ql         !lateral inflow through reach (m^3/sec)
        REAL    :: Ck         ! wave celerity (m/s)

!-- channel geometry and characteristics
        REAL    :: Bw         ! bottom width (meters)
        REAL    :: Cs         ! Channel side slope slope
        REAL    :: So         ! Channel bottom slope %
        REAL    :: dx         ! channel lngth (m)
        REAL    :: n          ! mannings coefficient
        REAL    :: Tw         ! top width at peak flow
        REAL    :: AREA       ! Cross sectional area m^2
        REAL    :: Z          ! trapezoid distance (m)
        REAL    :: R          ! Hydraulic radius
        REAL    :: WP         ! wetted perimmeter
        REAL    :: h          ! depth of flow
        REAL    :: h_0,h_1    ! secant method estimates
        REAL    :: Qj_0       ! secant method estimates
        REAL    :: Qj         ! intermediate flow estimate
        REAL    :: D,D1       ! diffusion coeff
        REAL    :: dtr        ! required timestep, minutes
        REAL    :: error
        REAL    :: hp         !courant, previous height
        INTEGER :: maxiter    !maximum number of iterations


!-- local variables.. needed if channel is sub-divded
        REAL    :: a,b,c
        INTEGER :: i,idx     !-- channel segment counter

!yw add
        goto 101
        C1   = 0
        C2   = 0
        C3   = 0
        C4   = 0
        Km   = 0
        X    = 0        !weighting factors 0<=X<=0.5
        Ck   = 0
        Tw   = 0       ! top width at peak flow
        AREA = 0      ! Cross sectional area m^2
        Z    = 0         ! trapezoid distance (m)
        R    = 0         ! Hydraulic radius
        WP   = 0         ! wetted perimmeter
         h   = 0        ! depth of flow
        h_0  = 0
        h_1  = 0       ! secant method estimates
        Qj_0 = 0      ! secant method estimates
        D    = 0
        D1   = 0        ! diffusion coeff
        dtr   = 0       ! required timestep, minutes
        error = 1.0
        hp    = 0        !courant, previous height
        maxiter = 0
        a = 0
101     continue
!end yw

        c = 0.52    !-- coefficnets for finding dx/Ckdt
        b = 1.15
        
        if(Cs .eq.0) then 
         z = 1.0 
        else 
         z = 1/Cs              !channel side distance (m)
        endif

        !qC = quc + ql !current upstream in reach

        if (n .le.0 .or. So .le. 0 .or. z .le. 0 .or. Bw .le. 0) then
          print*, "Error in channel coefficients -> Muskingum cunge",n,So,z,Bw
          call hydro_stop("In MUSKINGCUNGE() - Error in channel coefficients")
        end if

        error   = 1.0
        maxiter = 0
        a = 0.0

        if ((quc+ql) .lt. 100) then 
          b=5 
        else 
         b= 20
        endif

!-------------  Secant Method
        h    =  (a+b)/2  !- upper interval
        h_0  = 0.0       !- lower interval
        Qj_0 = 0.0       !- initial flow of lower interval

        do while ((error .gt. 0.05 .and. maxiter .le. 100 .and. h .gt. 0.01))  

          !----- lower interval  --------------------
           Tw = Bw + 2*z*h_0                    !--top width of the channel inflow
           Ck = (sqrt(So)/n)*(5./3.)*(h_0**0.667)   !-- pg 287 Chow, Mdt, Mays
           if(Ck .gt. 0.0) then
             Km = dx/Ck                       !-- seconds Muskingum Param
             if(Km .lt. dt) then           
               Km = dt
             endif
           else 
             Km = dt
           endif
 
           if(Tw*So*Ck*dx .eq. 0.0) then 
             X = 0.25
           else
             X = 0.5-(Qj_0/(2*Tw*So*Ck*dx))
           endif
  
           if(X .le. 0.0) then 
             X = 0.25
           elseif(X .gt. 0.35) then
             X = 0.35
           endif
  
           D = (Km*(1 - X) + dt/2)              !--seconds
            if(D .eq. 0.0) then 
              print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
              call hydro_stop("In MUSKINGCUNGE() - D is 0.")
           endif 
  
           C1 =  (Km*X + dt/2)/D
           C2 =  (dt/2 - Km*X)/D
           C3 =  (Km*(1-X)-dt/2)/D
           C4 =  (ql*dt)/D                       !-- ql already multipled by the dx length

           if(h_0 .le. 0.0) then 
             AREA= 0.0
             WP = 0.0
           else
            AREA = (Bw * h_0 + z * (h_0*h_0) )
            WP   = (Bw * h_0 + z * (h_0*h_0)) / (Bw + 2 * h_0 * sqrt(1+z*z))
           endif
           
           if(WP .le. 0.0) then 
              Qj_0 =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4)
           else 
              Qj_0 =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) - ((1/n) * AREA * (WP**(2./3.)) * sqrt(So)) !f(x)
           endif
           
           !--upper interval -----------
           Tw = Bw + 2*z*h                    !--top width of the channel inflow
           Ck = (sqrt(So)/n)*(5./3.)*(h**0.667)   !-- pg 287 Chow, Mdt, Mays
           if(Ck .gt. 0.0) then
             Km = dx/Ck                       !-- seconds Muskingum Param
             if(Km .lt. dt) then           
               Km = dt
             endif
           else 
             Km = dt
           endif
 
           if(Tw*So*Ck*dx .eq. 0.0) then 
             X = 0.25
           else
             X = 0.5-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2*Tw*So*Ck*dx))
           endif
  
           if(X .le. 0.0) then 
             X = 0.25
           elseif(X .gt. 0.35) then
             X = 0.35
           endif
  
           D = (Km*(1 - X) + dt/2)              !--seconds
            if(D .eq. 0.0) then 
              print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
              call hydro_stop("In MUSKINGCUNGE() - D is 0.")
           endif 
  
           C1 =  (Km*X + dt/2)/D
           C2 =  (dt/2 - Km*X)/D
           C3 =  (Km*(1-X)-dt/2)/D
           C4 =  (ql*dt)/D                       !-- ql already multipled by the dx length

           if(h .le. 0) then 
            AREA = 0.0
            WP   = 0.0
           else
            AREA = (Bw * h + z * (h*h) )
            WP   = (Bw * h + z * (h*h)) / (Bw + 2 * h * sqrt(1+z*z))
           endif 
           
           if(WP .le. 0.0) then 
             Qj =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4)
           else
             Qj =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) -((1/n) * AREA * (WP**(2./3.)) * sqrt(So))
           endif

           if(Qj_0-Qj .ne. 0.0) then
             h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) !update h, 3rd estimate
              if(h_1 .lt. 0.0) then
                h_1 = h
              endif
           else
             h_1 = h
           endif

           error = abs((h_1 - h)/h) !error is new estatimate and 2nd estimate

!           if(idx .eq. 626) then 
!             write(6,*) h_0,h,h_1,error
!           endif

           h_0  = h 
           h    = h_1
           maxiter = maxiter + 1

      end do

      if((maxiter .ge. 100 .and. error .gt. 0.05) .or. h .gt. 100) then 

         print*, "WARNING:"
         print*, "id,err,iters,h", idx, error, maxiter, h
         print*, "n,z,B,So,dx,X,dt,Km",n,z,Bw,So,dx,X,dt,Km
         print*, "qup,quc,qdp,ql", qup,quc,qdp,ql
         if(h.gt.100) then
              print*, "FATAL ERROR: Water Elevation Calculation is Diverging"
              call hydro_stop("In MUSKINGCUNGE() - Water Elevation Calculation is Diverging")
         endif
      endif

!      if(idx .eq. 626) then 
!         write(6,*) ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber
!      endif

!     MUSKINGCUNGE =  h 
!yw      MUSKINGCUNGE =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber

!yw added for test

      if(((C1*qup)+(C2*quc)+(C3*qdp) + C4) .lt. 0.0) then
         MUSKINGCUNGE =  MAX( ( (C1*qup)+(C2*quc) + C4),((C1*qup)+(C3*qdp) + C4) )
      else
          MUSKINGCUNGE =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber
      endif


! ----------------------------------------------------------------

  END FUNCTION MUSKINGCUNGE


  SUBROUTINE SUBMUSKINGCUNGE(qdc,vel,idx,qup,quc,qdp,ql,dt,So,dx,n,Cs,Bw)

        IMPLICIT NONE

        REAL, intent(IN)       :: dt         !routing period in  seconds
        REAL, intent(IN)       :: qup        !flow upstream previous timestep
        REAL, intent(IN)       :: quc        !flow upstream current timestep
        REAL, intent(IN)       :: qdp        !flow downstream previous timestep
        REAL, intent(INOUT)    :: qdc        !flow downstream current timestep
        REAL, intent(IN)       :: ql         !lateral inflow through reach (m^3/sec)
        REAL, intent(IN)       :: Bw         ! bottom width (meters)
        REAL, intent(IN)       :: Cs         ! Channel side slope slope
        REAL, intent(IN)       :: So         ! Channel bottom slope %
        REAL, intent(IN)       :: dx         ! channel lngth (m)
        REAL, intent(IN)       :: n          ! mannings coefficient
        REAL, intent(INOUT)    :: vel        ! mannings coefficient
        INTEGER, intent(IN)    :: idx        ! channel id

!--local variables
        REAL    :: C1, C2, C3, C4
        REAL    :: Km          !K travel time in hrs in reach
        REAL    :: X           !weighting factors 0<=X<=0.5
        REAL    :: Ck          ! wave celerity (m/s)

!-- channel geometry and characteristics
        REAL    :: Tw         ! top width at peak flow
        REAL    :: AREA       ! Cross sectional area m^2
        REAL    :: Z          ! trapezoid distance (m)
        REAL    :: R          ! Hydraulic radius
        REAL    :: WP         ! wetted perimmeter
        REAL    :: h          ! depth of flow
        REAL    :: h_0,h_1    ! secant method estimates
        REAL    :: Qj_0       ! secant method estimates
        REAL    :: Qj         ! intermediate flow estimate
        REAL    :: D,D1       ! diffusion coeff
        REAL    :: dtr        ! required timestep, minutes
        REAL    :: error
        REAL    :: hp         !courant, previous height
        INTEGER :: maxiter    !maximum number of iterations

!-- local variables.. needed if channel is sub-divded
        REAL    :: a,b,c
        INTEGER :: i          !-- channel segment counter

!yw add
        goto 101
        C1   = 0
        C2   = 0
        C3   = 0
        C4   = 0
        Km   = 0
        X    = 0        !weighting factors 0<=X<=0.5
        Ck   = 0
        Tw   = 0       ! top width at peak flow
        AREA = 0      ! Cross sectional area m^2
        Z    = 0         ! trapezoid distance (m)
        R    = 0         ! Hydraulic radius
        WP   = 0         ! wetted perimmeter
         h   = 0        ! depth of flow
        h_0  = 0
        h_1  = 0       ! secant method estimates
        Qj_0 = 0      ! secant method estimates
        D    = 0
        D1   = 0        ! diffusion coeff
        dtr   = 0       ! required timestep, minutes
        error = 1.0
        hp    = 0        !courant, previous height
        maxiter = 0
        a = 0
101     continue
!end yw

        c = 0.52    !-- coefficnets for finding dx/Ckdt
        b = 1.15
        
        if(Cs .eq.0) then 
         z = 1.0 
        else 
         z = 1/Cs              !channel side distance (m)
        endif

        !qC = quc + ql !current upstream in reach

        if (n .le.0 .or. So .le. 0 .or. z .le. 0 .or. Bw .le. 0) then
          print*, "Error in channel coefficients -> Muskingum cunge",n,So,z,Bw
          call hydro_stop("In MUSKINGCUNGE() - Error in channel coefficients")
        end if

        error   = 1.0
        maxiter = 0
        a = 0.0

        if ((quc+ql) .lt. 100) then 
          b=5 
        else 
         b= 20
        endif

!-------------  Secant Method
        h    =  (a+b)/2  !- upper interval
        h_0  = 0.0       !- lower interval
        Qj_0 = 0.0       !- initial flow of lower interval

        do while ((error .gt. 0.05 .and. maxiter .le. 100 .and. h .gt. 0.01))  

          !----- lower interval  --------------------
           Tw = Bw + 2*z*h_0                    !--top width of the channel inflow
           Ck = (sqrt(So)/n)*(5./3.)*(h_0**0.667)   !-- pg 287 Chow, Mdt, Mays
           if(Ck .gt. 0.0) then
             Km = dx/Ck                       !-- seconds Muskingum Param
             if(Km .lt. dt) then           
               Km = dt
             endif
           else 
             Km = dt
           endif
 
           if(Tw*So*Ck*dx .eq. 0.0) then 
             X = 0.25
           else
             X = 0.5-(Qj_0/(2*Tw*So*Ck*dx))
           endif
  
           if(X .le. 0.0) then 
             X = 0.25
           elseif(X .gt. 0.35) then
             X = 0.35
           endif
  
           D = (Km*(1 - X) + dt/2)              !--seconds
            if(D .eq. 0.0) then 
              print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
              call hydro_stop("In MUSKINGCUNGE() - D is 0.")
           endif 
  
           C1 =  (Km*X + dt/2)/D
           C2 =  (dt/2 - Km*X)/D
           C3 =  (Km*(1-X)-dt/2)/D
           C4 =  (ql*dt)/D                       !-- ql already multipled by the dx length

           if(h_0 .le. 0.0) then 
             AREA= 0.0
             WP = 0.0
           else
            AREA = (Bw * h_0 + z * (h_0*h_0) )
            WP   = (Bw * h_0 + z * (h_0*h_0)) / (Bw + 2 * h_0 * sqrt(1+z*z))
           endif
           
           if(WP .le. 0.0) then 
              Qj_0 =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4)
           else 
              Qj_0 =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) - ((1/n) * AREA * (WP**(2./3.)) * sqrt(So)) !f(x)
           endif
           
           !--upper interval -----------
           Tw = Bw + 2*z*h                    !--top width of the channel inflow
           Ck = (sqrt(So)/n)*(5./3.)*(h**0.667)   !-- pg 287 Chow, Mdt, Mays
           if(Ck .gt. 0.0) then
             Km = dx/Ck                       !-- seconds Muskingum Param
             if(Km .lt. dt) then           
               Km = dt
             endif
           else 
             Km = dt
           endif
 
           if(Tw*So*Ck*dx .eq. 0.0) then 
             X = 0.25
           else
             X = 0.5-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2*Tw*So*Ck*dx))
           endif
  
           if(X .le. 0.0) then 
             X = 0.25
           elseif(X .gt. 0.35) then
             X = 0.35
           endif
  
           D = (Km*(1 - X) + dt/2)              !--seconds
            if(D .eq. 0.0) then 
              print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
              call hydro_stop("In MUSKINGCUNGE() - D is 0.")
           endif 
  
           C1 =  (Km*X + dt/2)/D
           C2 =  (dt/2 - Km*X)/D
           C3 =  (Km*(1-X)-dt/2)/D
           C4 =  (ql*dt)/D                       !-- ql already multipled by the dx length

           if(h .le. 0) then 
            AREA = 0.0
            WP   = 0.0
           else
            AREA = (Bw * h + z * (h*h) )
            WP   = (Bw * h + z * (h*h)) / (Bw + 2 * h * sqrt(1+z*z))
           endif 
           
           if(WP .le. 0.0) then 
             Qj =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4)
           else
             Qj =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) -((1/n) * AREA * (WP**(2./3.)) * sqrt(So))
           endif

           if(Qj_0-Qj .ne. 0.0) then
             h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) !update h, 3rd estimate
              if(h_1 .lt. 0.0) then
                h_1 = h
              endif
           else
             h_1 = h
           endif

           error = abs((h_1 - h)/h) !error is new estatimate and 2nd estimate

!           if(idx .eq. 626) then 
!             write(6,*) h_0,h,h_1,error
!           endif

           h_0  = h 
           h    = h_1
           maxiter = maxiter + 1

      end do

      if((maxiter .ge. 100 .and. error .gt. 0.05) .or. h .gt. 100) then 

         print*, "WARNING:"
         print*, "id,err,iters,h", idx, error, maxiter, h
         print*, "n,z,B,So,dx,X,dt,Km",n,z,Bw,So,dx,X,dt,Km
         print*, "qup,quc,qdp,ql", qup,quc,qdp,ql
         if(h.gt.100) then
              print*, "FATAL ERROR: Water Elevation Calculation is Diverging"
              call hydro_stop("In MUSKINGCUNGE() - Water Elevation Calculation is Diverging")
         endif
      endif

!yw added for test
      if(((C1*qup)+(C2*quc)+(C3*qdp) + C4) .lt. 0.0) then
!       MUSKINGCUNGE =  MAX( ( (C1*qup)+(C2*quc) + C4),((C1*qup)+(C3*qdp) + C4) )
        qdc = MAX( ( (C1*qup)+(C2*quc) + C4),((C1*qup)+(C3*qdp) + C4) )

      else
!       MUSKINGCUNGE =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber
        qdc =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber

      endif

      Tw = Bw + (2*z*h)
      R = (h*(Bw + Tw) / 2) / (Bw + 2*(((Tw - Bw) / 2)**2 + h**2)**0.5)    
      vel =  (1./n) * (R **(2./3.)) * sqrt(So)  ! average velocity in m/s

! ----------------------------------------------------------------
!END FUNCTION MUSKINGCUNGE
END SUBROUTINE SUBMUSKINGCUNGE
! ----------------------------------------------------------------

! ------------------------------------------------
!   FUNCTION KINEMATIC
! ------------------------------------------------
	REAL FUNCTION KINEMATIC()

	IMPLICIT NONE

! -------- DECLARATIONS -----------------------
 
!	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)	:: OVRGH

        KINEMATIC = 1       
!----------------------------------------------------------------
  END FUNCTION KINEMATIC
!----------------------------------------------------------------


! ------------------------------------------------
!   SUBROUTINE drive_CHANNEL
! ------------------------------------------------
! ------------------------------------------------
     Subroutine drive_CHANNEL(latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, &
       QSUBRT, LAKEINFLORT, QSTRMVOLRT, TO_NODE, FROM_NODE, &
       TYPEL, ORDER, MAXORDER, NLINKS, CH_NETLNK, CH_NETRT, CH_LNKRT, &
       LAKE_MSKRT, DT, DTCT, DTRT_CH,MUSK, MUSX, QLINK, &
       HLINK, ELRT, CHANLEN, MannN, So, ChSSlp, Bw, &
       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  &
#ifdef MPP_LAND 
       , lake_index,link_location,mpp_nlinks,nlinks_index,yw_mpp_nlinks  &
       , LNLINKSL, LLINKID  &
       , gtoNode,toNodeInd,nToNodeInd &
#endif
       , CH_LNKRT_SL &
       ,gwBaseSwCRT, gwHead, qgw_chanrt, gwChanCondSw, gwChanCondConstIn, &
       gwChanCondConstOut)


       IMPLICIT NONE

! -------- DECLARATIONS ------------------------

        INTEGER, INTENT(IN) :: IXRT,JXRT,channel_option
        INTEGER, INTENT(IN) :: NLINKS,NLAKES, NLINKSL
        integer, INTENT(INOUT) :: KT   ! flag of cold start (1) or continue run.
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QSUBRT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QSTRMVOLRT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: LAKEINFLORT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: ELRT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QINFLOWBASE
        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETLNK

        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT
        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT_SL

       real , dimension(ixrt,jxrt):: latval,lonval

        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
        INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: ORDER, TYPEL !--link
        INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: TO_NODE, FROM_NODE
        INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: CHANXI, CHANYJ
        REAL,    INTENT(IN), DIMENSION(NLINKS)    :: ZELEV  !--elevation of nodes
        REAL, INTENT(INOUT), DIMENSION(NLINKS)    :: CVOL
        REAL, INTENT(IN), DIMENSION(NLINKS)       :: MUSK, MUSX
        REAL, INTENT(IN), DIMENSION(NLINKS)       :: CHANLEN
        REAL, INTENT(IN), DIMENSION(NLINKS)       :: So, MannN
        REAL, INTENT(IN), DIMENSION(NLINKS)       :: ChSSlp,Bw  !--properties of nodes or links
        REAL                                      :: Km, X
        REAL , INTENT(INOUT), DIMENSION(:,:) :: QLINK
        REAL ,  DIMENSION(NLINKS,2) :: tmpQLINK
        REAL , INTENT(INOUT), DIMENSION(NLINKS)   :: HLINK
        REAL, INTENT(IN)                          :: DT    !-- model timestep
        REAL, INTENT(IN)                          :: DTRT_CH  !-- routing timestep
        REAL, INTENT(INOUT)                       :: DTCT
        real                                      :: minDTCT !BF minimum routing timestep
        REAL                                      :: dist(ixrt,jxrt,9)
        REAL                                      :: RETDEP_CHAN
        INTEGER, INTENT(IN)                       :: MAXORDER, SUBRTSWCRT, &
                                                     gwBaseSwCRT, gwChanCondSw
        real, intent(in)                          :: gwChanCondConstIn, gwChanCondConstOut ! aquifer-channel conductivity constant from namelist                                             
        REAL , INTENT(IN), DIMENSION(NLINKS)   :: node_area

!DJG GW-chan coupling variables...
        REAL, DIMENSION(NLINKS)                   :: dzGwChanHead
        REAL, DIMENSION(NLINKS)                   :: Q_GW_CHAN_FLUX     !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state...
        REAL, DIMENSION(IXRT,JXRT)                :: ZWATTBLRT          !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state...
        REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: gwHead            !DJG !!! groundwater head from Fersch-2d gw implementation...units (m ASL)
        REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: qgw_chanrt         !DJG !!! Channel-gw flux as used in Fersch 2d gw implementation...units (m^3/s)...Change 'INTENT' to 'OUT' when ready to update groundwater state...
         


        !-- lake params
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: HRZAREA  !-- horizontal area (km^2)
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: LAKEMAXH !-- maximum lake depth  (m^2)
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: WEIRH    !-- lake depth  (m^2)
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: WEIRC    !-- weir coefficient
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: WEIRL    !-- weir length (m)
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: ORIFICEC !-- orrifice coefficient
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: ORIFICEA !-- orrifice area (m^2)
        REAL, INTENT(IN), DIMENSION(NLAKES)       :: ORIFICEE !-- orrifce elevation (m)

        REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: RESHT    !-- reservoir height (m)
        REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: QLAKEI   !-- lake inflow (cms)
        REAL,                DIMENSION(NLAKES)    :: QLAKEIP  !-- lake inflow previous timestep (cms)
        REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: QLAKEO   !-- outflow from lake used in diffusion scheme
        INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: LAKENODE !-- outflow from lake used in diffusion scheme
        INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: LINKID   !--  id of channel elements for linked scheme
        REAL, DIMENSION(NLINKS)                   :: QLateral !--lateral flow
        REAL, DIMENSION(NLINKS)                   :: QSUM     !--mass bal of node
        REAL, DIMENSION(NLAKES)                   :: QLLAKE   !-- lateral inflow to lake in diffusion scheme

!-- Local Variables
        INTEGER                     :: i,j,k,t,m,jj,kk,KRT,node
        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
        REAL                        :: AREA,WP                !--wetted area and perimiter for MuskingC. routing

        REAL ,DIMENSION(NLINKS)     :: HLINKTMP,CVOLTMP       !-- temporarily store head values and volume values
        REAL ,DIMENSION(NLINKS)     :: CD                     !-- critical depth
        real, DIMENSION(IXRT,JXRT)  :: tmp
        real, dimension(nlinks)     :: tmp2

#ifdef MPP_LAND
        integer lake_index(nlakes)
        integer nlinks_index(nlinks)
        integer mpp_nlinks, iyw, yw_mpp_nlinks
        integer link_location(ixrt,jxrt)
        real     ywtmp(ixrt,jxrt)
        integer LNLINKSL
        INTEGER, dimension(LNLINKSL) :: LLINKID
        real*8,  dimension(LNLINKSL) :: LQLateral
!        real*4,  dimension(LNLINKSL) :: LQLateral
        integer, dimension(:) ::  toNodeInd
        integer, dimension(:,:) ::  gtoNode
        integer  :: nToNodeInd
        real, dimension(nToNodeInd,2) :: gQLINK
#else
        REAL*8, DIMENSION(NLINKS)                   :: LQLateral !--lateral flow
#endif
        integer flag

        integer :: n, kk2, nt, nsteps  ! tmp 

        QLAKEIP = 0
        HLINKTMP = 0
        CVOLTMP = 0
        CD = 0  
        node = 1
        QLateral = 0
        QSUM     = 0
        QLLAKE   = 0


!yw      print *, "DRIVE_channel,option,nlinkl,nlinks!!", channel_option,NLINKSL,NLINKS
         
      dzGwChanHead = 0.

   IF(channel_option .ne. 3) then   !--muskingum methods ROUTE ON DT timestep, not DTRT!!

         nsteps = (DT+0.5)/DTRT_CH

#ifdef MPP_LAND
         LQLateral = 0          !-- initial lateral flow to 0 for this reach
         DO iyw = 1,yw_MPP_NLINKS
         jj = nlinks_index(iyw)
          !--------river grid points, convert depth in mm to rate across reach in m^3/sec
              if( .not. (  (CHANXI(jj) .eq. 1 .and. left_id .ge. 0) .or. &
                           (CHANXI(jj) .eq. ixrt .and. right_id .ge. 0) .or. &
                           (CHANYJ(jj) .eq. 1 .and. down_id .ge. 0) .or. &
                           (CHANYJ(jj) .eq. jxrt .and. up_id .ge. 0)      &
                   ) ) then
                  if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0) then
                     k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
                     LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 & 
                            *node_area(jj)/DT)
                  elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
                      k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
                      LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
                               *node_area(jj)/DT)
                  endif
              endif
         end do  ! jj


!   assign LQLATERAL to QLATERAL
       call updateLinkV(LQLateral, QLateral(1:NLINKSL))

#else
         LQLateral = 0          !-- initial lateral flow to 0 for this reach
         do jj = 1, NLINKS
          !--------river grid points, convert depth in mm to rate across reach in m^3/sec

                  if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0 ) then
                     k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
                     LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 & 
                            *node_area(jj)/DT)
                  elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
                      k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
                      LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
                               *node_area(jj)/DT)
                  endif   

          end do  ! jj
          QLateral = LQLateral
#endif

!       QLateral = QLateral / nsteps

   do nt = 1, nsteps

 
!----------  route order 1 reaches which have no upstream inflow
        do k=1, NLINKSL
           if (ORDER(k) .eq. 1) then  !-- first order stream has no headflow


              if(TYPEL(k) .eq. 1) then    !-- level pool route of reservoir
                  !CALL LEVELPOOL(1,0.0, 0.0, qd, QLINK(k,2), QLateral(k), &
                  ! DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), &
                  ! WEIRC(k), WEIRL(k), ORIFICEE(i), ORIFICEC(k), ORIFICEA(k) )
              elseif (channel_option .eq. 1) then
                   Km  = MUSK(k)
                   X   = MUSX(k)
                   QLINK(k,2) = MUSKING(k,0.0, QLateral(k), QLINK(k,1), DTRT_CH, Km, X) !--current outflow
              elseif (channel_option .eq. 2) then !-- upstream is assumed constant initial condition

                 ! HLINK(k) = MUSKINGCUNGE(k,0.0,0.0,QLINK(k,1), &
                   QLINK(k,2) = MUSKINGCUNGE(k,0.0,0.0,QLINK(k,1), &
                         QLateral(k), DTRT_CH, So(k), CHANLEN(k), &
                         MannN(k), ChSSlp(k), Bw(k))

                 ! AREA = (Bw(k) * HLINK(k) + 1/ChSSlp(k) * HLINK(k)**2)
                 ! WP   = (Bw(k) * HLINK(k) + 1/ChSSlp(k) * HLINK(k)**2) / (Bw(k) + 2 * HLINK(k) * sqrt(1+(1/ChSSlp(k))**2))
                 ! QLINK(k,2) = 1/MannN(k) * AREA * WP**(2./3.) * sqrt(So(k))

              else
                  print *, "FATAL ERROR: No channel option selected"
                  call hydro_stop("In drive_CHANNEL() -No channel option selected ") 
              endif
           endif
        end do

#ifdef MPP_LAND
       gQLINK = 0
       call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
       call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1))
#endif

      !---------- route other reaches, with upstream inflow
       tmpQlink = 0
       do k = 1,NLINKSL
          if (ORDER(k) .gt. 1 ) then  !-- exclude first order stream 
             Quc  = 0
             Qup  = 0

#ifdef MPP_LAND
!using mapping index
               do n = 1, gtoNODE(k,1)
                  m = gtoNODE(k,n+1)
!yw                  if (LINKID(k) .eq. m) then
                    Quc = Quc + gQLINK(m,2)  !--accum of upstream inflow of current timestep (2)
                    Qup = Qup + gQLINK(m,1)  !--accum of upstream inflow of previous timestep (1)

                      !     if(LINKID(k) .eq. 3259 .or. LINKID(k) .eq. 3316 .or. LINKID(k) .eq. 3219) then
                      !       write(6,*) "id,Uc,Up",LINKID(k),Quc,Qup
                      !       call flush(6)
                      !     endif

!yw                  endif
                end do ! do i

#else
               do m = 1, NLINKSL
                  if (LINKID(k) .eq. TO_NODE(m)) then
                    Quc = Quc + QLINK(m,2)  !--accum of upstream inflow of current timestep (2)
                    Qup = Qup + QLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
                  endif
               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
                       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

                    !HLINK(k) =   MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), &
                    tmpQLINK(k,2) = MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), &
                                    QLateral(k), DTRT_CH, So(k),  CHANLEN(k), &
                                    MannN(k), ChSSlp(k), Bw(k) )

                    ! AREA = (Bw(k) * HLINK(k) + 1/ChSSLP(k) * HLINK(k)**2)
                    ! WP   = (Bw(k) * HLINK(k) + 1/ChSSLP(k) * HLINK(k)**2) / (Bw(k) + 2 * HLINK(k) * sqrt(1+(1/ChSSLP(k))**2))
                    ! tmpQLINK(k,2) = ((1/MannN(k)) * AREA * WP**(2./3.) * sqrt(So(k)))

                   else
                    print *, "FATAL ERROR: no channel option selected"
                    call hydro_stop("In drive_CHANNEL() - no channel option selected") 
                   endif
            endif !!! order(1) .ne. 1
         end do       !--k links

!yw check
!        gQLINK = 0.0
!        call ReachLS_write_io(tmpQLINK(:,2), gQLINK(:,2))
!        call ReachLS_write_io(tmpQLINK(:,1), gQLINK(:,1))
!        write(6,*) " io_id = ", io_id
!        if(my_id .eq. io_id) then
!            write(71,*) gQLINK(:,1)
!            call flush(71)
!            call flush(72)
!        endif

          do k = 1, NLINKSL
            if(TYPEL(k) .ne. 1) 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
         end do

!#ifdef MPP_LAND 
!         call ReachLS_write_io(QLINK(:,2),buf1)
!         if(my_id .eq. IO_id) write(73,*) buf1
!#else
!         write(73,*) QLINK(1:NLINKSL,2)
!#endif

#ifdef HYDRO_D
          print *, "END OF ALL REACHES...",KRT,DT_STEPS
#endif

   end do  ! nsteps

!    END DO !-- krt timestep for muksingumcunge routing

   elseif(channel_option .eq. 3) then   !--- route using the diffusion scheme on nodes not links

#ifdef MPP_LAND
         call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINK,NLINKS,99)
         call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOL,NLINKS,99)
#endif

         KRT = 0                  !-- initialize the time counter
         minDTCT = 0.01           ! define minimum routing sub-timestep (s), simulation will end with smaller timestep
         DTCT = min(max(DTCT*2.0, minDTCT),DTRT_CH)
       
         HLINKTMP = HLINK         !-- temporary storage of the water elevations (m)
         CVOLTMP = CVOL           !-- temporary storage of the volume of water in channel (m^3)
         QLAKEIP = QLAKEI         !-- temporary lake inflow from previous timestep  (cms)

!        call check_channel(77,HLINKTMP,1,nlinks)
!        call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,ZELEV,NLINKS,99)
 crnt:    DO                      !-- loop on the courant condition
         QSUM   = 0              !-- initialize the total flow out of each cell to zero
         QLAKEI = 0              !-- set the lake inflow as zero
         QLLAKE = 0              !-- initialize each lake's lateral inflow to zero  
         DT_STEPS=INT(DT/DTCT)   !-- fix the timestep
         QLateral = 0. 
!DJG GW-chan coupling variables...
         if(gwBaseSwCRT == 3) then
	  Q_GW_CHAN_FLUX = 0.
	  qgw_chanrt     = 0.
         end if
         
!         ZWATTBLRT=1.0   !--HARDWIRE, remove this and pass in from subsfc/gw routing routines...


!-- vectorize
!--------------------- 
#ifdef MPP_LAND
         DO iyw = 1,yw_MPP_NLINKS
         i = nlinks_index(iyw)
#else
         DO i = 1,NLINKS
#endif
          
           if(node_area(i) .eq. 0) then
               write(6,*) "FATAL ERROR: node_area(i) is zero. i=", i
               call hydro_stop("In drive_CHANNEL() - Error node_area") 
           endif

           

nodeType:if((CH_NETRT(CHANXI(i), CHANYJ(i) ) .eq. 0) .and. &
              (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .lt.0) ) then !--a reg. node
              
gwOption:   if(gwBaseSwCRT == 3) then

             ! determine potential gradient between groundwater head and channel stage
             ! units in (m)
             dzGwChanHead(i) = gwHead(CHANXI(i),CHANYJ(i)) - (HLINK(i)+ZELEV(i)) 

             if(gwChanCondSw .eq. 0) then
	       
                qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
	        
             else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) > 0) then
	       
	       ! channel bed interface, units in (m^3/s), flux into channel...
	       ! BF todo: consider channel width
                qgw_chanrt(CHANXI(i),CHANYJ(i)) = gwChanCondConstIn * dzGwChanHead(i) &
                                                * CHANLEN(i) * 2. 

             else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) < 0) then
	       
	       ! channel bed interface, units in (m^3/s), flux out of channel...
	       ! BF todo: consider channel width
                qgw_chanrt(CHANXI(i),CHANYJ(i)) = max(-0.005, gwChanCondConstOut * dzGwChanHead(i) &
                                                * CHANLEN(i) * 2.)
!              else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then  TBD: exponential dependency
!              else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then
	       
             else
	       
	        qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
	        
             end if
             
             Q_GW_CHAN_FLUX(i) = qgw_chanrt(CHANXI(i),CHANYJ(i))
!             if ( i .eq. 1001 ) then
!                print *, Q_GW_CHAN_FLUX(i), dzGwChanHead(i), ELRT(CHANXI(i),CHANYJ(i)), HLINK(i), ZELEV(i)
!             end if
!              if ( Q_GW_CHAN_FLUX(i) .lt. 0. ) then   !-- temporary hardwire for only allowing flux into channel...REMOVE later...
!                 Q_GW_CHAN_FLUX(i) = 0.
! 	        qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
!              end if
            
            else
	      Q_GW_CHAN_FLUX(i) = 0.
	    end if gwOption


              QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) =  &
!DJG  awaiting gw-channel exchg...  Q_GW_CHAN_FLUX(i)+& ...obsolete-> ((QSUBRT(CHANXI(i),CHANYJ(i))+&
                Q_GW_CHAN_FLUX(i)+&
                ((QSTRMVOLRT(CHANXI(i),CHANYJ(i))+&
                 QINFLOWBASE(CHANXI(i),CHANYJ(i))) &
                   /DT_STEPS*node_area(i)/1000/DTCT)
	if((QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))).lt.0.) .and. (gwChanCondSw == 0)) then
#ifdef HYDRO_D
               print*, "i, CHANXI(i),CHANYJ(i) = ", i, CHANXI(i),CHANYJ(i)
               print *, "NEGATIVE Lat inflow...",QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))), &
                         QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)), &
                         QINFLOWBASE(CHANXI(i),CHANYJ(i))
#endif
        end if
            elseif(LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .gt. 0 .and. &
               (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .ne. -9999)) then !--a lake node
              QLLAKE(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) = &
                 QLLAKE(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) + &
                 (LAKEINFLORT(CHANXI(i),CHANYJ(i))+ &
                 QINFLOWBASE(CHANXI(i),CHANYJ(i)) &
                 /DT_STEPS*node_area(i)/1000/DTCT)
            elseif(CH_NETRT(CHANXI(i),CHANYJ(i)) .gt. 0) then  !pour out of lake
                 QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) =  &
                   QLAKEO(CH_NETRT(CHANXI(i),CHANYJ(i)))  !-- previous timestep
          endif nodeType
        ENDDO


#ifdef MPP_LAND
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLateral,NLINKS,99)
    if(NLAKES .gt. 0) then
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT   ,ixrt,jxrt,QLLAKE,NLAKES,99)
    endif
#endif

          !-- compute conveyances, with known depths (just assign to QLINK(,1)
          !--QLINK(,2) will not be used), QLINK is the flow across the node face
          !-- units should be m3/second.. consistent with QL (lateral flow)

#ifdef MPP_LAND
         DO iyw = 1,yw_MPP_NLINKS
         i = nlinks_index(iyw)
#else
           DO i = 1,NLINKS
#endif
           if (TYPEL(i) .eq. 0 .AND. HLINKTMP(FROM_NODE(i)) .gt. RETDEP_CHAN) then 
               if(from_node(i) .ne. to_node(i) .and. (to_node(i) .gt. 0) .and.(from_node(i) .gt. 0) ) &  ! added by Wei Yu
                   QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
                     HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
                     CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
            else !--  we are just computing critical depth for outflow points
               QLINK(i,1) =0.
            endif
          ENDDO

#ifdef MPP_LAND
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
#endif
 

          !-- compute total flow across face, into node
#ifdef MPP_LAND
         DO iyw = 1,yw_mpp_nlinks
         i = nlinks_index(iyw)
#else
          DO i = 1,NLINKS                                                 !-- inflow to node across each face
#endif
           if(TYPEL(i) .eq. 0) then                                       !-- only regular nodes have to attribute
              QSUM(TO_NODE(i)) = QSUM(TO_NODE(i)) + QLINK(i,1)
           endif
          END DO

#ifdef MPP_LAND
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,qsum,NLINKS,0)
#endif



#ifdef MPP_LAND
         DO iyw = 1,yw_mpp_nlinks
         i = nlinks_index(iyw)
#else
          DO i = 1,NLINKS                                                 !-- outflow from node across each face
#endif
            QSUM(FROM_NODE(i)) = QSUM(FROM_NODE(i)) - QLINK(i,1)
          END DO
#ifdef MPP_LAND
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,qsum,NLINKS,99)
#endif


         flag = 99


#ifdef MPP_LAND
         DO iyw = 1,yw_MPP_NLINKS
             i = nlinks_index(iyw)
#else
          DO i = 1, NLINKS                                                !--- compute volume and depth at each node
#endif
 
           if( TYPEL(i).eq.0 .and. CVOLTMP(i) .ge. 0.001 .and.(CVOLTMP(i)-QSUM(i)*DTCT)/CVOLTMP(i) .le. -0.01 ) then  
            flag = -99
#ifdef HYDRO_D
            write(6,*) "******* start diag ***************"
            write(6,*) "Unstable at node ",i, "i=",CHANXI(i),"j=",CHANYJ(i)
            write(6,*) "Unstatble at node ",i, "lat=",latval(CHANXI(i),CHANYJ(i)), "lon=",lonval(CHANXI(i),CHANYJ(i))
            write(6,*) "TYPEL, CVOLTMP, QSUM, QSUM*DTCT",TYPEL(i), CVOLTMP(i), QSUM(i), QSUM(i)*DTCT
            write(6,*) "qsubrt, qstrmvolrt,qlink",QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)),qlink(i,1),qlink(i,2)
!              write(6,*) "current nodes, z, h", ZELEV(FROM_NODE(i)),HLINKTMP(FROM_NODE(i))
!           if(TO_NODE(i) .gt. 0) then
!              write(6,*) "to nodes, z, h", ZELEV(TO_NODE(i)), HLINKTMP(TO_NODE(i))
!           else
!              write(6,*) "no to nodes   "
!           endif
               write(6,*) "CHANLEN(i), MannN(i), Bw(i), ChSSlp(i) ", CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)
            write(6,*) "*******end of  diag ***************"
#endif
            
            goto 999  
            endif 
          enddo 

999 continue
#ifdef MPP_LAND
        call mpp_same_int1(flag)
#endif


        if(flag < 0  .and. DTCT >0.1)   then   
             
             ! call smoth121(HLINK,nlinks,maxv_p,pnode,to_node)

             if(DTCT .gt. minDTCT) then                !-- timestep in seconds
              DTCT = max(DTCT/2 , minDTCT)             !-- 1/2 timestep
              KRT = 0                                  !-- restart counter
              HLINKTMP = HLINK                         !-- set head and vol to start value of timestep
              CVOLTMP = CVOL
              CYCLE crnt                               !-- start cycle over with smaller timestep
             else
              write(6,*) "Courant error with smallest routing timestep DTCT: ",DTCT
!              call hydro_stop("drive_CHANNEL")
              DTCT = 0.1
              HLINKTMP = HLINK                          !-- set head and volume to start values of timestep
              CVOLTMP  = CVOL
              goto 998  
             end if
        endif 

998 continue


#ifdef MPP_LAND
         DO iyw = 1,yw_MPP_NLINKS
            i = nlinks_index(iyw)
#else
         DO i = 1, NLINKS                                                !--- compute volume and depth at each node
#endif
 
            if(TYPEL(i) .eq. 0) then                   !--  regular channel grid point, compute volume
              CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) )* DTCT
              if((CVOLTMP(i) .lt. 0) .and. (gwChanCondSw == 0)) then 
#ifdef HYDRO_D
                print *, "WARNING! channel volume less than 0:i,CVOL,QSUM,QLat", &
                               i, CVOLTMP(i),QSUM(i),QLateral(i),HLINK(i)
#endif
                CVOLTMP(i) =0 
              endif

            elseif(TYPEL(i) .eq. 1) then               !-- pour point, critical depth downstream 

              if (QSUM(i)+QLateral(i) .lt. 0) then
              else

!DJG remove to have const. flux b.c....   CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
                  CD(i) = HLINKTMP(i)  !This is a temp hardwire for flow depth for the pour point...
               endif

               ! change in volume is inflow, lateral flow, and outflow 
               !yw DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*DXRT),HLINKTMP(i), &
                   CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
                       DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
                       CD(i),CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
            elseif (TYPEL(i) .eq. 2) then              !--- into a reservoir, assume critical depth
              if ((QSUM(i)+QLateral(i) .lt. 0) .and. (gwChanCondSw == 0)) then
#ifdef HYDRO_D
               print *, i, 'CrtDpth Qsum+QLat into lake< 0',QSUM(i), QLateral(i)
#endif
              else
!DJG remove to have const. flux b.c....    CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
               CD(i) = HLINKTMP(i)  !This is a temp hardwire for flow depth for the pour point...
              endif
 
              !-- compute volume in reach (m^3)
                   CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
                          DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
                             CD(i) ,CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
              !-- compute flow rate into lake from all contributing nodes (cms)
              QLAKEI(LAKENODE(i)) = QLAKEI(LAKENODE(i)) + QLINK(FROM_NODE(i),1)

            else
              print *, "FATAL ERROR: This node does not have a type.. error TYPEL =", TYPEL(i)
              call hydro_stop("In drive_CHANNEL() - error TYPEL") 
            endif
           
           if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
            HLINKTMP(i) = HEAD(i, CVOLTMP(i)/CHANLEN(i),Bw(i),1/ChSSlp(i))  !--updated depth 
           else
            HLINKTMP(i) = CD(i)  !!!   CRITICALDEPTH(i,QSUM(i)+QLateral(i), Bw(i), 1./ChSSlp(i)) !--critical depth is head
           endif 

           END DO  !--- done processing all the links


#ifdef MPP_LAND
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOLTMP,NLINKS,99)
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CD,NLINKS,99)
    if(NLAKES .gt. 0) then
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
    endif
    call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINKTMP,NLINKS,99)
#endif 
!   call check_channel(83,CVOLTMP,1,nlinks)
!   call check_channel(84,CD,1,nlinks)
!   call check_channel(85,HLINKTMP,1,nlinks)
!   call check_lake(86,QLAKEI,lake_index,nlakes)

      



           do i = 1, NLAKES !-- mass balances of lakes
#ifdef MPP_LAND
            if(lake_index(i) .gt. 0)  then
#endif
              CALL LEVELPOOL(i,QLAKEIP(i), QLAKEI(i), QLAKEO(i), QLLAKE(i), &
                DTCT, RESHT(i), HRZAREA(i), WEIRH(i), LAKEMAXH(i), WEIRC(i), &
                WEIRL(i), ORIFICEE(i), ORIFICEC(i), ORIFICEA(i))
                QLAKEIP(i) = QLAKEI(i)  !-- store total lake inflow for this timestep
#ifdef MPP_LAND
            endif
#endif
           enddo
#ifdef MPP_LAND
    if(NLAKES .gt. 0) then
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLLAKE,NLAKES,99)
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,RESHT,NLAKES,99)
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEO,NLAKES,99)
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEIP,NLAKES,99)
    endif
#endif


#ifdef MPP_LAND
         DO iyw = 1,yw_MPP_NLINKS
            i = nlinks_index(iyw)
#else
         DO i = 1, NLINKS                                                !--- compute volume and depth at each node
#endif
            if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
                   QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
                      HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
                      CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
            endif
         enddo

#ifdef MPP_LAND
          call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
#endif

           KRT = KRT + 1                     !-- iterate on the timestep
           IF(KRT .eq. DT_STEPS) EXIT crnt   !-- up to the maximum time in interval

          END DO crnt  !--- DTCT timestep of DT_STEPS
 
           HLINK = HLINKTMP                 !-- update head based on final solution in timestep
           CVOL  = CVOLTMP                  !-- update volume
        else                                !-- no channel option apparently selected
         print *, "FATAL ERROR: no channel option selected"
         call hydro_stop("In drive_CHANNEL() - no channel option selected") 
        endif

#ifdef HYDRO_D
         write(6,*) "finished call drive_CHANNEL"
#endif

        if (KT .eq. 1) KT = KT + 1
         

 END SUBROUTINE drive_CHANNEL
! ----------------------------------------------------------------

!-=======================================
     REAL FUNCTION AREAf(AREA,Bw,h,z)
     REAL :: AREA, Bw, z, h
       AREAf = (Bw+z*h)*h-AREA       !-- Flow area
     END FUNCTION AREAf

!-====critical depth function  ==========
     REAL FUNCTION CDf(Q,Bw,h,z)
     REAL :: Q, Bw, z, h
       if(h .le. 0) then
         print *, "FATAL ERROR: head is zero, will get division by zero error"
         call hydro_stop("In CDf() - head is zero") 
       else
       CDf = (Q/((Bw+z*h)*h))/(sqrt(9.81*(((Bw+z*h)*h)/(Bw+2*z*h))))-1  !--critical depth function
       endif
     END FUNCTION CDf

!=======find flow depth in channel with bisection Chapra pg. 131
    REAL FUNCTION HEAD(idx,AREA,Bw,z)  !-- find the water elevation given wetted area, 
                                         !--bottom widith and side channel.. index was for debuggin
     REAL :: Bw,z,AREA,test           
     REAL :: hl, hu, hr, hrold
     REAL :: fl, fr,error                !-- function evaluation
     INTEGER :: maxiter, idx

     error = 1.0
     maxiter = 0
     hl = 0.00001   !-- minimum depth is small
     hu = 30.  !-- assume maximum depth is 30 meters

    if (AREA .lt. 0.00001) then 
     hr = 0.
    else
      do while ((AREAf(AREA,BW,hl,z)*AREAf(AREA,BW,hu,z)).gt.0 .and. maxiter .lt. 100) 
       !-- allows for larger , smaller heads 
       if(AREA .lt. 1.) then
        hl=hl/2
       else
        hu = hu * 2
       endif
       maxiter = maxiter + 1
        
      end do

      maxiter =0
      hr = 0
      fl = AREAf(AREA,Bw,hl,z)
      do while (error .gt. 0.0001 .and. maxiter < 1000)
        hrold = hr
        hr = (hl+hu)/2
        fr =  AREAf(AREA,Bw,hr,z)
        maxiter = maxiter + 1
         if (hr .ne. 0) then
          error = abs((hr - hrold)/hr)
         endif
        test = fl * fr
         if (test.lt.0) then
           hu = hr
         elseif (test.gt.0) then
           hl=hr
           fl = fr
         else
           error = 0.0
         endif
      end do
     endif
     HEAD = hr

22   format("i,hl,hu,Area",i5,2x,f12.8,2x,f6.3,2x,f6.3,2x,f6.3,2x,f9.1,2x,i5)

    END FUNCTION HEAD
!=================================
     REAL FUNCTION MANNING(h1,n,Bw,Cs)

     REAL :: Bw,h1,Cs,n
     REAL :: z, AREA,R,Kd

     z=1/Cs
     R = ((Bw+z*h1)*h1)/(Bw+2*h1*sqrt(1+z*z)) !-- Hyd Radius
     AREA = (Bw+z*h1)*h1        !-- Flow area
     Kd = (1/n)*(R**(2./3.))*AREA     !-- convenyance
#ifdef HYDRO_D
     print *,"head, kd",  h1,Kd
#endif
     MANNING = Kd
     
     END FUNCTION MANNING

!=======find flow depth in channel with bisection Chapra pg. 131
     REAL FUNCTION CRITICALDEPTH(lnk,Q,Bw,z)  !-- find the critical depth
     REAL :: Bw,z,Q,test
     REAL :: hl, hu, hr, hrold
     REAL :: fl, fr,error   !-- function evaluation
     INTEGER :: maxiter
     INTEGER :: lnk

     error = 1.0
     maxiter = 0
     hl = 1e-5   !-- minimum depth is 0.00001 meters
!    hu = 35.       !-- assume maximum  critical depth 25 m
     hu = 100.       !-- assume maximum  critical depth 25 m

     if(CDf(Q,BW,hl,z)*CDf(Q,BW,hu,z) .gt. 0) then
      if(Q .gt. 0.001) then
#ifdef HYDRO_D
        print *, "interval won't work to find CD of lnk ", lnk
        print *, "Q, hl, hu", Q, hl, hu
        print *, "cd lwr, upr", CDf(Q,BW,hl,z), CDf(Q,BW,hu,z)
        ! call hydro_stop("In CRITICALDEPTH()") 
        CRITICALDEPTH = -9999
        return
#endif
      else
        Q = 0.0
      endif
     endif

     hr = 0.
     fl = CDf(Q,Bw,hl,z)

     if (Q .eq. 0.) then
       hr = 0.
     else
      do while (error .gt. 0.0001 .and. maxiter < 1000)
        hrold = hr
        hr = (hl+hu)/2
        fr =  CDf(Q,Bw,hr,z)
        maxiter = maxiter + 1
         if (hr .ne. 0) then
          error = abs((hr - hrold)/hr)
         endif
        test = fl * fr
         if (test.lt.0) then
           hu = hr
         elseif (test.gt.0) then
           hl=hr
           fl = fr
         else
           error = 0.0
         endif

       end do
      endif

     CRITICALDEPTH = hr

     END FUNCTION CRITICALDEPTH
!================================================
     REAL FUNCTION SGNf(val)  !-- function to return the sign of a number
     REAL:: val

     if (val .lt. 0) then
       SGNf= -1.
     elseif (val.gt.0) then
       SGNf= 1.
     else
       SGNf= 0.
     endif

     END FUNCTION SGNf
!================================================

     REAL FUNCTION fnDX(qp,Tw,So,Ck,dx,dt) !-- find channel sub-length for MK method
     REAL    :: qp,Tw,So,Ck,dx, dt,test
     REAL    :: dxl, dxu, dxr, dxrold
     REAL    :: fl, fr, error
     REAL    :: X
     INTEGER :: maxiter

     error = 1.0
     maxiter =0
     dxl = dx*0.9  !-- how to choose dxl???
     dxu = dx
     dxr=0

     do while (fnDXCDT(qp,Tw,So,Ck,dxl,dt)*fnDXCDT(qp,Tw,So,Ck,dxu,dt) .gt. 0 &
               .and. dxl .gt. 10)  !-- don't let dxl get too small
      dxl = dxl/1.1
     end do
     
      
     fl = fnDXCDT(qp,Tw,So,Ck,dxl,dt)
     do while (error .gt. 0.0001 .and. maxiter < 1000)
        dxrold = dxr
        dxr = (dxl+dxu)/2
        fr =  fnDXCDT(qp,Tw,So,Ck,dxr,dt)
        maxiter = maxiter + 1
         if (dxr .ne. 0) then
          error = abs((dxr - dxrold)/dxr)
         endif
        test = fl * fr
         if (test.lt.0) then
           dxu = dxr
         elseif (test.gt.0) then
           dxl=dxr
           fl = fr
         else
           error = 0.0
         endif
      end do
     FnDX = dxr

    END FUNCTION fnDX
!================================================
     REAL FUNCTION fnDXCDT(qp,Tw,So,Ck,dx,dt) !-- function to help find sub-length for MK method
     REAL    :: qp,Tw,So,Ck,dx,dt,X
     REAL    :: c,b  !-- coefficients on dx/cdt log approximation function
     
     c = 0.2407
     b = 1.16065
     X = 0.5-(qp/(2*Tw*So*Ck*dx))
     if (X .le.0) then 
      fnDXCDT = -1 !0.115
     else
      fnDXCDT = (dx/(Ck*dt)) - (c*LOG(X)+b)  !-- this function needs to converge to 0
     endif
     END FUNCTION fnDXCDT
! ----------------------------------------------------------------------

    subroutine check_lake(unit,cd,lake_index,nlakes)
         use module_RT_data, only: rt_domain
         implicit none 
         integer :: unit,nlakes,i,lake_index(nlakes)
         real cd(nlakes)
#ifdef MPP_LAND
         call write_lake_real(cd,lake_index,nlakes)
#endif
         write(unit,*) cd
          call flush(unit)
         return
    end subroutine check_lake

    subroutine check_channel(unit,cd,did,nlinks)
         use module_RT_data, only: rt_domain
#ifdef MPP_LAND
  USE module_mpp_land
#endif
         implicit none 
         integer :: unit,nlinks,i, did
         real cd(nlinks)
#ifdef MPP_LAND
         real g_cd(rt_domain(did)%gnlinks)
         call write_chanel_real(cd,rt_domain(did)%map_l2g,rt_domain(did)%gnlinks,nlinks,g_cd)
         if(my_id .eq. IO_id) then
            write(unit,*) "rt_domain(did)%gnlinks = ",rt_domain(did)%gnlinks
           write(unit,*) g_cd
         endif
#else
           write(unit,*) cd
#endif
          call flush(unit)
          close(unit)
         return
    end subroutine check_channel
    subroutine smoth121(var,nlinks,maxv_p,from_node,to_node)
        implicit none
        integer,intent(in) ::  nlinks, maxv_p
        integer, intent(in), dimension(nlinks):: to_node
        integer, intent(in), dimension(nlinks):: from_node(nlinks,maxv_p)
        real, intent(inout), dimension(nlinks) :: var
        real, dimension(nlinks) :: vartmp
        integer :: i,j  , k, from,to
        integer :: plen
              vartmp = 0
              do i = 1, nlinks
                 to = to_node(i)
                 plen = from_node(i,1)
                 if(plen .gt. 1) then 
                     do k = 1, plen-1 
                         from = from_node(i,k+1)
                         if(to .gt. 0) then
                            vartmp(i) = vartmp(i)+0.25*(var(from)+2.*var(i)+var(to))
                         else
                            vartmp(i) = vartmp(i)+(2.*var(i)+var(from))/3.0
                         endif
                     end do
                     vartmp(i) = vartmp(i) /(plen-1)
                 else
                         if(to .gt. 0) then
                            vartmp(i) = vartmp(i)+(2.*var(i)+var(to)/3.0)
                         else
                            vartmp(i) = var(i)
                         endif
                 endif
              end do
              var = vartmp 
        return
    end subroutine smoth121

!   SUBROUTINE drive_CHANNEL for NHDPLUS
! ------------------------------------------------

     Subroutine drive_CHANNEL_RSL(UDMP_OPT,KT, IXRT,JXRT,  &
        LAKEINFLORT, QSTRMVOLRT, TO_NODE, FROM_NODE, &
        TYPEL, ORDER, MAXORDER,   CH_LNKRT, &
        LAKE_MSKRT, DT, DTCT, DTRT_CH,MUSK, MUSX, QLINK, &
        CHANLEN, MannN, So, ChSSlp, Bw, &
        RESHT, HRZAREA, LAKEMAXH, WEIRH, WEIRC, WEIRL, ORIFICEC, ORIFICEA, &
        ORIFICEE,  CVOL, QLAKEI, QLAKEO, LAKENODE, &
        QINFLOWBASE, CHANXI, CHANYJ, channel_option,  &
        nlinks,NLINKSL, LINKID, node_area, qout_gwsubbas, &
        LAKEIDA, LAKEIDM, NLAKES, LAKEIDX, &
#ifdef MPP_LAND 
        nlinks_index,mpp_nlinks,yw_mpp_nlinks, &
        LNLINKSL, &
        gtoNode,toNodeInd,nToNodeInd,   &
#endif
         CH_LNKRT_SL, landRunOff  & 
#ifdef WRF_HYDRO_NUDGING
       , nudge &
#endif

       ,  accLndRunOff, accQLateral, accStrmvolrt, accBucket   &
       , QLateral, velocity &
       ,nsize , OVRTSWCRT, SUBRTSWCRT      )

       use module_UDMAP, only: LNUMRSL, LUDRSL

#ifdef WRF_HYDRO_NUDGING
       use module_stream_nudging,  only: setup_stream_nudging,  & 
                                         nudge_term_all,        &
                                         nudgeWAdvance
#endif 


       IMPLICIT NONE

! -------- DECLARATIONS ------------------------

        INTEGER, INTENT(IN) :: IXRT,JXRT,channel_option, OVRTSWCRT, SUBRTSWCRT
        INTEGER, INTENT(IN) :: NLAKES, NLINKSL, nlinks
        integer, INTENT(INOUT) :: KT   ! flag of cold start (1) or continue run.
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QSTRMVOLRT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: LAKEINFLORT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QINFLOWBASE
        real, dimension(ixrt,jxrt) :: landRunOff

        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT
        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT_SL

        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
        INTEGER, INTENT(IN), DIMENSION(:)         :: ORDER, TYPEL !--link
        INTEGER, INTENT(IN), DIMENSION(:)     :: TO_NODE, FROM_NODE
        INTEGER, INTENT(IN), DIMENSION(:)     :: CHANXI, CHANYJ
        REAL, INTENT(IN), DIMENSION(:)        :: MUSK, MUSX
        REAL, INTENT(IN), DIMENSION(:)        :: CHANLEN
        REAL, INTENT(IN), DIMENSION(:)        :: So, MannN
        REAL, INTENT(IN), DIMENSION(:)        :: ChSSlp,Bw  !--properties of nodes or links
        REAL                                      :: Km, X
        REAL , INTENT(INOUT), DIMENSION(:,:)  :: QLINK
#ifdef WRF_HYDRO_NUDGING
        real, intent(out),    dimension(:)    :: nudge
#endif
        REAL, DIMENSION(:), intent(out)       :: QLateral, velocity !--lateral flow
        real, dimension(:), intent(out)       :: accLndRunOff, accQLateral, accStrmvolrt, accBucket 

        REAL ,  DIMENSION(NLINKSL,2) :: tmpQLINK
        REAL, INTENT(IN)                          :: DT    !-- model timestep
        REAL, INTENT(IN)                          :: DTRT_CH  !-- routing timestep
        REAL, INTENT(INOUT)                       :: DTCT
        real                                      :: minDTCT !BF minimum routing timestep
        INTEGER, INTENT(IN)                       :: MAXORDER
        REAL , INTENT(IN), DIMENSION(:)   :: node_area

!DJG GW-chan coupling variables...
        REAL, DIMENSION(NLINKS)                   :: dzGwChanHead
        REAL, DIMENSION(NLINKS)                   :: Q_GW_CHAN_FLUX     !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state...
        REAL, DIMENSION(IXRT,JXRT)                :: ZWATTBLRT          !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state...

        !-- lake params

        REAL, INTENT(IN), DIMENSION(:)       :: HRZAREA  !-- horizontal area (km^2)
        REAL, INTENT(IN), DIMENSION(:)       :: LAKEMAXH !-- maximum lake depth  (m^2)
        REAL, INTENT(IN), DIMENSION(:)       :: WEIRH    !--  lake depth  (m^2)
        REAL, INTENT(IN), DIMENSION(:)       :: WEIRC    !-- weir coefficient
        REAL, INTENT(IN), DIMENSION(:)       :: WEIRL    !-- weir length (m)
        REAL, INTENT(IN), DIMENSION(:)       :: ORIFICEC !-- orrifice coefficient
        REAL, INTENT(IN), DIMENSION(:)       :: ORIFICEA !-- orrifice area (m^2)
        REAL, INTENT(IN), DIMENSION(:)       :: ORIFICEE !-- orrifce elevation (m)
        INTEGER, INTENT(IN), DIMENSION(:)    :: LAKEIDM  !-- NHDPLUS lakeid for lakes to be modeled

        REAL, INTENT(INOUT), DIMENSION(:)    :: RESHT    !-- reservoir height (m)
        REAL, INTENT(INOUT), DIMENSION(:)    :: QLAKEI   !-- lake inflow (cms)
        REAL,                DIMENSION(NLAKES)    :: QLAKEIP  !-- lake inflow previous timestep (cms)
        REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: QLAKEO   !-- outflow from lake used in diffusion scheme

        INTEGER, INTENT(IN), DIMENSION(:)    :: LAKENODE !-- outflow from lake used in diffusion scheme
        INTEGER, INTENT(IN), DIMENSION(:)   :: LINKID   !--  id of channel elements for linked scheme
        INTEGER, INTENT(IN), DIMENSION(:)   :: LAKEIDA  !--  (don't need) NHDPLUS lakeid for all lakes in domain
        INTEGER, INTENT(IN), DIMENSION(:)   :: LAKEIDX  !--  the sequential index of the lakes id by com id

        REAL, DIMENSION(NLINKS)                   :: QSUM     !--mass bal of node
        REAL, DIMENSION(NLAKES)                   :: QLLAKE   !-- lateral inflow to lake in diffusion scheme
        integer :: nsize

!-- Local Variables
        INTEGER                      :: i,j,k,t,m,jj,ii,lakeid, kk,KRT,node, UDMP_OPT
        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

        REAL ,DIMENSION(NLINKS)                          :: CD    !-- critical depth
        real, DIMENSION(IXRT,JXRT)                       :: tmp
        real, dimension(nlinks)                          :: tmp2
        REAL, INTENT(INOUT), DIMENSION(:)           :: CVOL

#ifdef MPP_LAND
        real*8,  dimension(LNLINKSL) :: LQLateral
        real*8,  dimension(LNLINKSL) :: tmpLQLateral
        real,  dimension(NLINKSL)    :: tmpQLateral
        integer nlinks_index(:)
        integer  iyw, yw_mpp_nlinks, mpp_nlinks
        real     ywtmp(ixrt,jxrt)
        integer LNLINKSL
        integer, dimension(:)         ::  toNodeInd
        integer, dimension(:,:)       ::  gtoNode
        integer  :: nToNodeInd
        real, dimension(nToNodeInd,2) :: gQLINK
#else
        real*8,  dimension(NLINKS) :: tmpLQLateral
        real,  dimension(NLINKSL) :: tmpQLateral
        real,  dimension(NLINKSL) :: LQLateral
#endif
        integer flag

        integer :: n, kk2, nt, nsteps  ! tmp 
        real, dimension(:) :: qout_gwsubbas
        real, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI, tmpRESHT
  

        real, dimension(NLINKS) ::  lcLndRunOff, lcQLateral, lcStrmvolrt, lcBucket  ! local variables




#ifdef MPP_LAND
        if(my_id .eq. io_id) then
#endif
            allocate(tmpQLAKEO(NLAKES))
            allocate(tmpQLAKEI(NLAKES))
            allocate(tmpRESHT(NLAKES))
#ifdef MPP_LAND
        endif
#endif

        QLAKEIP = 0
        CD = 0  
        node = 1
        QLateral = 0
        QSUM     = 0
        QLLAKE   = 0
        dzGwChanHead = 0.

#ifdef WRF_HYDRO_NUDGING
         !! Initialize nudging for the current timestep.
         !! This establishes the data structure used to solve the nudges. 
         call setup_stream_nudging(0)  !! always zero b/c at beginning of hydro timestep
#endif /* WRF_HYDRO_NUDGING */

         nsteps = (DT+0.5)/DTRT_CH
         LQLateral = 0          !-- initial lateral flow to 0 for this reach


        tmpLQLateral = 0
        tmpQLateral = 0

        ! NHDPLUS maping
    if(OVRTSWCRT .eq. 0)      then
        do k = 1, LNUMRSL
           ! get from land grid runoff
             do m = 1, LUDRSL(k)%ncell  
                ii =  LUDRSL(k)%cell_i(m)
                jj =  LUDRSL(k)%cell_j(m)
                LQLateral(k) = LQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 & 
                              *LUDRSL(k)%cellArea(m)/DT
                tmpLQLateral(k) = tmpLQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 & 
                              *LUDRSL(k)%cellArea(m)/DT
             end do
        end do
#ifdef MPP_LAND
        call updateLinkV(tmpLQLateral, tmpQLateral)
#endif
        if(NLINKSL .gt. 0) then
            accLndRunOff(1:NLINKSL) = accLndRunOff(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
        endif
        tmpLQLateral = 0
        tmpQLateral = 0
     endif

     if(OVRTSWCRT .ne. 0 .or. SUBRTSWCRT .ne. 0 ) then
        do k = 1, LNUMRSL
              ! get from channel grid
              do m = 1, LUDRSL(k)%ngrids
                  ii =  LUDRSL(k)%grid_i(m)
                  jj =  LUDRSL(k)%grid_j(m)
                  LQLateral(k) = LQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 & 
                          *LUDRSL(k)%nodeArea(m)/DT
                  tmpLQLateral(k) = tmpLQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 & 
                          *LUDRSL(k)%nodeArea(m)/DT
              end do
        end do
#ifdef MPP_LAND
        call updateLinkV(tmpLQLateral, tmpQLateral)
#endif
        if(NLINKSL .gt. 0) then
            accStrmvolrt(1:NLINKSL) = accStrmvolrt(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
        endif
     endif


#ifdef MPP_LAND
       call updateLinkV(LQLateral, QLateral(1:NLINKSL))
#else
       call hydro_stop("fatal error: NHDPlus only works for parallel now.")
       QLateral = LQLateral
#endif

     if(NLINKSL .gt. 0) then
        QLateral(1:NLINKSL) = QLateral(1:NLINKSL) + qout_gwsubbas(1:NLINKSL)
     endif

     ! accQLateral  = accLndRunOff + QLateral * DT 
     if(NLINKSL .gt. 0) then
        accQLateral(1:NLINKSL)  = accQLateral(1:NLINKSL) + QLateral(1:NLINKSL) * DT 
        accBucket(1:NLINKSL) = accBucket(1:NLINKSL) + qout_gwsubbas(1:NLINKSL) * DT
     endif

!       QLateral = QLateral / nsteps

   do nt = 1, nsteps
 
#ifdef MPP_LAND

       gQLINK = 0
       call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
       call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1)) 
      !---------- route other reaches, with upstream inflow
#endif

       tmpQlink = 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

        Quc  = 0
        Qup  = 0

        !process as standard link or a lake inflow link, or lake outflow link
        ! link flowing out of lake, accumulate all the inflows with the revised TO_NODEs
        ! TYPEL = -999 stnd; TYPEL=1 outflow from lake; TYPEL = 3 inflow to a lake

        if(TYPEL(k) .ne. 2) then ! don't process internal lake links only

#ifdef MPP_LAND
!using mapping index
           do n = 1, gtoNODE(k,1)
              m = gtoNODE(k,n+1)
                if(gQLINK(m,2) .gt. 0)   Quc = Quc + gQLINK(m,2)  !--accum of upstream inflow of current timestep (2)
                if(gQLINK(m,1) .gt. 0)   Qup = Qup + gQLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
           end do ! do i
#else
           do m = 1, NLINKSL

               if (LINKID(k) .eq. TO_NODE(m)) then
                 Quc = Quc + QLINK(m,2)  !--accum of upstream inflow of current timestep (2)
                 Qup = Qup + QLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
               endif
           end do ! do m
#endif
        endif !note that we won't process type 2 links, since they are internal to a lake


!yw ### process each link k,
!       There is a situation that different k point to the same LAKEIDX
!        if(TYPEL(k) .eq. 1 .and. LAKEIDX(k) .gt. 0) then   !--link is a reservoir
        if(TYPEL(k) .eq. 1 ) then   !--link is a reservoir
             
             lakeid = LAKEIDX(k)
           if(lakeid .ge. 0) then
             CALL LEVELPOOL(lakeid,Qup, Quc, tmpQLINK(k,2), &
               QLateral(k), DT, RESHT(lakeid), HRZAREA(lakeid), WEIRH(lakeid), LAKEMAXH(lakeid), &
               WEIRC(lakeid), WEIRL(lakeid),ORIFICEE(lakeid), ORIFICEC(lakeid), ORIFICEA(lakeid))
             
               QLAKEO(lakeid)  = tmpQLINK(k,2) !save outflow to lake
               QLAKEI(lakeid)  = Quc           !save inflow to lake
           endif
105  continue
            

        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, don't process internal lake nodes TYP=2
!              tmpQLINK(k,2) = MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), &
!                  QLateral(k), DTRT_CH, So(k),  CHANLEN(k), &
!                  MannN(k), ChSSlp(k), Bw(k) )

              CALL SUBMUSKINGCUNGE(tmpQLINK(k,2),velocity(k), k,Qup, Quc, QLINK(k,1), &
                   QLateral(k), DTRT_CH, So(k), CHANLEN(k), &
                   MannN(k), ChSSlp(k), Bw(k) )
                
        else
#ifdef HYDRO_D
                    print *, " no channel option selected"
#endif
                    call hydro_stop("drive_CHANNEL") 
        endif

       END DO        !--k links

#ifdef MPP_LAND
       call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO)
       call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI)
       call updateLake_seq(RESHT,nlakes,tmpRESHT)
#endif

       do k = 1, NLINKSL !tmpQLINK?
          if(TYPEL(k) .ne. 2) then !only the internal lake nodes don't have info.. but need to save QLINK of lake out too
             QLINK(k,2) = tmpQLINK(k,2)
          endif
            QLINK(k,1) = QLINK(k,2)    !assigng link flow of current to be previous for next time step
       end do


#ifdef WRF_HYDRO_NUDGING         
         if(.not. nudgeWAdvance) call nudge_term_all(qlink, nudge, int(nt*dtrt_ch))
#endif /* WRF_HYDRO_NUDGING */


!#ifdef HYDRO_D
!          print *, "END OF ALL REACHES...",KRT,DT_STEPS
!#endif

    end do  ! nsteps

    if (KT .eq. 1) KT = KT + 1

#ifdef MPP_LAND
    if(my_id .eq. io_id)      then 
       if(allocated(tmpQLAKEO))  deallocate(tmpQLAKEO)
       if(allocated(tmpQLAKEI))  deallocate(tmpQLAKEI)
       if(allocated(tmpRESHT))  deallocate(tmpRESHT)
    endif
#endif        

   if (KT .eq. 1) KT = KT + 1

 END SUBROUTINE drive_CHANNEL_RSL

! ----------------------------------------------------------------

END MODULE module_channel_routing

#ifdef MPP_LAND
 subroutine checkReach(ii,  inVar)
   use module_mpp_land
   use module_RT_data, only: rt_domain
   use MODULE_mpp_ReachLS, only : updatelinkv,                   &
                                 ReachLS_write_io, gbcastvalue, &
                                 gbcastreal2
   implicit none
   integer :: ii
   real,dimension(rt_domain(1)%nlinksl) :: inVar
   real:: g_var(rt_domain(1)%gnlinksl)
   call ReachLS_write_io(inVar, g_var)
   if(my_id .eq. io_id) then
      write(ii,*) g_var
      call flush(ii)
   endif
 end subroutine checkReach
#endif