hydro/Routing/module_channel_routing.F90

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 MODULE module_channel_routing
  22 #ifdef MPP_LAND
  23 use module_mpp_land
  24 use MODULE_mpp_ReachLS, only : updatelinkv,                   &
  25                                ReachLS_write_io, gbcastvalue, &
  26                                gbcastreal2,      linkls_s
  27
  28
  29 use module_reservoir, only: reservoir
  30 use module_RT_data, only: rt_domain
  31 use module_hydro_stop, only: HYDRO_stop
  32
  33 use iso_fortran_env, only: int64
  34
  35 #endif
  36 implicit none
  37
  38 contains
  39
  40 ! ------------------------------------------------
  41 !   FUNCTION MUSKING
  42 ! ------------------------------------------------
  43 real function MUSKING(idx,qup,quc,qdp,dt,Km,X)
  44
  45 implicit none
  46
  47 !--local variables
  48 real    :: C1, C2, C3
  49 real    :: Km         !K travel time in hrs in reach
  50 real    :: X          !weighting factors 0<=X<=0.5
  51 real    :: dt         !routing period in hrs
  52 real    :: avgbf      !average base flow for initial condition
  53 real    :: qup        !inflow from previous timestep
  54 real    :: quc        !inflow  of current timestep
  55 real    :: qdp        !outflow of previous timestep
  56 real    :: dth        !timestep in hours
  57 integer :: idx        ! index
  58
  59 dth = dt/3600.    !hours in timestep
  60 C1 = (dth - 2.0 *Km*X)/(2.0*Km*(1.0-X)+dth)
  61 C2 = (dth+2.0*Km*X)/(2.0*Km*(1.0-X)+dth)
  62 C3 = (2.0*Km*(1.0-X)-dth)/(2.0*Km*(1.0-X)+dth)
  63 MUSKING = (C1*quc)+(C2*qup)+(C3*qdp)
  64
  65 ! ----------------------------------------------------------------
  66 end function MUSKING
  67 ! ----------------------------------------------------------------
  68
  69
  70 ! ------------------------------------------------
  71 !   FUNCTION Diffusive wave
  72 ! ------------------------------------------------
  73 real function DIFFUSION(nod,z1,z20,h1,h2,dx,n, &
  74      Bw, Cs)
  75 implicit none
  76 !-- channel geometry and characteristics
  77 real    :: Bw         !-bottom width (meters)
  78 real    :: Cs         !-Channel side slope slope
  79 real    :: dx         !-channel lngth (m)
  80 real,intent(in)    :: n          !-mannings coefficient
  81 real    :: R          !-Hydraulic radius
  82 real    :: AREA       !- wetted area
  83 real    :: h1,h2      !-tmp height variables
  84 real    :: z1,z2      !-z1 is 'from', z2 is 'to' elevations
  85 real    :: z          !-channel side distance
  86 real    :: w          !-upstream weight
  87 real    :: Ku,Kd      !-upstream and downstream conveyance
  88 real    :: Kf         !-final face conveyance
  89 real    :: Sf         !-friction slope
  90 real    :: sgn        !-0 or 1
  91 integer :: nod         !- node
  92 real ::  z20, dzx
  93
  94 ! added by Wei Yu for bad data.
  95
  96 dzx = (z1 - z20)/dx
  97 if(dzx .lt. 0.002) then
  98    z2 = z1 - dx*0.002
  99 else
 100    z2 = z20
 101 endif
 102 !end
 103
 104 if (n.le.0.0.or.Cs.le.0.or.Bw.le.0) then
 105    print *, "Error in Diffusion function ->channel coefficients"
 106    print *, "nod, n, Cs, Bw", nod, n, Cs, Bw
 107    call hydro_stop("In DIFFUSION() - Error channel coefficients.")
 108 endif
 109
 110 !        Sf = ((z1+h1)-(z2+h2))/dx  !-- compute the friction slope
 111 !if(z1 .eq. z2) then
 112 ! Sf = ((z1-(z2-0.01))+(h1-h2))/dx  !-- compute the friction slope
 113 !else
 114 !         Sf = ((z1-z2)+(h1-h2))/dx  !-- compute the friction slope
 115 !endif
 116
 117 !modifieed by Wei Yu for false geography data
 118 if(abs(z1-z2) .gt. 1.0E5) then
 119 #ifdef HYDRO_D
 120    print*, "WARNING: huge slope rest to 0 for channel grid.", z1,z2
 121 #endif
 122    Sf = ((h1-h2))/dx  !-- compute the friction slope
 123 else
 124    Sf = ((z1-z2)+(h1-h2))/dx  !-- compute the friction slope
 125 endif
 126 !end  modfication
 127
 128 sgn = SGNf(Sf)             !-- establish sign
 129
 130 w = 0.5*(sgn + 1.)         !-- compute upstream or downstream weighting
 131
 132 z = 1.0/Cs                   !--channel side distance (m)
 133 R = ((Bw + z*h1)*h1) / (Bw+ 2.0*h1*sqrt(1.0 + z*z)) !-- Hyd Radius
 134 AREA = (Bw + z*h1)*h1        !-- Flow area
 135 Ku = (1.0/n)*(R**(2./3.))*AREA     !-- convenyance
 136
 137 R = ((Bw + z*h2)*h2)/(Bw + 2.0*h2*sqrt(1.0 + z*z)) !-- Hyd Radius
 138 AREA = (Bw + z*h2)*h2        !-- Flow area
 139 Kd = (1.0/n)*(R**(2.0/3.0))*AREA     !-- convenyance
 140
 141 Kf =  (1.0-w)*Kd + w*Ku      !-- conveyance
 142 DIFFUSION = Kf * sqrt(abs(Sf))*sgn
 143
 144
 145 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)
 146
 147 end function DIFFUSION
 148 ! ----------------------------------------------------------------
 149
 150 subroutine SUBMUSKINGCUNGE(    &
 151      qdc, vel, qloss, idx, qup,  quc, &
 152      qdp, ql,   dt,  So,   dx, &
 153      n,   Cs,   Bw,  Tw, TwCC, &
 154      nCC, depth, ChannK        )
 155
 156 #ifdef HYDRO_D
 157 use module_RT_data, only: rt_domain  !! JLM: this is only used in a c3 paramter diagnostic print
 158 #endif
 159
 160         IMPLICIT NONE
 161
 162         REAL, intent(IN)       :: dt         ! routing period in  seconds
 163         REAL, intent(IN)       :: qup        ! flow upstream previous timestep
 164         REAL, intent(IN)       :: quc        ! flow upstream current timestep
 165         REAL, intent(IN)       :: qdp        ! flow downstream previous timestep
 166         REAL, intent(INOUT)    :: qdc        ! flow downstream current timestep
 167         REAL, intent(IN)       :: ql         ! lateral inflow through reach (m^3/sec)
 168         REAL, intent(IN)       :: Bw         ! bottom width (meters)
 169         REAL, intent(IN)       :: Tw         ! top width before bankfull (meters)
 170         REAL, intent(IN)       :: TwCC       ! top width of Compund (meters)
 171         REAL, intent(IN)       :: nCC        ! mannings of compund
 172         REAL, intent(IN)       :: ChannK     ! Channel conductivity (m/s)
 173         REAL, intent(IN)       :: Cs         ! Channel side slope slope
 174         REAL, intent(IN)       :: So         ! Channel bottom slope %
 175         REAL, intent(IN)       :: dx         ! channel lngth (m)
 176         REAL, intent(IN)       :: n          ! mannings coefficient
 177         REAL, intent(INOUT)    :: vel        ! velocity (m/s)
 178         REAL, intent(INOUT)    :: qloss      ! channel infiltration (m^3/sec)
 179         integer(kind=int64), intent(IN)    :: idx        ! channel id
 180         REAL, intent(INOUT)    :: depth      ! depth of flow in channel
 181
 182 !--local variables
 183         REAL    :: C1, C2, C3, C4
 184         REAL    :: Km             !K travel time in hrs in reach
 185         REAL    :: X              !weighting factors 0<=X<=0.5
 186         REAL    :: Ck             ! wave celerity (m/s)
 187
 188 !-- channel geometry and characteristics, local variables
 189         REAL    :: Twl            ! top width at simulated flow (m)
 190         REAL    :: AREA,AREAC     ! Cross sectional area channel and compound m^2
 191         REAL    :: Z              ! trapezoid distance (m)
 192         REAL    :: R,RC           ! Hydraulic radius of channel and compound
 193         REAL    :: WP,WPC         ! wetted perimmeter of channel and compound
 194         REAL    :: h              ! depth of flow in channel
 195         REAL    :: h_0,h_1        ! secant method estimates
 196         REAL    :: bfd            ! bankfull depth (m)
 197         REAL    :: Qj_0           ! secant method estimates
 198         REAL    :: Qj             ! intermediate flow estimate
 199         REAL    :: D,D1           ! diffusion coeff
 200         REAL    :: dtr            ! required timestep, minutes
 201         REAL    :: aerror,rerror  ! absolute and relative error
 202         REAL    :: hp             ! courant, previous height
 203         INTEGER :: iter, maxiter  ! maximum number of iterations
 204
 205 !-- additional channel infiltration parameters
 206         REAL    :: WPk        ! KINEROS2 modified wetted perimeter
 207         REAL    :: modK       ! KINEROS2 constant (set to 0.15, line )
 208
 209 !-- local variables.. needed if channel is sub-divded
 210         REAL    :: a,b,c,F
 211         REAL    :: mindepth   !  minimum depth in channel
 212         INTEGER :: i,tries    !-- channel segment counter
 213
 214         !TML-- Hard code KINEROS2 Modification Parameter
 215         !TML-- Set to 0.15, consistent with USDA-ARS
 216         modK = 0.15
 217
 218         c = 0.52           !-- coefficnets for finding dx/Ckdt
 219         b = 1.15
 220         a = 0.0
 221         maxiter  = 100
 222         mindepth = 0.01
 223         aerror = 0.01
 224         rerror = 1.0
 225         tries = 0
 226
 227 !-------------  locals
 228         if(Cs .eq. 0.00000000) then
 229          z = 1.0
 230         else
 231          z = 1.0/Cs          !channel side distance (m)
 232         endif
 233
 234         if(Bw .gt. Tw) then   !effectively infinite deep bankful
 235            bfd = Bw/0.00001
 236         elseif (Bw .eq. Tw) then
 237           bfd =  Bw/(2.0*z)  !bankfull depth is effectively
 238         else
 239           bfd =  (Tw - Bw)/(2.0*z)  !bankfull depth (m)
 240         endif
 241         !qC = quc + ql !current upstream in reach
 242
 243         if (n .le. 0.0 .or. So .le. 0.0 .or. z .le. 0.0 .or. Bw .le. 0.0) then
 244           print*, "Error in channel coefficients -> Muskingum cunge", n, So, z, Bw
 245           call hydro_stop("In MUSKINGCUNGE() - Error in channel coefficients")
 246         end if
 247
 248 !-------------  Secant Method
 249         depth = max(depth, 0.0)
 250         h     = (depth * 1.33) + mindepth !1.50 of  depth
 251         h_0   = (depth * 0.67)            !0.50 of depth
 252
 253         if(ql .gt. 0.0 .or. qup .gt. 0.0 .or. qdp .gt. 0.0) then  !only solve if there's water to flux
 254
 255 110     continue
 256
 257         Qj_0  = 0.0                       !- initial flow of lower interval
 258         iter   = 0
 259
 260         do while (rerror .gt. 0.01 .and. aerror .ge. mindepth .and. iter .le. maxiter)
 261
 262            WPC    = 0.0
 263            AREAC  = 0.0
 264
 265           !----- lower interval  --------------------
 266
 267            Twl = Bw + 2.0*z*h_0      !--top surface water width of the channel inflow
 268
 269             if ( (h_0 .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) ) then !water outside of defined channel
 270              AREA =  (Bw + bfd * z) * bfd
 271              AREAC = (TwCC * (h_0 -bfd)) !assume compound component is rect. chan, that's 3 times the Tw
 272              WP = (Bw + 2.0 * bfd * sqrt(1.0 + z*z))
 273              WPC = TwCC + (2.0 * (h_0-bfd)) !WPC is 2 times the Tw
 274              WPk = WP + WPC*MIN((h_0/(modK*SQRT(Bw))),1.0)  ! KINEROS2 Mod.
 275              R   = (AREA + AREAC)/(WP +WPC)  ! hydraulic radius
 276             else
 277               AREA = (Bw + h_0 * z ) * h_0
 278               WP = (Bw + 2.0 * h_0 * sqrt(1.0 + z*z))
 279               WPk = WP*MIN((h_0/(modK*SQRT(Bw))),1.0)  ! KINEROS2 Mod.
 280               if(WP .gt. 0.0) then
 281                R = AREA/ WP
 282               else
 283                R = 0.0
 284               endif
 285
 286             endif
 287
 288           if ( (h_0 .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) ) then !water outside of defined channel
 289             !weight the celerity by the contributing area, and assume that the mannings
 290             !of the spills is 2x the manning of the channel
 291                 Ck = max(0.0,((sqrt(So)/n)*((5./3.)*R**(2./3.) - &
 292                 ((2./3.)*R**(5./3.)*(2.0*sqrt(1.0 + z*z)/(Bw+2.0*bfd*z))))*AREA &
 293                 + ((sqrt(So)/(nCC))*(5./3.)*(h_0-bfd)**(2./3.))*AREAC)/(AREA+AREAC))
 294           else
 295                if(h_0 .gt. 0.0) then
 296                  Ck = max(0.0,(sqrt(So)/n)*((5./3.)*R**(2./3.) - &
 297                  ((2./3.)*R**(5./3.)*(2.0*sqrt(1.0 + z*z)/(Bw+2.0*h_0*z)))))
 298                 else
 299                  Ck = 0.0
 300                 endif
 301           endif
 302
 303           if(Ck .gt. 0.000000) then
 304             Km = max(dt,dx/Ck)
 305           else
 306             Km = dt
 307           endif
 308
 309           if ( (h_0 .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) .and. (Ck .gt. 0.0) ) then
 310           !water outside of defined channel
 311              X = min(0.5,max(0.0,0.5*(1-(Qj_0/(2.0*TwCC*So*Ck*dx)))))
 312           else
 313             if(Ck .gt. 0.0) then
 314               X = min(0.5,max(0.0,0.5*(1-(Qj_0/(2.0*Twl*So*Ck*dx)))))
 315             else
 316               X = 0.5
 317             endif
 318           endif
 319
 320            D = (Km*(1.000 - X) + dt/2.0000)              !--seconds
 321             if(D .eq. 0.0) then
 322               print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
 323               call hydro_stop("In MUSKINGCUNGE() - D is 0.")
 324            endif
 325
 326            C1 =  (Km*X + dt/2.000000)/D
 327            C2 =  (dt/2.0000 - Km*X)/D
 328            C3 =  (Km*(1.00000000-X)-dt/2.000000)/D
 329 !          C1 =  max(0.0,min(1.0,1.0-C3))
 330            if(ChannK .le. 0.0) then ! Disable channel loss for zero ChannK
 331              C4 =  (ql*dt)/D
 332            else
 333              !C4 =  (ql*dt)/D - (ChannK * dx * WPk)  !-- DY & LKR lat inflow + channel loss
 334              C4 =  ((ql - (ChannK * dx * WPk))*dt)/D  !-- TML: Loss adjusted by dt/D, closes water balance
 335            endif
 336
 337            if((WP+WPC) .gt. 0.0) then  !avoid divide by zero
 338 ! Another sanity check, prevent channel loss from exceeding flow
 339              if( (C4 .lt. 0.0) .and. (abs(C4) .gt. (C1*qup)+(C2*quc)+(C3*qdp)) )  then
 340               C4 = -((C1*qup)+(C2*quc)+(C3*qdp))
 341              endif
 342              Qj_0 =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) - ((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * &
 343                     (AREA+AREAC) * (R**(2./3.)) * sqrt(So)) !f0(x)
 344            endif
 345 !           WPk = WP*MIN((h/(modK*SQRT(Bw))),1.0)  ! KINEROS2 Mod. This shouldn't be HERE.
 346
 347            !--upper interval -----------
 348            WPC    = 0.0
 349            AREAC  = 0.0
 350
 351            Twl = Bw + 2.0*z*h                    !--top width of the channel inflow
 352
 353            if ( (h .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) ) then !water outside of defined channel
 354              AREA =  (Bw + bfd * z) * bfd
 355              AREAC = (TwCC * (h-bfd)) !assume compound component is rect. chan, that's 3 times the Tw
 356              WP = (Bw + 2.0 * bfd * sqrt(1.0 + z*z))
 357              WPC = TwCC + (2.0*(h-bfd)) !the additional wetted perimeter
 358              WPk = WP + WPC*MIN((h/(modK*SQRT(Bw))),1.0)  ! KINEROS2 Mod.
 359              R   = (AREA + AREAC)/(WP +WPC)
 360             ! RC  = AREAC/WPC
 361            else
 362               AREA = (Bw + h * z ) * h
 363               WP = (Bw + 2.0 * h * sqrt(1.000000 + z*z))
 364               WPk = WP*MIN((h/(modK*SQRT(Bw))),1.0)  ! KINEROS2 Mod.
 365               if(WP .gt. 0.0) then
 366                R = AREA/WP
 367               else
 368                R = 0.0
 369               endif
 370            endif
 371
 372           if ( (h .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) ) then
 373           !water outside of defined channel, assumed rectangular, 3x TW and n = 3x
 374                 Ck = max(0.0,((sqrt(So)/n)*((5./3.)*R**(2./3.) - &
 375                 ((2./3.)*R**(5./3.)*(2.0*sqrt(1.0 + z*z)/(Bw + 2.0*bfd*z))))*AREA &
 376                 + ((sqrt(So)/(nCC))*(5./3.)*(h-bfd)**(2./3.))*AREAC)/(AREA+AREAC))
 377           else
 378                if(h .gt. 0.0) then !avoid divide by zero
 379                  Ck = max(0.0,(sqrt(So)/n)*((5./3.)*R**(2./3.) - &
 380                  ((2./3.)*R**(5./3.)*(2.0 * sqrt(1.0 + z*z)/(Bw + 2.0*h*z)))))
 381                else
 382                  Ck = 0.0
 383                endif
 384           endif
 385
 386            if(Ck .gt. 0.0) then
 387             Km = max(dt,dx/Ck)
 388            else
 389             Km = dt
 390            endif
 391
 392           if ( (h .gt. bfd) .and. (TwCC .gt. 0.0) .and. (nCC .gt. 0.0) .and. (Ck .gt. 0.0) ) then
 393           !water outside of defined channel
 394             X = min(0.5,max(0.25,0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*TwCC*So*Ck*dx)))))
 395           else
 396             if(Ck .gt. 0.0) then
 397              X = min(0.5,max(0.25,0.5*(1-(((C1*qup)+(C2*quc)+(C3*qdp) + C4)/(2.0*Twl*So*Ck*dx)))))
 398             else
 399              X = 0.5
 400             endif
 401           endif
 402
 403            D = (Km*(1 - X) + dt/2)              !--seconds
 404             if(D .eq. 0.0) then
 405               print *, "FATAL ERROR: D is 0 in MUSKINGCUNGE", Km, X, dt,D
 406               call hydro_stop("In MUSKINGCUNGE() - D is 0.")
 407            endif
 408
 409            C1 =  (Km*X + dt/2.000000)/D
 410            C2 =  (dt/2.000000 - Km*X)/D
 411            C3 =  (Km*(1.000000-X)-dt/2.000000)/D
 412 !          C1 =  max(0.0,min(1.0,1.0-C3)) !! eliminate influence of upstream current
 413            if(ChannK .le. 0.0) then ! Disable channel loss for zero ChannK
 414              C4 =  (ql*dt)/D
 415            else
 416              !C4 =  (ql*dt)/D  - (ChannK * dx * WPk)  !-- (loss units: m/s * m * m)
 417              C4 =  ((ql - (ChannK * dx * WPk))*dt)/D  !-- TML: Loss adjusted by dt/D, closes water balance
 418            endif
 419
 420            if( (C4 .lt. 0.0) .and. (abs(C4) .gt. (C1*qup)+(C2*quc)+(C3*qdp)))  then
 421             C4 = -((C1*qup)+(C2*quc)+(C3*qdp))
 422            endif
 423
 424            if((WP+WPC) .gt. 0.0) then
 425             if( (C4 .lt. 0.0) .and. (abs(C4) .gt. (C1*qup)+(C2*quc)+(C3*qdp)))  then
 426                C4 = -((C1*qup)+(C2*quc)+(C3*qdp))
 427             endif
 428             Qj =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) - ((1.0000000/(((WP*n)+(WPC*nCC))/(WP+WPC))) * &
 429                     (AREA+AREAC) * (R**(2./3.)) * sqrt(So)) !f(x)
 430            endif
 431
 432            if(Qj_0-Qj .ne. 0.0) then
 433              h_1 = h - ((Qj * (h_0 - h))/(Qj_0 - Qj)) !update h, 3rd estimate
 434               if(h_1 .lt. 0.0) then
 435                 h_1 = h
 436               endif
 437            else
 438              h_1 = h
 439            endif
 440
 441            if(h .gt. 0.0) then
 442              rerror = abs((h_1 - h)/h) !relative error is new estatimate and 2nd estimate
 443              aerror = abs(h_1 -h)      !absolute error
 444            else
 445              rerror = 0.0
 446              aerror = 0.9
 447            endif
 448
 449 !          if(idx .eq. 6276407) then
 450 !            print*, "err,itr,hs", rerror, iter, depth, h_0, h, h_1
 451 !          endif
 452
 453            h_0  = max(0.0,h)
 454            h    = max(0.0,h_1)
 455            iter = iter + 1
 456
 457            if( h .lt. mindepth) then  ! exit loop if depth is very small
 458              goto 111
 459            endif
 460
 461          end do
 462
 463 !          if(idx .eq.  6276407) then
 464 !            print*, "id,itr,err,h", idx, iter, rerror, h
 465 !          endif
 466
 467 111      continue
 468
 469          if(iter .ge. maxiter) then
 470
 471            tries = tries + 1
 472            if(tries .le. 4) then  ! expand the search space
 473              h     =  h * 1.33
 474              h_0   =  h_0 * 0.67
 475              maxiter = maxiter + 25 !and increase the number of allowable iterations
 476             goto 110
 477            endif
 478
 479            print*, "Musk Cunge WARNING: Failure to converge"
 480            print*, 'RouteLink index:', idx + linkls_s(my_id+1) - 1
 481            print*, "id,err,iters,tries", idx, rerror, iter, tries
 482            print*, "Ck,X,dt,Km",Ck,X,dt,Km
 483            print*, "So,dx,h",So,dx,h
 484            print*, "qup,quc,qdp,ql", qup,quc,qdp,ql
 485            print*, "bfd,Bw,Tw,Twl", bfd,Bw,Tw,Twl
 486            print*, "Qmc,Qmn", (C1*qup)+(C2*quc)+(C3*qdp) + C4,((1/(((WP*n)+(WPC*nCC))/(WP+WPC))) * &
 487                     (AREA+AREAC) * (R**(2./3.)) * sqrt(So))
 488          endif
 489
 490
 491 #ifdef HYDRO_D
 492       !! Getting information on c3.
 493       !if(rt_domain(1)%gages(idx + linkls_s(my_id+1) - 1) .ne. rt_domain(1)%gageMiss) then
 494       !   print*, rt_domain(1)%gages(idx+linkls_s(my_id+1)-1)
 495 !     if(rt_domain(1)%gages(idx) .ne. rt_domain(1)%gageMiss) then
 496 !        print*, rt_domain(1)%gages(idx)
 497 !        print*,'JLM submuskcunge idx, C3, C, D:', idx + linkls_s(my_id+1) - 1, C3, dt/(2*Km), X
 498 !     end if
 499 #endif
 500
 501 !yw added for test
 502       !DY and LKR Added to update for channel loss
 503         if(((C1*qup)+(C2*quc)+(C3*qdp) + C4) .lt. 0.0) then
 504 !       MUSKINGCUNGE =  MAX( ( (C1*qup)+(C2*quc) + C4),((C1*qup)+(C3*qdp) + C4) )
 505            if( (C4 .lt. 0.0) .and. (abs(C4) .gt. (C1*qup)+(C2*quc)+(C3*qdp)) )  then ! channel loss greater than water in chan
 506              qdc = 0.0
 507            else
 508              qdc = MAX( ( (C1*qup)+(C2*quc) + C4),((C1*qup)+(C3*qdp) + C4) )
 509            endif
 510         else
 511 !       MUSKINGCUNGE =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber
 512           qdc =  ((C1*qup)+(C2*quc)+(C3*qdp) + C4) !-- pg 295 Bedient huber
 513
 514         endif
 515
 516         Twl = Bw + (2.0*z*h)
 517         R = (h*(Bw + Twl) / 2.0) / (Bw + 2.0*(((Twl - Bw) / 2.0)**2.0 + h**2)**0.5)
 518         vel =  (1./n) * (R **(2.0/3.0)) * sqrt(So)  ! average velocity in m/s
 519         depth = h
 520
 521       else   ! no flow to route
 522        qdc = 0.0
 523        depth = 0.0
 524      endif
 525
 526       qloss = (ChannK * dx * WPk) ! TML -- Compute qloss to pass
 527
 528       !TML:Added print statement to test qlos function;
 529       !comment out to prevent excessive file sizes when running model
 530       !print*, "qloss,dx,WP,WPk,depth,ChannK,qdc,ql,dt,D", qloss,dx,WP,WPk,depth,ChannK,qdc,ql,dt,D
 531       if((qloss*dt)/D > ((ql*dt)/D - C4)) then
 532            qloss = ql - C4*(D/dt)
 533            if (qloss < 0) then
 534               print*, 'WARNING CHANNEL LOSS IS NEGATIVE',qloss
 535            endif
 536       endif
 537
 538 ! ----------------------------------------------------------------
 539 END SUBROUTINE SUBMUSKINGCUNGE
 540 ! ----------------------------------------------------------------
 541
 542 ! ------------------------------------------------
 543 !   FUNCTION KINEMATIC
 544 ! ------------------------------------------------
 545         REAL FUNCTION KINEMATIC()
 546
 547         IMPLICIT NONE
 548
 549 ! -------- DECLARATIONS -----------------------
 550
 551 !       REAL, INTENT(OUT), DIMENSION(IXRT,JXRT) :: OVRGH
 552
 553         KINEMATIC = 1
 554 !----------------------------------------------------------------
 555   END FUNCTION KINEMATIC
 556 !----------------------------------------------------------------
 557
 558
 559 ! ------------------------------------------------
 560 !   SUBROUTINE drive_CHANNEL
 561 ! ------------------------------------------------
 562 ! ------------------------------------------------
 563      Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, &
 564        QSUBRT, LAKEINFLORT, QSTRMVOLRT, TO_NODE, FROM_NODE, &
 565        TYPEL, ORDER, MAXORDER, NLINKS, CH_NETLNK, CH_NETRT, CH_LNKRT, &
 566        LAKE_MSKRT, DT, DTCT, DTRT_CH,MUSK, MUSX, QLINK, &
 567        QLateral, &
 568        HLINK, ELRT, CHANLEN, MannN, So, ChSSlp, Bw, Tw, Tw_CC, n_CC,  &
 569        ChannK, RESHT, &
 570        ZELEV, CVOL, NLAKES, QLAKEI, QLAKEO, LAKENODE, &
 571        dist, QINFLOWBASE, CHANXI, CHANYJ, channel_option, RETDEP_CHAN, &
 572        NLINKSL, LINKID, node_area  &
 573 #ifdef MPP_LAND
 574        , lake_index,link_location,mpp_nlinks,nlinks_index,yw_mpp_nlinks  &
 575        , LNLINKSL  &
 576        , gtoNode,toNodeInd,nToNodeInd &
 577 #endif
 578        , CH_LNKRT_SL &
 579        ,gwBaseSwCRT, gwHead, qgw_chanrt, gwChanCondSw, gwChanCondConstIn, &
 580        gwChanCondConstOut, velocity, qloss)
 581
 582
 583        IMPLICIT NONE
 584
 585         ! -------- DECLARATIONS ------------------------
 586         INTEGER, INTENT(IN) :: did, IXRT,JXRT,channel_option
 587         INTEGER, INTENT(IN) :: NLINKS,NLAKES, NLINKSL
 588         integer, INTENT(INOUT) :: KT   ! flag of cold start (1) or continue run.
 589         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QSUBRT
 590         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QSTRMVOLRT
 591         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: LAKEINFLORT
 592         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: ELRT
 593         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: QINFLOWBASE
 594         INTEGER(kind=int64), INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETLNK
 595
 596         INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
 597         INTEGER(kind=int64), INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT
 598         INTEGER(kind=int64), INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_LNKRT_SL
 599
 600        real , dimension(ixrt,jxrt):: latval,lonval
 601
 602         INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
 603         INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: ORDER, TYPEL !--link
 604         INTEGER(kind=int64), INTENT(IN), DIMENSION(NLINKS)    :: TO_NODE, FROM_NODE
 605         INTEGER, INTENT(IN), DIMENSION(NLINKS)    :: CHANXI, CHANYJ
 606         REAL,    INTENT(IN), DIMENSION(NLINKS)    :: ZELEV  !--elevation of nodes
 607         REAL, INTENT(INOUT), DIMENSION(NLINKS)    :: CVOL
 608         REAL, INTENT(IN), DIMENSION(NLINKS)       :: MUSK, MUSX
 609         REAL, INTENT(IN), DIMENSION(NLINKS)       :: CHANLEN
 610         REAL, INTENT(IN), DIMENSION(NLINKS)       :: So, MannN
 611         REAL, INTENT(IN), DIMENSION(NLINKS)       :: ChSSlp,Bw,Tw  !--properties of nodes or links
 612         REAL, INTENT(IN), DIMENSION(NLINKS)       :: Tw_CC, n_CC  !properties of compund channel
 613         REAL, INTENT(IN), DIMENSION(NLINKS)       :: ChannK  !--Channel Infiltration
 614         REAL                                      :: Km, X
 615         REAL , INTENT(INOUT), DIMENSION(:,:) :: QLINK
 616         REAL ,  DIMENSION(NLINKS,2) :: tmpQLINK
 617         REAL , INTENT(INOUT), DIMENSION(NLINKS)   :: HLINK
 618         REAL, dimension(NLINKS), intent(inout)    :: QLateral !--lateral flow
 619         REAL, INTENT(IN)                          :: DT    !-- model timestep
 620         REAL, INTENT(IN)                          :: DTRT_CH  !-- routing timestep
 621         REAL, INTENT(INOUT)                       :: DTCT
 622         real                                      :: minDTCT !BF minimum routing timestep
 623         REAL                                      :: dist(ixrt,jxrt,9)
 624         REAL                                      :: RETDEP_CHAN
 625         INTEGER, INTENT(IN)                       :: MAXORDER, SUBRTSWCRT, &
 626                                                      gwBaseSwCRT, gwChanCondSw
 627         real, intent(in)                          :: gwChanCondConstIn, gwChanCondConstOut ! aquifer-channel conductivity constant from namelist
 628         REAL , INTENT(IN), DIMENSION(NLINKS)      :: node_area
 629         real, dimension(:), INTENT(inout)         :: velocity, qloss
 630
 631
 632 !DJG GW-chan coupling variables...
 633         REAL, DIMENSION(NLINKS)                   :: dzGwChanHead
 634         REAL, DIMENSION(NLINKS)                   :: Q_GW_CHAN_FLUX     !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state...
 635         REAL, DIMENSION(IXRT,JXRT)                :: ZWATTBLRT          !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state...
 636         REAL, INTENT(INOUT), DIMENSION(:,:), allocatable :: gwHead      !DJG !!! groundwater head from Fersch-2d gw implementation...units (m ASL)
 637         REAL, INTENT(INOUT), DIMENSION(:,:), allocatable :: 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...
 638
 639
 640
 641         !-- lake params
 642         REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: RESHT    !-- reservoir height (m)
 643         REAL*8,  DIMENSION(NLAKES)    :: QLAKEI8   !-- lake inflow (cms)
 644         REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: QLAKEI   !-- lake inflow (cms)
 645         REAL,                DIMENSION(NLAKES)    :: QLAKEIP  !-- lake inflow previous timestep (cms)
 646         REAL, INTENT(INOUT), DIMENSION(NLAKES)    :: QLAKEO   !-- outflow from lake used in diffusion scheme
 647         INTEGER(kind=int64), INTENT(IN), DIMENSION(NLINKS)    :: LAKENODE !-- outflow from lake used in diffusion scheme
 648         INTEGER(kind=int64), INTENT(IN), DIMENSION(NLINKS)    :: LINKID   !--  id of channel elements for linked scheme
 649         !REAL, DIMENSION(NLINKS)                   :: QLateral !--lateral flow
 650         REAL, DIMENSION(NLINKS)                   :: QSUM     !--mass bal of node
 651         REAL*8, DIMENSION(NLINKS)                   :: QSUM8     !--mass bal of node
 652         REAL, DIMENSION(NLAKES)                   :: QLLAKE   !-- lateral inflow to lake in diffusion scheme
 653         REAL*8, DIMENSION(NLAKES)                   :: QLLAKE8   !-- lateral inflow to lake in diffusion scheme
 654
 655 !-- Local Variables
 656         INTEGER                     :: i,j,k,t,m,jj,kk,KRT,node
 657         INTEGER                     :: DT_STEPS               !-- number of timestep in routing
 658         REAL                        :: Qup,Quc                !--Q upstream Previous, Q Upstream Current, downstream Previous
 659         REAL                        :: bo                     !--critical depth, bnd outflow just for testing
 660         REAL                        :: AREA,WP                !--wetted area and perimiter for MuskingC. routing
 661
 662         REAL ,DIMENSION(NLINKS)     :: HLINKTMP,CVOLTMP       !-- temporarily store head values and volume values
 663         REAL ,DIMENSION(NLINKS)     :: CD                     !-- critical depth
 664         real, DIMENSION(IXRT,JXRT)  :: tmp
 665         real, dimension(nlinks)     :: tmp2
 666
 667 #ifdef MPP_LAND
 668         integer lake_index(nlakes)
 669         integer nlinks_index(nlinks)
 670         integer mpp_nlinks, iyw, yw_mpp_nlinks
 671         integer(kind=int64) link_location(ixrt,jxrt)
 672         real     ywtmp(ixrt,jxrt)
 673         integer LNLINKSL
 674         real*8,  dimension(LNLINKSL) :: LQLateral
 675 !        real*4,  dimension(LNLINKSL) :: LQLateral
 676         integer, dimension(:) ::  toNodeInd
 677         integer(kind=int64), dimension(:,:) ::  gtoNode
 678         integer  :: nToNodeInd
 679         real, dimension(nToNodeInd,2) :: gQLINK
 680 #else
 681         real*8, dimension(NLINKS)                   :: LQLateral !--lateral flow
 682 #endif
 683         integer flag
 684
 685         integer :: n, kk2, nt, nsteps  ! tmp
 686
 687         QLAKEIP = 0
 688         QLAKEI8 = 0
 689         HLINKTMP = 0
 690         CVOLTMP = 0
 691         CD = 0
 692         node = 1
 693         QLateral = 0
 694         QSUM     = 0
 695         QLLAKE   = 0
 696
 697
 698 !yw      print *, "DRIVE_channel,option,nlinkl,nlinks!!", channel_option,NLINKSL,NLINKS
 699 !      print *, "DRIVE_channel, RESHT", RESHT
 700
 701
 702       dzGwChanHead = 0.
 703
 704    IF(channel_option .ne. 3) then   !--muskingum methods ROUTE ON DT timestep, not DTRT!!
 705
 706          nsteps = (DT+0.5)/DTRT_CH
 707
 708 #ifdef MPP_LAND
 709          LQLateral = 0          !-- initial lateral flow to 0 for this reach
 710          DO iyw = 1,yw_MPP_NLINKS
 711          jj = nlinks_index(iyw)
 712           !--------river grid points, convert depth in mm to rate across reach in m^3/sec
 713               if( .not. (  (CHANXI(jj) .eq. 1 .and. left_id .ge. 0) .or. &
 714                            (CHANXI(jj) .eq. ixrt .and. right_id .ge. 0) .or. &
 715                            (CHANYJ(jj) .eq. 1 .and. down_id .ge. 0) .or. &
 716                            (CHANYJ(jj) .eq. jxrt .and. up_id .ge. 0)      &
 717                    ) ) then
 718                   if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0) then
 719                      k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
 720                      LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
 721                             *node_area(jj)/DT)
 722                   elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
 723                       k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
 724                       LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
 725                                *node_area(jj)/DT)
 726                   endif
 727               endif
 728          end do  ! jj
 729
 730
 731 !   assign LQLATERAL to QLATERAL
 732        call updateLinkV(LQLateral, QLateral(1:NLINKSL))
 733
 734 #else
 735
 736          LQLateral = 0          !-- initial lateral flow to 0 for this reach
 737          do jj = 1, NLINKS
 738           !--------river grid points, convert depth in mm to rate across reach in m^3/sec
 739
 740                   if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0 ) then
 741                      k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
 742                      LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
 743                             *node_area(jj)/DT)
 744                   elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
 745                       k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
 746                       LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
 747                                *node_area(jj)/DT)
 748                   endif
 749
 750           end do  ! jj
 751           QLateral = LQLateral
 752
 753 #endif
 754
 755 !       QLateral = QLateral / nsteps
 756
 757    do nt = 1, nsteps
 758
 759 !Per Yates, this check is not needed. Commenting out for now.
 760 !----------  route order 1 reaches which have no upstream inflow
 761 !        do k=1, NLINKSL
 762 !           if (ORDER(k) .eq. 1) then  !-- first order stream has no headflow
 763
 764
 765 !              if(TYPEL(k) .eq. 1) then    !-- level pool route of reservoir
 766 !                  !CALL LEVELPOOL(1,0.0, 0.0, qd, QLINK(k,2), QLateral(k), &
 767 !                  ! DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), &
 768 !                  ! WEIRC(k), WEIRL(k), ORIFICEE(i), ORIFICEC(k), ORIFICEA(k) )
 769 !              elseif (channel_option .eq. 1) then
 770 !                   Km  = MUSK(k)
 771 !                   X   = MUSX(k)
 772 !                   QLINK(k,2) = MUSKING(k,0.0, QLateral(k), QLINK(k,1), DTRT_CH, Km, X) !--current outflow
 773 !              elseif (channel_option .eq. 2) then !-- upstream is assumed constant initial condition
 774
 775 !                   call SUBMUSKINGCUNGE(QLINK(k,2), velocity(k), k,  &
 776 !                    0.0,0.0, QLINK(k,1), QLateral(k),   DTRT_CH, So(k), &
 777 !                    CHANLEN(k), MannN(k), ChSSlp(k), Bw(k), Tw(k) )
 778
 779 !              else
 780 !                  print *, "FATAL ERROR: No channel option selected"
 781 !                  call hydro_stop("In drive_CHANNEL() -No channel option selected ")
 782 !              endif
 783 !           endif
 784 !        end do
 785
 786 #ifdef MPP_LAND
 787        gQLINK = 0.000
 788        call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
 789        call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1))
 790 #endif
 791
 792       !---------- route other reaches, with upstream inflow
 793        tmpQlink = 0.0
 794        do k = 1,NLINKSL
 795 !          if (ORDER(k) .gt. 1 ) then  !-- exclude first order stream
 796              Quc  = 0.0
 797              Qup  = 0.0
 798
 799 #ifdef MPP_LAND
 800 !using mapping index
 801                do n = 1, gtoNODE(k,1)
 802                   m = gtoNODE(k,n+1)
 803 !yw                  if (LINKID(k) .eq. m) then
 804                     Quc = Quc + gQLINK(m,2)  !--accum of upstream inflow of current timestep (2)
 805                     Qup = Qup + gQLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
 806
 807                       !     if(LINKID(k) .eq. 3259 .or. LINKID(k) .eq. 3316 .or. LINKID(k) .eq. 3219) then
 808                       !       write(6,*) "id,Uc,Up",LINKID(k),Quc,Qup
 809                       !       call flush(6)
 810                       !     endif
 811
 812 !yw                  endif
 813                 end do ! do i
 814
 815 #else
 816                do m = 1, NLINKSL
 817                   if (LINKID(k) .eq. TO_NODE(m)) then
 818                     Quc = Quc + QLINK(m,2)  !--accum of upstream inflow of current timestep (2)
 819                     Qup = Qup + QLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
 820                   endif
 821                end do ! do m
 822 #endif
 823
 824                 if(TYPEL(k) .eq. 1) then   !--link is a reservoir
 825
 826                    ! CALL LEVELPOOL(1,QLINK(k,1), Qup, QLINK(k,1), QLINK(k,2), &
 827                    !  QLateral(k), DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), &
 828                    !  WEIRC(k), WEIRL(k),ORIFICEE(k),  ORIFICEC(k), ORIFICEA(k))
 829
 830                    elseif (channel_option .eq. 1) then  !muskingum routing
 831                        Km = MUSK(k)
 832                        X = MUSX(k)
 833                        tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow
 834                    elseif (channel_option .eq. 2) then ! muskingum cunge
 835
 836                    call SUBMUSKINGCUNGE(tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k),  &
 837                     Qup,Quc, QLINK(k,1), QLateral(k),   DTRT_CH, So(k), &
 838                     CHANLEN(k), MannN(k), ChSSlp(k), Bw(k), Tw(k),Tw_CC(k), n_CC(k),  HLINK(k), ChannK(k) )
 839
 840                    else
 841                     print *, "FATAL ERROR: no channel option selected"
 842                     call hydro_stop("In drive_CHANNEL() - no channel option selected")
 843                    endif
 844 !            endif !!! order(1) .ne. 1
 845          end do       !--k links
 846
 847 !yw check
 848 !        gQLINK = 0.0
 849 !        call ReachLS_write_io(tmpQLINK(:,2), gQLINK(:,2))
 850 !        call ReachLS_write_io(tmpQLINK(:,1), gQLINK(:,1))
 851 !        write(6,*) " io_id = ", io_id
 852 !        if(my_id .eq. io_id) then
 853 !            write(71,*) gQLINK(:,1)
 854 !            call flush(71)
 855 !            call flush(72)
 856 !        endif
 857
 858           do k = 1, NLINKSL
 859             if(TYPEL(k) .ne. 1) then
 860                QLINK(k,2) = tmpQLINK(k,2)
 861             endif
 862             QLINK(k,1) = QLINK(k,2)    !assing link flow of current to be previous for next time step
 863          end do
 864
 865    end do ! nsteps
 866
 867 #ifdef HYDRO_D
 868           print *, "END OF ALL REACHES...",KRT,DT_STEPS
 869 #endif
 870
 871 !    END DO !-- krt timestep for muksingumcunge routing
 872
 873    elseif(channel_option .eq. 3) then   !--- route using the diffusion scheme on nodes not links
 874
 875 #ifdef MPP_LAND
 876          call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINK,NLINKS,99)
 877          call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOL,NLINKS,99)
 878 #endif
 879
 880          KRT = 0                  !-- initialize the time counter
 881          minDTCT = 0.01           ! define minimum routing sub-timestep (s), simulation will end with smaller timestep
 882          DTCT = min(max(DTCT*2.0, minDTCT),DTRT_CH)
 883
 884          HLINKTMP = HLINK         !-- temporary storage of the water elevations (m)
 885          CVOLTMP = CVOL           !-- temporary storage of the volume of water in channel (m^3)
 886          QLAKEIP = QLAKEI         !-- temporary lake inflow from previous timestep  (cms)
 887
 888 !        call check_channel(77,HLINKTMP,1,nlinks)
 889 !        call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,ZELEV,NLINKS,99)
 890  crnt:   DO                      !-- loop on the courant condition
 891           QSUM     = 0.              !-- initialize the total flow out of each cell to zero
 892           QSUM8     = 0.              !-- initialize the total flow out of each cell to zero
 893           QLAKEI8  = 0.              !-- set the lake inflow as zero
 894           QLAKEI   = 0.              !-- set the lake inflow as zero
 895           QLLAKE   = 0.              !-- initialize each lake's lateral inflow to zero
 896           QLLAKE8   = 0.              !-- initialize each lake's lateral inflow to zero
 897           DT_STEPS = INT(DT/DTCT)   !-- fix the timestep
 898           QLateral = 0.
 899 !DJG GW-chan coupling variables...
 900           if(gwBaseSwCRT == 3) then
 901           Q_GW_CHAN_FLUX = 0.
 902           qgw_chanrt     = 0.
 903           end if
 904
 905 !         ZWATTBLRT=1.0   !--HARDWIRE, remove this and pass in from subsfc/gw routing routines...
 906
 907
 908 !-- vectorize
 909 !---------------------
 910 #ifdef MPP_LAND
 911          DO iyw = 1,yw_MPP_NLINKS
 912          i = nlinks_index(iyw)
 913 #else
 914          DO i = 1,NLINKS
 915 #endif
 916
 917            if(node_area(i) .eq. 0) then
 918                write(6,*) "FATAL ERROR: node_area(i) is zero. i=", i
 919                call hydro_stop("In drive_CHANNEL() - Error node_area")
 920            endif
 921
 922
 923
 924 nodeType:if((CH_NETRT(CHANXI(i), CHANYJ(i) ) .eq. 0) .and. &
 925               (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .lt.0) ) then !--a reg. node
 926
 927 gwOption:   if(gwBaseSwCRT == 3) then
 928
 929              ! determine potential gradient between groundwater head and channel stage
 930              ! units in (m)
 931              dzGwChanHead(i) = gwHead(CHANXI(i),CHANYJ(i)) - (HLINK(i)+ZELEV(i))
 932
 933              if(gwChanCondSw .eq. 0) then
 934
 935                 qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
 936
 937              else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) > 0) then
 938
 939                ! channel bed interface, units in (m^3/s), flux into channel...
 940                ! BF todo: consider channel width
 941                 qgw_chanrt(CHANXI(i),CHANYJ(i)) = gwChanCondConstIn * dzGwChanHead(i) &
 942                                                 * CHANLEN(i) * 2.
 943
 944              else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) < 0) then
 945
 946                ! channel bed interface, units in (m^3/s), flux out of channel...
 947                ! BF todo: consider channel width
 948                 qgw_chanrt(CHANXI(i),CHANYJ(i)) = max(-0.005, gwChanCondConstOut * dzGwChanHead(i) &
 949                                                 * CHANLEN(i) * 2.)
 950 !              else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then  TBD: exponential dependency
 951 !              else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then
 952
 953              else
 954
 955                 qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
 956
 957              end if
 958
 959              Q_GW_CHAN_FLUX(i) = qgw_chanrt(CHANXI(i),CHANYJ(i))
 960 !             if ( i .eq. 1001 ) then
 961 !                print *, Q_GW_CHAN_FLUX(i), dzGwChanHead(i), ELRT(CHANXI(i),CHANYJ(i)), HLINK(i), ZELEV(i)
 962 !             end if
 963 !              if ( Q_GW_CHAN_FLUX(i) .lt. 0. ) then   !-- temporary hardwire for only allowing flux into channel...REMOVE later...
 964 !                 Q_GW_CHAN_FLUX(i) = 0.
 965 !               qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
 966 !              end if
 967
 968             else
 969               Q_GW_CHAN_FLUX(i) = 0.
 970             end if gwOption
 971
 972
 973               QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) =  &
 974 !DJG  awaiting gw-channel exchg...  Q_GW_CHAN_FLUX(i)+& ...obsolete-> ((QSUBRT(CHANXI(i),CHANYJ(i))+&
 975                 Q_GW_CHAN_FLUX(i)+&
 976                 ((QSTRMVOLRT(CHANXI(i),CHANYJ(i))+&
 977                  QINFLOWBASE(CHANXI(i),CHANYJ(i))) &
 978                    /DT_STEPS*node_area(i)/1000/DTCT)
 979                if((QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))).lt.0.) .and. (gwChanCondSw == 0)) then
 980 #ifdef HYDRO_D
 981                print*, "i, CHANXI(i),CHANYJ(i) = ", i, CHANXI(i),CHANYJ(i)
 982                print *, "NEGATIVE Lat inflow...",QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))), &
 983                          QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)), &
 984                          QINFLOWBASE(CHANXI(i),CHANYJ(i))
 985 #endif
 986                elseif (QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) .gt. 1.0) then
 987 !#ifdef HYDRO_D
 988 !               print *, "LatIn(Ql,Qsub,Qstrmvol)..",i,QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))), &
 989 !                          QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i))
 990 !#endif
 991                end if
 992
 993          elseif(LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .gt. 0 .and. &
 994 !               (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .ne. -9999)) then !--a lake node
 995                 (CH_NETRT(CHANXI(i),CHANYJ(i)) .le. 0)) then !--a lake node
 996               QLLAKE8(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) = &
 997                  QLLAKE8(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) + &
 998                  (LAKEINFLORT(CHANXI(i),CHANYJ(i))+ &
 999                  QINFLOWBASE(CHANXI(i),CHANYJ(i))) &
1000                  /DT_STEPS*node_area(i)/1000/DTCT
1001          elseif(CH_NETRT(CHANXI(i),CHANYJ(i)) .gt. 0) then  !pour out of lake
1002                  QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) =  &
1003                    QLAKEO(CH_NETRT(CHANXI(i),CHANYJ(i)))  !-- previous timestep
1004          endif nodeType
1005
1006         ENDDO
1007
1008
1009 #ifdef MPP_LAND
1010     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLateral,NLINKS,99)
1011     if(NLAKES .gt. 0) then
1012        !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT   ,ixrt,jxrt,QLLAKE,NLAKES,99)
1013        call sum_real8(QLLAKE8,NLAKES)
1014        QLLAKE = QLLAKE8
1015     endif
1016 #endif
1017
1018           !-- compute conveyances, with known depths (just assign to QLINK(,1)
1019           !--QLINK(,2) will not be used), QLINK is the flow across the node face
1020           !-- units should be m3/second.. consistent with QL (lateral flow)
1021
1022 #ifdef MPP_LAND
1023          DO iyw = 1,yw_MPP_NLINKS
1024          i = nlinks_index(iyw)
1025 #else
1026            DO i = 1,NLINKS
1027 #endif
1028            if (TYPEL(i) .eq. 0 .AND. HLINKTMP(FROM_NODE(i)) .gt. RETDEP_CHAN) then
1029                if(from_node(i) .ne. to_node(i) .and. (to_node(i) .gt. 0) .and.(from_node(i) .gt. 0) ) &  ! added by Wei Yu
1030                    QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
1031                      HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
1032                      CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
1033             else !--  we are just computing critical depth for outflow points
1034                QLINK(i,1) =0.
1035             endif
1036           ENDDO
1037
1038 #ifdef MPP_LAND
1039     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
1040 #endif
1041
1042
1043           !-- compute total flow across face, into node
1044 #ifdef MPP_LAND
1045          DO iyw = 1,yw_mpp_nlinks
1046          i = nlinks_index(iyw)
1047 #else
1048           DO i = 1,NLINKS                                                 !-- inflow to node across each face
1049 #endif
1050            if(TYPEL(i) .eq. 0) then                                       !-- only regular nodes have to attribute
1051               QSUM8(TO_NODE(i)) = QSUM8(TO_NODE(i)) + QLINK(i,1)
1052            endif
1053           END DO
1054
1055 #ifdef MPP_LAND
1056     call MPP_CHANNEL_COM_REAL8(Link_location,ixrt,jxrt,qsum8,NLINKS,0)
1057 #endif
1058     qsum = qsum8
1059
1060
1061
1062 #ifdef MPP_LAND
1063          do iyw = 1,yw_mpp_nlinks
1064          i = nlinks_index(iyw)
1065 #else
1066          do i = 1,NLINKS                                                 !-- outflow from node across each face
1067 #endif
1068             QSUM(FROM_NODE(i)) = QSUM(FROM_NODE(i)) - QLINK(i,1)
1069          end do
1070 #ifdef MPP_LAND
1071     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,qsum,NLINKS,99)
1072 #endif
1073
1074
1075          flag = 99
1076
1077
1078 #ifdef MPP_LAND
1079          do iyw = 1,yw_MPP_NLINKS
1080              i = nlinks_index(iyw)
1081 #else
1082          do i = 1, NLINKS                                                !--- compute volume and depth at each node
1083 #endif
1084
1085            if( TYPEL(i).eq.0 .and. CVOLTMP(i) .ge. 0.001 .and.(CVOLTMP(i)-QSUM(i)*DTCT)/CVOLTMP(i) .le. -0.01 ) then
1086             flag = -99
1087 #ifdef HYDRO_D
1088             write(6,*) "******* start diag ***************"
1089             write(6,*) "Unstable at node ",i, "i=",CHANXI(i),"j=",CHANYJ(i)
1090             write(6,*) "Unstatble at node ",i, "lat=",latval(CHANXI(i),CHANYJ(i)), "lon=",lonval(CHANXI(i),CHANYJ(i))
1091             write(6,*) "TYPEL, CVOLTMP, QSUM, QSUM*DTCT",TYPEL(i), CVOLTMP(i), QSUM(i), QSUM(i)*DTCT
1092             write(6,*) "qsubrt, qstrmvolrt,qlink",QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)),qlink(i,1),qlink(i,2)
1093 !              write(6,*) "current nodes, z, h", ZELEV(FROM_NODE(i)),HLINKTMP(FROM_NODE(i))
1094 !           if(TO_NODE(i) .gt. 0) then
1095 !              write(6,*) "to nodes, z, h", ZELEV(TO_NODE(i)), HLINKTMP(TO_NODE(i))
1096 !           else
1097 !              write(6,*) "no to nodes   "
1098 !           endif
1099                write(6,*) "CHANLEN(i), MannN(i), Bw(i), ChSSlp(i) ", CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)
1100             write(6,*) "*******end of  diag ***************"
1101 #endif
1102
1103             goto 999
1104             endif
1105           enddo
1106
1107 999 continue
1108 #ifdef MPP_LAND
1109         call mpp_same_int1(flag)
1110 #endif
1111
1112
1113         if(flag < 0  .and. DTCT >0.1)   then
1114
1115              ! call smoth121(HLINK,nlinks,maxv_p,pnode,to_node)
1116
1117              if(DTCT .gt. minDTCT) then                !-- timestep in seconds
1118               DTCT = max(DTCT/2 , minDTCT)             !-- 1/2 timestep
1119               KRT = 0                                  !-- restart counter
1120               HLINKTMP = HLINK                         !-- set head and vol to start value of timestep
1121               CVOLTMP = CVOL
1122               CYCLE crnt                               !-- start cycle over with smaller timestep
1123              else
1124               write(6,*) "Courant error with smallest routing timestep DTCT: ",DTCT
1125 !              call hydro_stop("drive_CHANNEL")
1126               DTCT = 0.1
1127               HLINKTMP = HLINK                          !-- set head and volume to start values of timestep
1128               CVOLTMP  = CVOL
1129               goto 998
1130              end if
1131         endif
1132
1133 998 continue
1134
1135
1136 #ifdef MPP_LAND
1137         do iyw = 1,yw_MPP_NLINKS
1138             i = nlinks_index(iyw)
1139 #else
1140         do i = 1, NLINKS                                                !--- compute volume and depth at each node
1141 #endif
1142
1143            if(TYPEL(i) .eq. 0) then                   !--  regular channel grid point, compute volume
1144               CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) )* DTCT
1145               if((CVOLTMP(i) .lt. 0) .and. (gwChanCondSw == 0)) then
1146 #ifdef HYDRO_D
1147                 print *, "WARNING! channel volume less than 0:i,CVOL,QSUM,QLat", &
1148                                i, CVOLTMP(i),QSUM(i),QLateral(i),HLINK(i)
1149 #endif
1150                 CVOLTMP(i) =0
1151               endif
1152
1153            elseif(TYPEL(i) .eq. 1) then               !-- pour point, critical depth downstream
1154
1155               if (QSUM(i)+QLateral(i) .lt. 0) then
1156               else
1157
1158 !DJG remove to have const. flux b.c....   CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
1159                   CD(i) = HLINKTMP(i)  !This is a temp hardwire for flow depth for the pour point...
1160               endif
1161
1162                ! change in volume is inflow, lateral flow, and outflow
1163                !yw DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*DXRT),HLINKTMP(i), &
1164                    CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
1165                        DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
1166                        CD(i),CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
1167           elseif (TYPEL(i) .eq. 2) then              !--- into a reservoir, assume critical depth
1168               if ((QSUM(i)+QLateral(i) .lt. 0) .and. (gwChanCondSw == 0)) then
1169 #ifdef HYDRO_D
1170                print *, i, 'CrtDpth Qsum+QLat into lake< 0',QSUM(i), QLateral(i)
1171 #endif
1172              else
1173 !DJG remove to have const. flux b.c....    CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
1174                CD(i) = HLINKTMP(i)  !This is a temp hardwire for flow depth for the pour point...
1175              endif
1176
1177               !-- compute volume in reach (m^3)
1178                    CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
1179                           DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
1180                              CD(i) ,CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
1181           else
1182               print *, "FATAL ERROR: This node does not have a type.. error TYPEL =", TYPEL(i)
1183               call hydro_stop("In drive_CHANNEL() - error TYPEL")
1184           endif
1185
1186           if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
1187               HLINKTMP(i) = HEAD(i, CVOLTMP(i)/CHANLEN(i),Bw(i),1/ChSSlp(i))  !--updated depth
1188           else
1189               HLINKTMP(i) = CD(i)  !!!   CRITICALDEPTH(i,QSUM(i)+QLateral(i), Bw(i), 1./ChSSlp(i)) !--critical depth is head
1190           endif
1191
1192           if(TO_NODE(i) .gt. 0) then
1193              if(LAKENODE(TO_NODE(i)) .gt. 0) then
1194                   QLAKEI8(LAKENODE(TO_NODE(i))) = QLAKEI8(LAKENODE(TO_NODE(i))) + QLINK(i,1)
1195              endif
1196           endif
1197
1198         END DO  !--- done processing all the links
1199
1200 #ifdef MPP_LAND
1201     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOLTMP,NLINKS,99)
1202     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CD,NLINKS,99)
1203     if(NLAKES .gt. 0) then
1204 !       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
1205         call sum_real8(QLAKEI8,NLAKES)
1206         QLAKEI = QLAKEI8
1207     endif
1208     call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINKTMP,NLINKS,99)
1209 #endif
1210 !   call check_channel(83,CVOLTMP,1,nlinks)
1211 !   call check_channel(84,CD,1,nlinks)
1212 !   call check_channel(85,HLINKTMP,1,nlinks)
1213 !   call check_lake(86,QLAKEI,lake_index,nlakes)
1214
1215            do i = 1, NLAKES !-- mass balances of lakes
1216 #ifdef MPP_LAND
1217             if(lake_index(i) .gt. 0)  then
1218 #endif
1219
1220               ! Calls run for a single reservoir depending on its type as
1221               ! in whether it uses level pool, machine learning, or persistence.
1222               ! Inflow, lateral inflow, water elevation, and the routing period
1223               ! are passed in. Updated outflow and water elevation are returned.
1224               call rt_domain(did)%reservoirs(i)%ptr%run(QLAKEIP(i), QLAKEI(i), &
1225                    QLLAKE(i), RESHT(i), QLAKEO(i), DTCT, rt_domain(did)%final_reservoir_type(i), &
1226                    rt_domain(did)%reservoir_assimilated_value(i), rt_domain(did)%reservoir_assimilated_source_file(i))
1227
1228               ! TODO: Encapsulate the lake state variable for water elevation (RESHT(i))
1229               !       inside the reservoir module, but it requires a redesign of the lake
1230               !       MPI communication model.
1231
1232               QLAKEIP(i) = QLAKEI(i)  !-- store total lake inflow for this timestep
1233
1234
1235 #ifdef MPP_LAND
1236             endif
1237 #endif
1238            enddo
1239 #ifdef MPP_LAND
1240     if(NLAKES .gt. 0) then
1241 !yw       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLLAKE,NLAKES,99)
1242 !yw       call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,RESHT,NLAKES,99)
1243 !yw      call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEO,NLAKES,99)
1244 !yw      call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
1245 !yw      call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEIP,NLAKES,99)
1246
1247          ! TODO: Change updateLake_grid calls below to updated distributed reservoir
1248          !       objects and not arrays as currently implemented.
1249          call updateLake_grid(QLLAKE, nlakes,lake_index)
1250          call updateLake_grid(RESHT,  nlakes,lake_index)
1251          call updateLake_grid(QLAKEO, nlakes,lake_index)
1252          call updateLake_grid(QLAKEI, nlakes,lake_index)
1253          call updateLake_grid(QLAKEIP,nlakes,lake_index)
1254     endif
1255 #endif
1256
1257
1258 #ifdef MPP_LAND
1259          do iyw = 1,yw_MPP_NLINKS
1260             i = nlinks_index(iyw)
1261 #else
1262          do i = 1, NLINKS                                                !--- compute volume and depth at each node
1263 #endif
1264             if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
1265                    QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
1266                       HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
1267                       CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
1268             endif
1269          enddo
1270
1271 #ifdef MPP_LAND
1272           call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
1273 #endif
1274
1275            KRT = KRT + 1                     !-- iterate on the timestep
1276            if(KRT .eq. DT_STEPS) exit crnt   !-- up to the maximum time in interval
1277
1278           end do crnt   !--- DTCT timestep of DT_STEPS
1279
1280            HLINK = HLINKTMP                 !-- update head based on final solution in timestep
1281            CVOL  = CVOLTMP                  !-- update volume
1282         else                                !-- no channel option apparently selected
1283          print *, "FATAL ERROR: no channel option selected"
1284          call hydro_stop("In drive_CHANNEL() - no channel option selected")
1285         endif
1286
1287 #ifdef HYDRO_D
1288          write(6,*) "finished call drive_CHANNEL"
1289 #endif
1290
1291         if (KT .eq. 1) KT = KT + 1
1292
1293
1294 end subroutine drive_CHANNEL
1295 ! ----------------------------------------------------------------
1296
1297 !-=======================================
1298      REAL FUNCTION AREAf(AREA,Bw,h,z)
1299      REAL :: AREA, Bw, z, h
1300        AREAf = (Bw+z*h)*h-AREA       !-- Flow area
1301      END FUNCTION AREAf
1302
1303 !-====critical depth function  ==========
1304      REAL FUNCTION CDf(Q,Bw,h,z)
1305      REAL :: Q, Bw, z, h
1306        if(h .le. 0) then
1307          print *, "FATAL ERROR: head is zero, will get division by zero error"
1308          call hydro_stop("In CDf() - head is zero")
1309        else
1310        CDf = (Q/((Bw+z*h)*h))/(sqrt(9.81*(((Bw+z*h)*h)/(Bw + 2.0*z*h)))) - 1.0  !--critical depth function
1311        endif
1312      END FUNCTION CDf
1313
1314 !=======find flow depth in channel with bisection Chapra pg. 131
1315     REAL FUNCTION HEAD(idx,AREA,Bw,z)  !-- find the water elevation given wetted area,
1316                                          !--bottom widith and side channel.. index was for debuggin
1317      REAL :: Bw,z,AREA,test
1318      REAL :: hl, hu, hr, hrold
1319      REAL :: fl, fr,error                !-- function evaluation
1320      INTEGER :: maxiter, idx
1321
1322      error = 1.0
1323      maxiter = 0
1324      hl = 0.00001   !-- minimum depth is small
1325      hu = 30.  !-- assume maximum depth is 30 meters
1326
1327     if (AREA .lt. 0.00001) then
1328      hr = 0.
1329     else
1330        ! JLM: .gt. "0" is somewhat alarming. We really should have a constant like zero_r4
1331       do while ((AREAf(AREA,BW,hl,z)*AREAf(AREA,BW,hu,z)) .gt. 0.0 .and. maxiter .lt. 100)
1332        !-- allows for larger , smaller heads
1333        if(AREA .lt. 1.) then
1334         hl=hl/2
1335        else
1336         hu = hu * 2
1337        endif
1338        maxiter = maxiter + 1
1339
1340       end do
1341
1342       maxiter =0
1343       hr = 0
1344       fl = AREAf(AREA,Bw,hl,z)
1345       do while (error .gt. 0.0001 .and. maxiter < 1000)
1346         hrold = hr
1347         hr = (hl+hu)/2
1348         fr =  AREAf(AREA,Bw,hr,z)
1349         maxiter = maxiter + 1
1350          if (hr .ne. 0) then
1351           error = abs((hr - hrold)/hr)
1352          endif
1353         test = fl * fr
1354          if (test.lt.0) then
1355            hu = hr
1356          elseif (test.gt.0) then
1357            hl=hr
1358            fl = fr
1359          else
1360            error = 0.0
1361          endif
1362       end do
1363      endif
1364      HEAD = hr
1365
1366 22   format("i,hl,hu,Area",i5,2x,f12.8,2x,f6.3,2x,f6.3,2x,f6.3,2x,f9.1,2x,i5)
1367
1368     END FUNCTION HEAD
1369 !=================================
1370      REAL FUNCTION MANNING(h1,n,Bw,Cs)
1371
1372      REAL :: Bw,h1,Cs,n
1373      REAL :: z, AREA,R,Kd
1374
1375      z = 1.0/Cs
1376      R = ((Bw + z*h1)*h1)/(Bw + 2.0*h1*sqrt(1.0 + z*z)) !-- Hyd Radius
1377      AREA = (Bw + z*h1)*h1        !-- Flow area
1378      Kd = (1.0/n) * (R**(2.0/3.0))*AREA     !-- convenyance
1379 #ifdef HYDRO_D
1380      print *,"head, kd",  h1,Kd
1381 #endif
1382      MANNING = Kd
1383
1384      END FUNCTION MANNING
1385
1386 !=======find flow depth in channel with bisection Chapra pg. 131
1387      REAL FUNCTION CRITICALDEPTH(lnk, Q, Bw, z)  !-- find the critical depth
1388      REAL :: Bw, z, Q, test
1389      REAL :: hl, hu, hr, hrold
1390      REAL :: fl, fr, error   !-- function evaluation
1391      INTEGER :: maxiter
1392      INTEGER :: lnk
1393
1394      error = 1.0
1395      maxiter = 0
1396      hl = 1.0e-5   !-- minimum depth is 0.00001 meters
1397 !    hu = 35.       !-- assume maximum  critical depth 25 m
1398      hu = 100.       !-- assume maximum  critical depth 25 m
1399
1400      if(CDf(Q, BW, hl, z) * CDf(Q, BW, hu, z) .gt. 0.0) then
1401       if(Q .gt. 0.001) then
1402 #ifdef HYDRO_D
1403         print *, "interval won't work to find CD of lnk ", lnk
1404         print *, "Q, hl, hu", Q, hl, hu
1405         print *, "cd lwr, upr", CDf(Q, BW, hl, z), CDf(Q, BW, hu, z)
1406         ! call hydro_stop("In CRITICALDEPTH()")
1407         CRITICALDEPTH = -9999
1408         return
1409 #endif
1410       else
1411         Q = 0.0
1412       endif
1413      endif
1414
1415      hr = 0.
1416      fl = CDf(Q, Bw, hl, z)
1417
1418      if (Q .eq. 0.0) then
1419        hr = 0.
1420      else
1421       do while (error .gt. 0.0001 .and. maxiter < 1000)
1422         hrold = hr
1423         hr = (hl+hu)/2.0
1424         fr =  CDf(Q, Bw, hr, z)
1425         maxiter = maxiter + 1
1426          if (hr .ne. 0.0) then
1427           error = abs((hr - hrold)/hr)
1428          endif
1429         test = fl * fr
1430          if (test .lt. 0) then
1431            hu = hr
1432          elseif (test .gt. 0) then
1433            hl=hr
1434            fl = fr
1435          else
1436            error = 0.0
1437          endif
1438
1439        end do
1440       endif
1441
1442      CRITICALDEPTH = hr
1443      END FUNCTION CRITICALDEPTH
1444
1445
1446      REAL FUNCTION SGNf(val)  !-- function to return the sign of a number
1447      REAL:: val
1448
1449      if (val .lt. 0) then
1450        SGNf= -1.
1451      elseif (val.gt.0) then
1452        SGNf= 1.
1453      else
1454        SGNf= 0.
1455      endif
1456
1457      END FUNCTION SGNf
1458 !================================================
1459
1460      REAL FUNCTION fnDX(qp,Tw,So,Ck,dx,dt) !-- find channel sub-length for MK method
1461      REAL    :: qp,Tw,So,Ck,dx, dt,test
1462      REAL    :: dxl, dxu, dxr, dxrold
1463      REAL    :: fl, fr, error
1464      REAL    :: X
1465      INTEGER :: maxiter
1466
1467      error = 1.0
1468      maxiter =0
1469      dxl = dx*0.9  !-- how to choose dxl???
1470      dxu = dx
1471      dxr=0
1472
1473      do while (fnDXCDT(qp,Tw,So,Ck,dxl,dt)*fnDXCDT(qp,Tw,So,Ck,dxu,dt) .gt. 0 &
1474                .and. dxl .gt. 10)  !-- don't let dxl get too small
1475       dxl = dxl/1.1
1476      end do
1477
1478
1479      fl = fnDXCDT(qp,Tw,So,Ck,dxl,dt)
1480      do while (error .gt. 0.0001 .and. maxiter < 1000)
1481         dxrold = dxr
1482         dxr = (dxl+dxu)/2
1483         fr =  fnDXCDT(qp,Tw,So,Ck,dxr,dt)
1484         maxiter = maxiter + 1
1485          if (dxr .ne. 0) then
1486           error = abs((dxr - dxrold)/dxr)
1487          endif
1488         test = fl * fr
1489          if (test.lt.0) then
1490            dxu = dxr
1491          elseif (test.gt.0) then
1492            dxl=dxr
1493            fl = fr
1494          else
1495            error = 0.0
1496          endif
1497       end do
1498      FnDX = dxr
1499
1500     END FUNCTION fnDX
1501 !================================================
1502      REAL FUNCTION fnDXCDT(qp,Tw,So,Ck,dx,dt) !-- function to help find sub-length for MK method
1503      REAL    :: qp,Tw,So,Ck,dx,dt,X
1504      REAL    :: c,b  !-- coefficients on dx/cdt log approximation function
1505
1506      c = 0.2407
1507      b = 1.16065
1508      X = 0.5-(qp/(2.0*Tw*So*Ck*dx))
1509      if (X .le. 0.0) then
1510       fnDXCDT = -1.0 !0.115
1511      else
1512       fnDXCDT = (dx/(Ck*dt)) - (c*LOG(X)+b)  !-- this function needs to converge to 0
1513      endif
1514      END FUNCTION fnDXCDT
1515 ! ----------------------------------------------------------------------
1516
1517     subroutine check_lake(unit,cd,lake_index,nlakes)
1518          use module_RT_data, only: rt_domain
1519          implicit none
1520          integer :: unit,nlakes,i,lake_index(nlakes)
1521          real cd(nlakes)
1522 #ifdef MPP_LAND
1523          call write_lake_real(cd,lake_index,nlakes)
1524 #endif
1525          write(unit,*) cd
1526           call flush(unit)
1527          return
1528     end subroutine check_lake
1529
1530     subroutine check_channel(unit,cd,did,nlinks)
1531          use module_RT_data, only: rt_domain
1532 #ifdef MPP_LAND
1533   USE module_mpp_land
1534 #endif
1535          implicit none
1536          integer :: unit,nlinks,i, did
1537          real cd(nlinks)
1538 #ifdef MPP_LAND
1539          real g_cd(rt_domain(did)%gnlinks)
1540          call write_chanel_real(cd,rt_domain(did)%map_l2g,rt_domain(did)%gnlinks,nlinks,g_cd)
1541          if(my_id .eq. IO_id) then
1542             write(unit,*) "rt_domain(did)%gnlinks = ",rt_domain(did)%gnlinks
1543            write(unit,*) g_cd
1544          endif
1545 #else
1546            write(unit,*) cd
1547 #endif
1548           call flush(unit)
1549           close(unit)
1550          return
1551     end subroutine check_channel
1552     subroutine smoth121(var,nlinks,maxv_p,from_node,to_node)
1553         implicit none
1554         integer,intent(in) ::  nlinks, maxv_p
1555         integer, intent(in), dimension(nlinks):: to_node
1556         integer, intent(in), dimension(nlinks):: from_node(nlinks,maxv_p)
1557         real, intent(inout), dimension(nlinks) :: var
1558         real, dimension(nlinks) :: vartmp
1559         integer :: i,j  , k, from,to
1560         integer :: plen
1561               vartmp = 0
1562               do i = 1, nlinks
1563                  to = to_node(i)
1564                  plen = from_node(i,1)
1565                  if(plen .gt. 1) then
1566                      do k = 1, plen-1
1567                          from = from_node(i,k+1)
1568                          if(to .gt. 0) then
1569                             vartmp(i) = vartmp(i)+0.25*(var(from)+2.*var(i)+var(to))
1570                          else
1571                             vartmp(i) = vartmp(i)+(2.*var(i)+var(from))/3.0
1572                          endif
1573                      end do
1574                      vartmp(i) = vartmp(i) /(plen-1)
1575                  else
1576                          if(to .gt. 0) then
1577                             vartmp(i) = vartmp(i)+(2.*var(i)+var(to)/3.0)
1578                          else
1579                             vartmp(i) = var(i)
1580                          endif
1581                  endif
1582               end do
1583               var = vartmp
1584         return
1585     end subroutine smoth121
1586
1587 !   SUBROUTINE drive_CHANNEL for NHDPLUS
1588 ! ------------------------------------------------
1589
1590      subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT,  &
1591         LAKEINFLORT, QSTRMVOLRT, TO_NODE, FROM_NODE, &
1592         TYPEL, ORDER, MAXORDER,   CH_LNKRT, &
1593         LAKE_MSKRT, DT, DTCT, DTRT_CH,MUSK, MUSX, QLINK, &
1594         CHANLEN, MannN, So, ChSSlp, Bw, Tw, Tw_CC, n_CC, &
1595         ChannK, RESHT, CVOL, QLAKEI, QLAKEO, LAKENODE, &
1596         QINFLOWBASE, CHANXI, CHANYJ, channel_option,  &
1597         nlinks,NLINKSL, LINKID, node_area, qout_gwsubbas, &
1598         LAKEIDA, LAKEIDM, NLAKES, LAKEIDX, &
1599 #ifdef MPP_LAND
1600         nlinks_index,mpp_nlinks,yw_mpp_nlinks, &
1601         LNLINKSL, &
1602         gtoNode,toNodeInd,nToNodeInd,   &
1603 #endif
1604          CH_LNKRT_SL, landRunOff  &
1605 #ifdef WRF_HYDRO_NUDGING
1606        , nudge &
1607 #endif
1608        , accSfcLatRunoff, accBucket                  &
1609        , qSfcLatRunoff,     qBucket                  &
1610        , QLateral, velocity, qloss                   &
1611        , HLINK                                       &
1612        , nsize , OVRTSWCRT, SUBRTSWCRT, channel_only, channelBucket_only, &
1613        channel_bypass)
1614
1615        use module_UDMAP, only: LNUMRSL, LUDRSL
1616        use config_base, only: nlst
1617
1618 #ifdef WRF_HYDRO_NUDGING
1619        use module_RT_data, only: rt_domain
1620        use module_stream_nudging,  only: setup_stream_nudging,               &
1621                                          nudge_term_all,                     &
1622                                          nudgeWAdvance,                      &
1623                                          nudge_apply_upstream_muskingumCunge
1624 #endif
1625
1626
1627        implicit none
1628
1629 ! -------- DECLARATIONS ------------------------
1630        integer, intent(IN) :: did, IXRT,JXRT,channel_option, OVRTSWCRT, SUBRTSWCRT
1631        integer, intent(IN) :: NLAKES, NLINKSL, nlinks
1632        integer, intent(INOUT) :: KT   ! flag of cold start (1) or continue run.
1633        real, intent(IN), dimension(IXRT,JXRT)    :: QSTRMVOLRT
1634        real, intent(IN), dimension(IXRT,JXRT)    :: LAKEINFLORT
1635        real, intent(IN), dimension(IXRT,JXRT)    :: QINFLOWBASE
1636        real, dimension(ixrt,jxrt) :: landRunOff
1637
1638        integer(kind=int64), intent(IN), dimension(IXRT,JXRT) :: CH_LNKRT
1639        integer(kind=int64), intent(IN), dimension(IXRT,JXRT) :: CH_LNKRT_SL
1640
1641        integer, intent(IN), dimension(IXRT,JXRT) :: LAKE_MSKRT
1642        integer, intent(IN), dimension(:)         :: ORDER, TYPEL !--link
1643        integer(kind=int64), intent(in), dimension(:)     :: TO_NODE, FROM_NODE
1644        integer, intent(IN), dimension(:)     :: CHANXI, CHANYJ
1645        real, intent(IN), dimension(:)        :: MUSK, MUSX
1646        real, intent(IN), dimension(:)        :: CHANLEN
1647        real, intent(IN), dimension(:)        :: So, MannN
1648        real, intent(IN), dimension(:)        :: ChSSlp,Bw,Tw  !--properties of nodes or links
1649        real, intent(IN), dimension(:)        :: Tw_CC, n_CC   ! properties of compound channel
1650        real, intent(IN), dimension(:)        :: ChannK  !--channel infiltration
1651        real                                      :: Km, X
1652        real , intent(INOUT), dimension(:,:)  :: QLINK
1653        real , intent(INOUT), dimension(:)    :: HLINK
1654
1655        logical, intent(in) :: channel_bypass
1656
1657 #ifdef WRF_HYDRO_NUDGING
1658        !! inout for applying previous nudge to upstream components of flow at gages
1659        real, intent(inout),    dimension(:)    :: nudge
1660 #endif
1661
1662        real, dimension(:), intent(inout)     :: QLateral !--lateral flow
1663        real, dimension(:), intent(out)       :: velocity, qloss
1664        real*8, dimension(:), intent(inout)     :: accSfcLatRunoff, accBucket
1665        real  , dimension(:), intent(out)     :: qSfcLatRunoff  ,   qBucket
1666
1667        real ,  dimension(NLINKSL,2) :: tmpQLINK
1668        real, intent(IN)                          :: DT    !-- model timestep
1669        real, intent(IN)                          :: DTRT_CH  !-- routing timestep
1670        real, intent(INOUT)                       :: DTCT
1671        real                                      :: minDTCT !BF minimum routing timestep
1672        integer, intent(IN)                       :: MAXORDER
1673        real , intent(IN), dimension(:)   :: node_area
1674
1675        !DJG GW-chan coupling variables...
1676        real, dimension(NLINKS)                   :: dzGwChanHead
1677        real, dimension(NLINKS)                   :: Q_GW_CHAN_FLUX     !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state...
1678        real, dimension(IXRT,JXRT)                :: ZWATTBLRT          !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state...
1679
1680        !-- lake params
1681        integer(kind=int64), intent(IN), dimension(:)    :: LAKEIDM  !-- NHDPLUS lakeid for lakes to be modeled
1682
1683        real, intent(INOUT), dimension(:)    :: RESHT    !-- reservoir height (m)
1684        real, intent(INOUT), dimension(:)    :: QLAKEI   !-- lake inflow (cms)
1685        real,                dimension(NLAKES)    :: QLAKEIP  !-- lake inflow previous timestep (cms)
1686        real, intent(INOUT), dimension(NLAKES)    :: QLAKEO   !-- outflow from lake used in diffusion scheme
1687
1688        integer(kind=int64), intent(IN), dimension(:)    :: LAKENODE !-- outflow from lake used in diffusion scheme
1689        integer(kind=int64), intent(IN), dimension(:)   :: LINKID   !--  id of channel elements for linked scheme
1690        integer(kind=int64), intent(IN), dimension(:)   :: LAKEIDA  !--  (don't need) NHDPLUS lakeid for all lakes in domain
1691        integer, intent(IN), dimension(:)   :: LAKEIDX  !--  the sequential index of the lakes id by com id
1692
1693        real, dimension(NLINKS)                   :: QSUM     !--mass bal of node
1694        real, dimension(NLAKES)                   :: QLLAKE   !-- lateral inflow to lake in diffusion scheme
1695        integer :: nsize
1696
1697        !-- Local Variables
1698        integer                      :: i,j,k,t,m,jj,ii,lakeid, kk,KRT,node, UDMP_OPT
1699        integer                      :: DT_STEPS               !-- number of timestep in routing
1700        real                         :: Qup,Quc                !--Q upstream Previous, Q Upstream Current, downstream Previous
1701        real                         :: bo                     !--critical depth, bnd outflow just for testing
1702
1703        real ,dimension(NLINKS)                          :: CD    !-- critical depth
1704        real, dimension(IXRT,JXRT)                       :: tmp
1705        real, dimension(nlinks)                          :: tmp2
1706        real, intent(INOUT), dimension(:)           :: CVOL
1707
1708 #ifdef MPP_LAND
1709        real*8,  dimension(LNLINKSL) :: LQLateral
1710        real*8,  dimension(LNLINKSL) :: tmpLQLateral
1711        real,  dimension(NLINKSL)    :: tmpQLateral
1712
1713        integer nlinks_index(:)
1714        integer  iyw, yw_mpp_nlinks, mpp_nlinks
1715        real     ywtmp(ixrt,jxrt)
1716        integer LNLINKSL
1717        integer, dimension(:)         ::  toNodeInd
1718        integer(kind=int64), dimension(:,:)       ::  gtoNode
1719        integer  :: nToNodeInd
1720        real, dimension(nToNodeInd,2) :: gQLINK
1721 #else
1722        real*8,  dimension(NLINKS) :: tmpLQLateral
1723        real,  dimension(NLINKSL) :: tmpQLateral
1724        real,  dimension(NLINKSL) :: LQLateral
1725 #endif
1726        integer flag
1727        integer, intent(in) :: channel_only, channelBucket_only
1728
1729        integer :: n, kk2, nt, nsteps  ! tmp
1730        real, intent(in), dimension(:) :: qout_gwsubbas
1731        real, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI, tmpRESHT
1732        integer, allocatable, dimension(:) :: tmpFinalResType
1733        real, allocatable,dimension(:) :: tmpAssimilatedValue
1734        character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile
1735
1736 #ifdef MPP_LAND
1737        if(my_id .eq. io_id) then
1738 #endif
1739             allocate(tmpQLAKEO(NLAKES))
1740             allocate(tmpQLAKEI(NLAKES))
1741             allocate(tmpRESHT(NLAKES))
1742             allocate(tmpFinalResType(nlakes))
1743 #ifdef MPP_LAND
1744         endif
1745 #endif
1746
1747         QLAKEIP = 0
1748         CD = 0
1749         node = 1
1750         QSUM     = 0
1751         QLLAKE   = 0
1752         dzGwChanHead = 0.
1753         nsteps = (DT+0.5)/DTRT_CH
1754
1755 #ifdef WRF_HYDRO_NUDGING
1756          !! Initialize nudging for the current timestep.
1757          !! This establishes the data structure used to solve the nudges.
1758          call setup_stream_nudging(0)  !! always zero b/c at beginning of hydro timestep
1759 #endif /* WRF_HYDRO_NUDGING */
1760
1761 !---------------------------------------------
1762 if(channel_only .eq. 1 .or. channelBucket_only .eq. 1) then
1763 #ifdef HYDRO_D
1764    write(6,*)  "channel_only or channelBucket_only is not zero. Special flux calculations."
1765    call flush(6)
1766 #endif /* HYDRO_D */
1767
1768 !   if(nlst_rt(1)%output_channelBucket_influx .eq. 1 .or. &
1769 !        nlst_rt(1)%output_channelBucket_influx .eq. 2       ) &
1770 !        !! qScfLatRunoff = qLateral - qBucket
1771 !        qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL) - qout_gwsubbas(1:NLINKSL)
1772
1773    if(nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1774       nlst(1)%output_channelBucket_influx .eq. 2       ) then
1775
1776       if(channel_only .eq. 1) &
1777         !! qScfLatRunoff = qLateral - qBucket
1778         qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL) - qout_gwsubbas(1:NLINKSL)
1779
1780       if(channelBucket_only .eq. 1) &
1781         !! qScfLatRunoff = qLateral - qBucket
1782         qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL)
1783
1784    end if
1785
1786    if(nlst(1)%output_channelBucket_influx .eq. 3) &
1787         accSfcLatRunoff(1:NLINKSL) = qSfcLatRunoff * DT
1788
1789 else
1790
1791    QLateral = 0 !! the variable solved in this section. Channel only knows this already.
1792    LQLateral = 0          !-- initial lateral flow to 0 for this reach
1793
1794    tmpQLateral = 0  !! WHY DOES THIS tmp variable EXIST?? Only for accumulations??
1795    tmpLQLateral = 0
1796
1797    ! NHDPLUS maping
1798    if(OVRTSWCRT .eq. 0)      then
1799       do k = 1, LNUMRSL
1800          ! get from land grid runoff
1801          do m = 1, LUDRSL(k)%ncell
1802             ii =  LUDRSL(k)%cell_i(m)
1803             jj =  LUDRSL(k)%cell_j(m)
1804             LQLateral(k) = LQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 &
1805                  *LUDRSL(k)%cellArea(m)/DT
1806             tmpLQLateral(k) = tmpLQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 &
1807                  *LUDRSL(k)%cellArea(m)/DT
1808          end do
1809       end do
1810
1811 #ifdef MPP_LAND
1812       call updateLinkV(tmpLQLateral, tmpQLateral)
1813 #endif
1814       if(NLINKSL .gt. 0) then
1815          if (nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1816              nlst(1)%output_channelBucket_influx .eq. 2      ) &
1817                qSfcLatRunoff(1:NLINKSL) = tmpQLateral(1:NLINKSL)
1818          if (nlst(1)%output_channelBucket_influx .eq. 3) &
1819               accSfcLatRunoff(1:NLINKSL) = accSfcLatRunoff(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
1820       endif
1821       tmpLQLateral = 0  !! JLM:: These lines imply that squeege runoff does not count towards
1822       tmpQLateral = 0   !! JLM:: accumulated runoff to be output but it does for internal QLateral?
1823       !! JLM: But then the next accumulation is added to the amt before zeroing? result
1824       !! JLM: should be identical to LQLateral.... I'm totally mystified.
1825    endif
1826
1827    !! JLM:: if ovrtswcrt=0 and subrtswcrt=1, then this accumulation is calculated twice for LQLateral???
1828    !! This impiles that if overland routing is off and subsurface routing is on, that
1829    !! qstrmvolrt represents only the subsurface contribution to the channel.
1830    if(OVRTSWCRT .ne. 0 .or. SUBRTSWCRT .ne. 0 ) then
1831       do k = 1, LNUMRSL
1832          ! get from channel grid
1833          do m = 1, LUDRSL(k)%ngrids
1834             ii =  LUDRSL(k)%grid_i(m)
1835             jj =  LUDRSL(k)%grid_j(m)
1836             LQLateral(k) = LQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 &
1837                  *LUDRSL(k)%nodeArea(m)/DT
1838             tmpLQLateral(k) = tmpLQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 &
1839                  *LUDRSL(k)%nodeArea(m)/DT
1840          end do
1841       end do
1842
1843 #ifdef MPP_LAND
1844       call updateLinkV(tmpLQLateral, tmpQLateral)
1845 #endif
1846
1847       !! JLM:: again why output in this conditional ?? why not just output QLateral
1848       !! after this section ????
1849       if(NLINKSL .gt. 0) then
1850          if(nlst(1)%output_channelBucket_influx .eq. 1 .OR. &
1851             nlst(1)%output_channelBucket_influx .eq. 2       ) &
1852               qSfcLatRunoff(1:NLINKSL) = tmpQLateral(1:NLINKSL)
1853          if(nlst(1)%output_channelBucket_influx .eq. 3) &
1854               accSfcLatRunoff(1:NLINKSL) = accSfcLatRunoff(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
1855       end if
1856
1857    endif
1858
1859 #ifdef MPP_LAND
1860    call updateLinkV(LQLateral, QLateral(1:NLINKSL))
1861 #else
1862    call hydro_stop("fatal error: NHDPlus only works for parallel now.")
1863    QLateral = LQLateral
1864 #endif
1865 endif !! (channel_only .eq. 1 .or. channelBucket_only .eq. 1) then; else; endif
1866
1867
1868 !---------------------------------------------
1869 !! If not running channelOnly, here is where the bucket model is picked up
1870 if(channel_only .eq. 1) then
1871 #ifdef HYDRO_D
1872    write(6,*)  "channel_only is not zero. No bucket."
1873    call flush(6)
1874 #endif /* HYDRO_D */
1875 else
1876    !! REQUIRE BUCKET MODEL ON HERE?
1877    if(NLINKSL .gt. 0) QLateral(1:NLINKSL) = QLateral(1:NLINKSL) + qout_gwsubbas(1:NLINKSL)
1878 endif  !! if(channel_only .eq. 1) then; else; endif
1879
1880 if(nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1881    nlst(1)%output_channelBucket_influx .eq. 2       ) &
1882       qBucket(1:NLINKSL) = qout_gwsubbas(1:NLINKSL)
1883
1884 if(nlst(1)%output_channelBucket_influx .eq. 3) &
1885      accBucket(1:NLINKSL) = accBucket(1:NLINKSL) + qout_gwsubbas(1:NLINKSL) * DT
1886
1887 ! Skip this section if we are NOT running any actual channel routing
1888 if (.not. channel_bypass) then
1889
1890 !---------------------------------------------
1891 !       QLateral = QLateral / nsteps
1892 do nt = 1, nsteps
1893
1894 #ifdef MPP_LAND
1895    gQLINK = 0
1896    call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
1897    call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1))
1898    !---------- route other reaches, with upstream inflow
1899 #endif
1900
1901    tmpQlink = 0
1902 #ifdef MPP_LAND
1903    if(my_id .eq. io_id) then
1904 #endif
1905       tmpQLAKEO = QLAKEO
1906       tmpQLAKEI = QLAKEI
1907       tmpRESHT = RESHT
1908       tmpFinalResType = rt_domain(did)%final_reservoir_type
1909       tmpAssimilatedValue = rt_domain(did)%reservoir_assimilated_value
1910       tmpAssimilatedSourceFile = rt_domain(did)%reservoir_assimilated_source_file
1911
1912 #ifdef MPP_LAND
1913    endif
1914 #endif
1915
1916    do k = 1,NLINKSL
1917
1918       Quc  = 0
1919       Qup  = 0
1920
1921       !process as standard link or a lake inflow link, or lake outflow link
1922       ! link flowing out of lake, accumulate all the inflows with the revised TO_NODEs
1923       ! TYPEL = -999 stnd; TYPEL=1 outflow from lake; TYPEL = 3 inflow to a lake
1924
1925       if(TYPEL(k) .ne. 2) then ! don't process internal lake links only
1926
1927 #ifdef MPP_LAND
1928          !using mapping index
1929          do n = 1, gtoNODE(k,1)
1930             m = gtoNODE(k,n+1)
1931             !! JLM - I think gQLINK(,2) is actually previous. Global array never sees current. Seeing
1932             !! current would require global communication at the end of each loop through k
1933             !! (=kth reach). Additionally, how do you synchronize to make sure the upstream are all
1934             !! done before doing the downstream?
1935             if(gQLINK(m,2) .gt. 0)   Quc = Quc + gQLINK(m,2)  !--accum of upstream inflow of current timestep (2)
1936             if(gQLINK(m,1) .gt. 0)   Qup = Qup + gQLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
1937          end do ! do i
1938 #else
1939          do m = 1, NLINKSL
1940             if (LINKID(k) .eq. TO_NODE(m)) then
1941                Quc = Quc + QLINK(m,2)  !--accum of upstream inflow of current timestep (2)
1942                Qup = Qup + QLINK(m,1)  !--accum of upstream inflow of previous timestep (1)
1943             endif
1944          end do ! do m
1945 #endif
1946       endif !note that we won't process type 2 links, since they are internal to a lake
1947
1948
1949       !yw ### process each link k,
1950       !       There is a situation that different k point to the same LAKEIDX
1951       !        if(TYPEL(k) .eq. 1 .and. LAKEIDX(k) .gt. 0) then   !--link is a reservoir
1952       if(TYPEL(k) .eq. 1 ) then   !--link is a reservoir
1953
1954          lakeid = LAKEIDX(k)
1955          if(lakeid .ge. 0) then
1956
1957             ! Calls run for a single reservoir depending on its type as
1958             ! in whether it uses level pool, machine learning, or persistence.
1959             ! Inflow, lateral inflow, water elevation, and the routing period
1960             ! are passed in. Updated outflow and water elevation are returned.
1961             call rt_domain(did)%reservoirs(lakeid)%ptr%run(Qup, Quc, 0.0, &
1962                  RESHT(lakeid), tmpQLINK(k,2), DTRT_CH, rt_domain(did)%final_reservoir_type(lakeid), &
1963                  rt_domain(did)%reservoir_assimilated_value(lakeid), rt_domain(did)%reservoir_assimilated_source_file(lakeid))
1964
1965             ! TODO: Encapsulate the lake state variable for water elevation (RESHT(lakeid))
1966             !       inside the reservoir module, but it requires a redesign of the lake
1967             !       MPI communication model.
1968
1969             QLAKEO(lakeid)  = tmpQLINK(k,2) !save outflow to lake
1970             QLAKEI(lakeid)  = Quc           !save inflow to lake
1971          endif
1972 105      continue
1973
1974
1975       elseif (channel_option .eq. 1) then  !muskingum routing
1976          Km = MUSK(k)
1977          X = MUSX(k)
1978          tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow
1979
1980       elseif (channel_option .eq. 2) then ! muskingum cunge, don't process internal lake nodes TYP=2
1981          !              tmpQLINK(k,2) = MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), &
1982          !                  QLateral(k), DTRT_CH, So(k),  CHANLEN(k), &
1983          !                  MannN(k), ChSSlp(k), Bw(k), Tw(k) )
1984
1985 #ifdef WRF_HYDRO_NUDGING
1986          call nudge_apply_upstream_muskingumCunge( Qup,  Quc,  nudge(k),  k )
1987 #endif
1988
1989          call SUBMUSKINGCUNGE(&
1990               tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k),     Qup,        Quc, QLINK(k,1), &
1991               QLateral(k),   DTRT_CH,     So(k), CHANLEN(k),                  &
1992               MannN(k),      ChSSlp(k),   Bw(k), Tw(k), Tw_CC(k), n_CC(k), HLINK(k), ChannK(k) )
1993
1994       else
1995 #ifdef HYDRO_D
1996          print *, " no channel option selected"
1997 #endif
1998          call hydro_stop("drive_CHANNEL")
1999       endif
2000
2001    end do  !--k links
2002
2003
2004 #ifdef MPP_LAND
2005    call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO)
2006    call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI)
2007    call updateLake_seq(RESHT,nlakes,tmpRESHT)
2008    call updateLake_seqInt(rt_domain(did)%final_reservoir_type, nlakes, tmpFinalResType)
2009    call updateLake_seq(rt_domain(did)%reservoir_assimilated_value, nlakes, tmpAssimilatedValue)
2010    !call updateLake_seq_char(rt_domain(did)%reservoir_assimilated_source_file, nlakes, tmpAssimilatedSourceFile)
2011 #endif
2012
2013    do k = 1, NLINKSL !tmpQLINK?
2014       if(TYPEL(k) .ne. 2) then   !only the internal lake nodes don't have info.. but need to save QLINK of lake out too
2015          QLINK(k,2) = tmpQLINK(k,2)
2016       endif
2017       QLINK(k,1) = QLINK(k,2)    !assigng link flow of current to be previous for next time step
2018    end do
2019
2020
2021 #ifdef WRF_HYDRO_NUDGING
2022    if(.not. nudgeWAdvance) call nudge_term_all(qlink, nudge, int(nt*dtrt_ch))
2023 #endif /* WRF_HYDRO_NUDGING */
2024
2025
2026 !#ifdef HYDRO_D
2027 !   print *, "END OF ALL REACHES...",KRT,DT_STEPS
2028 !#endif
2029
2030 end do  ! nsteps
2031
2032 endif ! channel_bypass
2033
2034 if (KT .eq. 1) KT = KT + 1
2035
2036 #ifdef MPP_LAND
2037 if(my_id .eq. io_id)      then
2038    if(allocated(tmpQLAKEO))  deallocate(tmpQLAKEO)
2039    if(allocated(tmpQLAKEI))  deallocate(tmpQLAKEI)
2040    if(allocated(tmpRESHT))  deallocate(tmpRESHT)
2041    if(allocated(tmpFinalResType)) deallocate(tmpFinalResType)
2042 endif
2043 #endif
2044
2045 if (KT .eq. 1) KT = KT + 1  ! redundant?
2046
2047 end subroutine drive_CHANNEL_RSL
2048
2049 ! ----------------------------------------------------------------
2050
2051 end module module_channel_routing
2052
2053 !! Is this outside the module scope on purpose?
2054 #ifdef MPP_LAND
2055  subroutine checkReach(ii,  inVar)
2056    use module_mpp_land
2057    use module_RT_data, only: rt_domain
2058    use MODULE_mpp_ReachLS, only : updatelinkv,                   &
2059                                  ReachLS_write_io, gbcastvalue, &
2060                                  gbcastreal2
2061    implicit none
2062    integer :: ii
2063    real,dimension(rt_domain(1)%nlinksl) :: inVar
2064    real:: g_var(rt_domain(1)%gnlinksl)
2065    call ReachLS_write_io(inVar, g_var)
2066    if(my_id .eq. io_id) then
2067       write(ii,*) g_var
2068       call flush(ii)
2069    endif
2070  end subroutine checkReach
2071 #endif