| blob | ff26557e4c9bd8b87673b102e476f37dd684f8e3 |
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>
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
29 use module_reservoir, only: reservoir
30 use module_RT_data, only: rt_domain
31 use module_hydro_stop, only: HYDRO_stop
33 use iso_fortran_env, only: int64
35 #endif
36 implicit none
38 contains
40 ! ------------------------------------------------
41 ! FUNCTION MUSKING
42 ! ------------------------------------------------
43 real function MUSKING(idx,qup,quc,qdp,dt,Km,X)
45 implicit none
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
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)
65 ! ----------------------------------------------------------------
66 end function MUSKING
67 ! ----------------------------------------------------------------
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
94 ! added by Wei Yu for bad data.
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
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
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
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
128 sgn = SGNf(Sf) !-- establish sign
130 w = 0.5*(sgn + 1.) !-- compute upstream or downstream weighting
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
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
141 Kf = (1.0-w)*Kd + w*Ku !-- conveyance
142 DIFFUSION = Kf * sqrt(abs(Sf))*sgn
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)
147 end function DIFFUSION
148 ! ----------------------------------------------------------------
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 )
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
160 IMPLICIT NONE
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
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)
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
205 !-- additional channel infiltration parameters
206 REAL :: WPk ! KINEROS2 modified wetted perimeter
207 REAL :: modK ! KINEROS2 constant (set to 0.15, line )
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
214 !TML-- Hard code KINEROS2 Modification Parameter
215 !TML-- Set to 0.15, consistent with USDA-ARS
216 modK = 0.15
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
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
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
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
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
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
255 110 continue
257 Qj_0 = 0.0 !- initial flow of lower interval
258 iter = 0
260 do while (rerror .gt. 0.01 .and. aerror .ge. mindepth .and. iter .le. maxiter)
262 WPC = 0.0
263 AREAC = 0.0
265 !----- lower interval --------------------
267 Twl = Bw + 2.0*z*h_0 !--top surface water width of the channel inflow
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
286 endif
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
303 if(Ck .gt. 0.000000) then
304 Km = max(dt,dx/Ck)
305 else
306 Km = dt
307 endif
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
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
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
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.
347 !--upper interval -----------
348 WPC = 0.0
349 AREAC = 0.0
351 Twl = Bw + 2.0*z*h !--top width of the channel inflow
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
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
386 if(Ck .gt. 0.0) then
387 Km = max(dt,dx/Ck)
388 else
389 Km = dt
390 endif
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
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
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
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
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
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
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
449 ! if(idx .eq. 6276407) then
450 ! print*, "err,itr,hs", rerror, iter, depth, h_0, h, h_1
451 ! endif
453 h_0 = max(0.0,h)
454 h = max(0.0,h_1)
455 iter = iter + 1
457 if( h .lt. mindepth) then ! exit loop if depth is very small
458 goto 111
459 endif
461 end do
463 ! if(idx .eq. 6276407) then
464 ! print*, "id,itr,err,h", idx, iter, rerror, h
465 ! endif
467 111 continue
469 if(iter .ge. maxiter) then
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
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
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
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
514 endif
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
521 else ! no flow to route
522 qdc = 0.0
523 depth = 0.0
524 endif
526 qloss = (ChannK * dx * WPk) ! TML -- Compute qloss to pass
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
538 ! ----------------------------------------------------------------
539 END SUBROUTINE SUBMUSKINGCUNGE
540 ! ----------------------------------------------------------------
542 ! ------------------------------------------------
543 ! FUNCTION KINEMATIC
544 ! ------------------------------------------------
545 REAL FUNCTION KINEMATIC()
547 IMPLICIT NONE
549 ! -------- DECLARATIONS -----------------------
551 ! REAL, INTENT(OUT), DIMENSION(IXRT,JXRT) :: OVRGH
553 KINEMATIC = 1
554 !----------------------------------------------------------------
555 END FUNCTION KINEMATIC
556 !----------------------------------------------------------------
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, LLINKID &
576 , gtoNode,toNodeInd,nToNodeInd &
577 #endif
578 , CH_LNKRT_SL &
579 ,gwBaseSwCRT, gwHead, qgw_chanrt, gwChanCondSw, gwChanCondConstIn, &
580 gwChanCondConstOut, velocity, qloss)
583 IMPLICIT NONE
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
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
600 real , dimension(ixrt,jxrt):: latval,lonval
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
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(IXRT,JXRT) :: gwHead !DJG !!! groundwater head from Fersch-2d gw implementation...units (m ASL)
637 REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: qgw_chanrt !DJG !!! Channel-gw flux as used in Fersch 2d gw implementation...units (m^3/s)...Change 'INTENT' to 'OUT' when ready to update groundwater state...
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
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
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
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 integer(kind=int64), dimension(LNLINKSL) :: LLINKID
675 real*8, dimension(LNLINKSL) :: LQLateral
676 ! real*4, dimension(LNLINKSL) :: LQLateral
677 integer, dimension(:) :: toNodeInd
678 integer(kind=int64), dimension(:,:) :: gtoNode
679 integer :: nToNodeInd
680 real, dimension(nToNodeInd,2) :: gQLINK
681 #else
682 real*8, dimension(NLINKS) :: LQLateral !--lateral flow
683 #endif
684 integer flag
686 integer :: n, kk2, nt, nsteps ! tmp
688 QLAKEIP = 0
689 QLAKEI8 = 0
690 HLINKTMP = 0
691 CVOLTMP = 0
692 CD = 0
693 node = 1
694 QLateral = 0
695 QSUM = 0
696 QLLAKE = 0
699 !yw print *, "DRIVE_channel,option,nlinkl,nlinks!!", channel_option,NLINKSL,NLINKS
700 ! print *, "DRIVE_channel, RESHT", RESHT
703 dzGwChanHead = 0.
705 IF(channel_option .ne. 3) then !--muskingum methods ROUTE ON DT timestep, not DTRT!!
707 nsteps = (DT+0.5)/DTRT_CH
709 #ifdef MPP_LAND
710 LQLateral = 0 !-- initial lateral flow to 0 for this reach
711 DO iyw = 1,yw_MPP_NLINKS
712 jj = nlinks_index(iyw)
713 !--------river grid points, convert depth in mm to rate across reach in m^3/sec
714 if( .not. ( (CHANXI(jj) .eq. 1 .and. left_id .ge. 0) .or. &
715 (CHANXI(jj) .eq. ixrt .and. right_id .ge. 0) .or. &
716 (CHANYJ(jj) .eq. 1 .and. down_id .ge. 0) .or. &
717 (CHANYJ(jj) .eq. jxrt .and. up_id .ge. 0) &
718 ) ) then
719 if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0) then
720 k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
721 LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
722 *node_area(jj)/DT)
723 elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
724 k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
725 LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
726 *node_area(jj)/DT)
727 endif
728 endif
729 end do ! jj
732 ! assign LQLATERAL to QLATERAL
733 call updateLinkV(LQLateral, QLateral(1:NLINKSL))
735 #else
737 LQLateral = 0 !-- initial lateral flow to 0 for this reach
738 do jj = 1, NLINKS
739 !--------river grid points, convert depth in mm to rate across reach in m^3/sec
741 if (CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj)) .gt. 0 ) then
742 k = CH_LNKRT_SL(CHANXI(jj),CHANYJ(jj))
743 LQLateral(k) = LQLateral(k)+((QSTRMVOLRT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
744 *node_area(jj)/DT)
745 elseif ( (LAKE_MSKRT(CHANXI(jj),CHANYJ(jj)) .gt. 0)) then !-lake grid
746 k = LAKE_MSKRT(CHANXI(jj),CHANYJ(jj))
747 LQLateral(k) = LQLateral(k) +((LAKEINFLORT(CHANXI(jj),CHANYJ(jj))+QINFLOWBASE(CHANXI(jj),CHANYJ(jj)))/1000 &
748 *node_area(jj)/DT)
749 endif
751 end do ! jj
752 QLateral = LQLateral
754 #endif
756 ! QLateral = QLateral / nsteps
758 do nt = 1, nsteps
760 !Per Yates, this check is not needed. Commenting out for now.
761 !---------- route order 1 reaches which have no upstream inflow
762 ! do k=1, NLINKSL
763 ! if (ORDER(k) .eq. 1) then !-- first order stream has no headflow
766 ! if(TYPEL(k) .eq. 1) then !-- level pool route of reservoir
767 ! !CALL LEVELPOOL(1,0.0, 0.0, qd, QLINK(k,2), QLateral(k), &
768 ! ! DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), &
769 ! ! WEIRC(k), WEIRL(k), ORIFICEE(i), ORIFICEC(k), ORIFICEA(k) )
770 ! elseif (channel_option .eq. 1) then
771 ! Km = MUSK(k)
772 ! X = MUSX(k)
773 ! QLINK(k,2) = MUSKING(k,0.0, QLateral(k), QLINK(k,1), DTRT_CH, Km, X) !--current outflow
774 ! elseif (channel_option .eq. 2) then !-- upstream is assumed constant initial condition
776 ! call SUBMUSKINGCUNGE(QLINK(k,2), velocity(k), k, &
777 ! 0.0,0.0, QLINK(k,1), QLateral(k), DTRT_CH, So(k), &
778 ! CHANLEN(k), MannN(k), ChSSlp(k), Bw(k), Tw(k) )
780 ! else
781 ! print *, "FATAL ERROR: No channel option selected"
782 ! call hydro_stop("In drive_CHANNEL() -No channel option selected ")
783 ! endif
784 ! endif
785 ! end do
787 #ifdef MPP_LAND
788 gQLINK = 0.000
789 call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
790 call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1))
791 #endif
793 !---------- route other reaches, with upstream inflow
794 tmpQlink = 0.0
795 do k = 1,NLINKSL
796 ! if (ORDER(k) .gt. 1 ) then !-- exclude first order stream
797 Quc = 0.0
798 Qup = 0.0
800 #ifdef MPP_LAND
801 !using mapping index
802 do n = 1, gtoNODE(k,1)
803 m = gtoNODE(k,n+1)
804 !yw if (LINKID(k) .eq. m) then
805 Quc = Quc + gQLINK(m,2) !--accum of upstream inflow of current timestep (2)
806 Qup = Qup + gQLINK(m,1) !--accum of upstream inflow of previous timestep (1)
808 ! if(LINKID(k) .eq. 3259 .or. LINKID(k) .eq. 3316 .or. LINKID(k) .eq. 3219) then
809 ! write(6,*) "id,Uc,Up",LINKID(k),Quc,Qup
810 ! call flush(6)
811 ! endif
813 !yw endif
814 end do ! do i
816 #else
817 do m = 1, NLINKSL
818 if (LINKID(k) .eq. TO_NODE(m)) then
819 Quc = Quc + QLINK(m,2) !--accum of upstream inflow of current timestep (2)
820 Qup = Qup + QLINK(m,1) !--accum of upstream inflow of previous timestep (1)
821 endif
822 end do ! do m
823 #endif
825 if(TYPEL(k) .eq. 1) then !--link is a reservoir
827 ! CALL LEVELPOOL(1,QLINK(k,1), Qup, QLINK(k,1), QLINK(k,2), &
828 ! QLateral(k), DT, RESHT(k), HRZAREA(k), LAKEMAXH(k), &
829 ! WEIRC(k), WEIRL(k),ORIFICEE(k), ORIFICEC(k), ORIFICEA(k))
831 elseif (channel_option .eq. 1) then !muskingum routing
832 Km = MUSK(k)
833 X = MUSX(k)
834 tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow
835 elseif (channel_option .eq. 2) then ! muskingum cunge
837 call SUBMUSKINGCUNGE(tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k), &
838 Qup,Quc, QLINK(k,1), QLateral(k), DTRT_CH, So(k), &
839 CHANLEN(k), MannN(k), ChSSlp(k), Bw(k), Tw(k),Tw_CC(k), n_CC(k), HLINK(k), ChannK(k) )
841 else
842 print *, "FATAL ERROR: no channel option selected"
843 call hydro_stop("In drive_CHANNEL() - no channel option selected")
844 endif
845 ! endif !!! order(1) .ne. 1
846 end do !--k links
848 !yw check
849 ! gQLINK = 0.0
850 ! call ReachLS_write_io(tmpQLINK(:,2), gQLINK(:,2))
851 ! call ReachLS_write_io(tmpQLINK(:,1), gQLINK(:,1))
852 ! write(6,*) " io_id = ", io_id
853 ! if(my_id .eq. io_id) then
854 ! write(71,*) gQLINK(:,1)
855 ! call flush(71)
856 ! call flush(72)
857 ! endif
859 do k = 1, NLINKSL
860 if(TYPEL(k) .ne. 1) then
861 QLINK(k,2) = tmpQLINK(k,2)
862 endif
863 QLINK(k,1) = QLINK(k,2) !assing link flow of current to be previous for next time step
864 end do
866 end do ! nsteps
868 #ifdef HYDRO_D
869 print *, "END OF ALL REACHES...",KRT,DT_STEPS
870 #endif
872 ! END DO !-- krt timestep for muksingumcunge routing
874 elseif(channel_option .eq. 3) then !--- route using the diffusion scheme on nodes not links
876 #ifdef MPP_LAND
877 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINK,NLINKS,99)
878 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOL,NLINKS,99)
879 #endif
881 KRT = 0 !-- initialize the time counter
882 minDTCT = 0.01 ! define minimum routing sub-timestep (s), simulation will end with smaller timestep
883 DTCT = min(max(DTCT*2.0, minDTCT),DTRT_CH)
885 HLINKTMP = HLINK !-- temporary storage of the water elevations (m)
886 CVOLTMP = CVOL !-- temporary storage of the volume of water in channel (m^3)
887 QLAKEIP = QLAKEI !-- temporary lake inflow from previous timestep (cms)
889 ! call check_channel(77,HLINKTMP,1,nlinks)
890 ! call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,ZELEV,NLINKS,99)
891 crnt: DO !-- loop on the courant condition
892 QSUM = 0. !-- initialize the total flow out of each cell to zero
893 QSUM8 = 0. !-- initialize the total flow out of each cell to zero
894 QLAKEI8 = 0. !-- set the lake inflow as zero
895 QLAKEI = 0. !-- set the lake inflow as zero
896 QLLAKE = 0. !-- initialize each lake's lateral inflow to zero
897 QLLAKE8 = 0. !-- initialize each lake's lateral inflow to zero
898 DT_STEPS = INT(DT/DTCT) !-- fix the timestep
899 QLateral = 0.
900 !DJG GW-chan coupling variables...
901 if(gwBaseSwCRT == 3) then
902 Q_GW_CHAN_FLUX = 0.
903 qgw_chanrt = 0.
904 end if
906 ! ZWATTBLRT=1.0 !--HARDWIRE, remove this and pass in from subsfc/gw routing routines...
909 !-- vectorize
910 !---------------------
911 #ifdef MPP_LAND
912 DO iyw = 1,yw_MPP_NLINKS
913 i = nlinks_index(iyw)
914 #else
915 DO i = 1,NLINKS
916 #endif
918 if(node_area(i) .eq. 0) then
919 write(6,*) "FATAL ERROR: node_area(i) is zero. i=", i
920 call hydro_stop("In drive_CHANNEL() - Error node_area")
921 endif
925 nodeType:if((CH_NETRT(CHANXI(i), CHANYJ(i) ) .eq. 0) .and. &
926 (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .lt.0) ) then !--a reg. node
928 gwOption: if(gwBaseSwCRT == 3) then
930 ! determine potential gradient between groundwater head and channel stage
931 ! units in (m)
932 dzGwChanHead(i) = gwHead(CHANXI(i),CHANYJ(i)) - (HLINK(i)+ZELEV(i))
934 if(gwChanCondSw .eq. 0) then
936 qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
938 else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) > 0) then
940 ! channel bed interface, units in (m^3/s), flux into channel...
941 ! BF todo: consider channel width
942 qgw_chanrt(CHANXI(i),CHANYJ(i)) = gwChanCondConstIn * dzGwChanHead(i) &
943 * CHANLEN(i) * 2.
945 else if(gwChanCondSw .eq. 1 .and. dzGwChanHead(i) < 0) then
947 ! channel bed interface, units in (m^3/s), flux out of channel...
948 ! BF todo: consider channel width
949 qgw_chanrt(CHANXI(i),CHANYJ(i)) = max(-0.005, gwChanCondConstOut * dzGwChanHead(i) &
950 * CHANLEN(i) * 2.)
951 ! else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then TBD: exponential dependency
952 ! else if(gwChanCondSw .eq. 2 .and. dzGwChanHead(i) > 0) then
954 else
956 qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
958 end if
960 Q_GW_CHAN_FLUX(i) = qgw_chanrt(CHANXI(i),CHANYJ(i))
961 ! if ( i .eq. 1001 ) then
962 ! print *, Q_GW_CHAN_FLUX(i), dzGwChanHead(i), ELRT(CHANXI(i),CHANYJ(i)), HLINK(i), ZELEV(i)
963 ! end if
964 ! if ( Q_GW_CHAN_FLUX(i) .lt. 0. ) then !-- temporary hardwire for only allowing flux into channel...REMOVE later...
965 ! Q_GW_CHAN_FLUX(i) = 0.
966 ! qgw_chanrt(CHANXI(i),CHANYJ(i)) = 0.
967 ! end if
969 else
970 Q_GW_CHAN_FLUX(i) = 0.
971 end if gwOption
974 QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) = &
975 !DJG awaiting gw-channel exchg... Q_GW_CHAN_FLUX(i)+& ...obsolete-> ((QSUBRT(CHANXI(i),CHANYJ(i))+&
976 Q_GW_CHAN_FLUX(i)+&
977 ((QSTRMVOLRT(CHANXI(i),CHANYJ(i))+&
978 QINFLOWBASE(CHANXI(i),CHANYJ(i))) &
979 /DT_STEPS*node_area(i)/1000/DTCT)
980 if((QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))).lt.0.) .and. (gwChanCondSw == 0)) then
981 #ifdef HYDRO_D
982 print*, "i, CHANXI(i),CHANYJ(i) = ", i, CHANXI(i),CHANYJ(i)
983 print *, "NEGATIVE Lat inflow...",QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))), &
984 QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)), &
985 QINFLOWBASE(CHANXI(i),CHANYJ(i))
986 #endif
987 elseif (QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) .gt. 1.0) then
988 !#ifdef HYDRO_D
989 ! print *, "LatIn(Ql,Qsub,Qstrmvol)..",i,QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))), &
990 ! QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i))
991 !#endif
992 end if
994 elseif(LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .gt. 0 .and. &
995 ! (LAKE_MSKRT(CHANXI(i),CHANYJ(i)) .ne. -9999)) then !--a lake node
996 (CH_NETRT(CHANXI(i),CHANYJ(i)) .le. 0)) then !--a lake node
997 QLLAKE8(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) = &
998 QLLAKE8(LAKE_MSKRT(CHANXI(i),CHANYJ(i))) + &
999 (LAKEINFLORT(CHANXI(i),CHANYJ(i))+ &
1000 QINFLOWBASE(CHANXI(i),CHANYJ(i))) &
1001 /DT_STEPS*node_area(i)/1000/DTCT
1002 elseif(CH_NETRT(CHANXI(i),CHANYJ(i)) .gt. 0) then !pour out of lake
1003 QLateral(CH_NETLNK(CHANXI(i),CHANYJ(i))) = &
1004 QLAKEO(CH_NETRT(CHANXI(i),CHANYJ(i))) !-- previous timestep
1005 endif nodeType
1007 ENDDO
1010 #ifdef MPP_LAND
1011 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLateral,NLINKS,99)
1012 if(NLAKES .gt. 0) then
1013 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT ,ixrt,jxrt,QLLAKE,NLAKES,99)
1014 call sum_real8(QLLAKE8,NLAKES)
1015 QLLAKE = QLLAKE8
1016 endif
1017 #endif
1019 !-- compute conveyances, with known depths (just assign to QLINK(,1)
1020 !--QLINK(,2) will not be used), QLINK is the flow across the node face
1021 !-- units should be m3/second.. consistent with QL (lateral flow)
1023 #ifdef MPP_LAND
1024 DO iyw = 1,yw_MPP_NLINKS
1025 i = nlinks_index(iyw)
1026 #else
1027 DO i = 1,NLINKS
1028 #endif
1029 if (TYPEL(i) .eq. 0 .AND. HLINKTMP(FROM_NODE(i)) .gt. RETDEP_CHAN) then
1030 if(from_node(i) .ne. to_node(i) .and. (to_node(i) .gt. 0) .and.(from_node(i) .gt. 0) ) & ! added by Wei Yu
1031 QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
1032 HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
1033 CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
1034 else !-- we are just computing critical depth for outflow points
1035 QLINK(i,1) =0.
1036 endif
1037 ENDDO
1039 #ifdef MPP_LAND
1040 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
1041 #endif
1044 !-- compute total flow across face, into node
1045 #ifdef MPP_LAND
1046 DO iyw = 1,yw_mpp_nlinks
1047 i = nlinks_index(iyw)
1048 #else
1049 DO i = 1,NLINKS !-- inflow to node across each face
1050 #endif
1051 if(TYPEL(i) .eq. 0) then !-- only regular nodes have to attribute
1052 QSUM8(TO_NODE(i)) = QSUM8(TO_NODE(i)) + QLINK(i,1)
1053 endif
1054 END DO
1056 #ifdef MPP_LAND
1057 call MPP_CHANNEL_COM_REAL8(Link_location,ixrt,jxrt,qsum8,NLINKS,0)
1058 #endif
1059 qsum = qsum8
1063 #ifdef MPP_LAND
1064 do iyw = 1,yw_mpp_nlinks
1065 i = nlinks_index(iyw)
1066 #else
1067 do i = 1,NLINKS !-- outflow from node across each face
1068 #endif
1069 QSUM(FROM_NODE(i)) = QSUM(FROM_NODE(i)) - QLINK(i,1)
1070 end do
1071 #ifdef MPP_LAND
1072 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,qsum,NLINKS,99)
1073 #endif
1076 flag = 99
1079 #ifdef MPP_LAND
1080 do iyw = 1,yw_MPP_NLINKS
1081 i = nlinks_index(iyw)
1082 #else
1083 do i = 1, NLINKS !--- compute volume and depth at each node
1084 #endif
1086 if( TYPEL(i).eq.0 .and. CVOLTMP(i) .ge. 0.001 .and.(CVOLTMP(i)-QSUM(i)*DTCT)/CVOLTMP(i) .le. -0.01 ) then
1087 flag = -99
1088 #ifdef HYDRO_D
1089 write(6,*) "******* start diag ***************"
1090 write(6,*) "Unstable at node ",i, "i=",CHANXI(i),"j=",CHANYJ(i)
1091 write(6,*) "Unstatble at node ",i, "lat=",latval(CHANXI(i),CHANYJ(i)), "lon=",lonval(CHANXI(i),CHANYJ(i))
1092 write(6,*) "TYPEL, CVOLTMP, QSUM, QSUM*DTCT",TYPEL(i), CVOLTMP(i), QSUM(i), QSUM(i)*DTCT
1093 write(6,*) "qsubrt, qstrmvolrt,qlink",QSUBRT(CHANXI(i),CHANYJ(i)),QSTRMVOLRT(CHANXI(i),CHANYJ(i)),qlink(i,1),qlink(i,2)
1094 ! write(6,*) "current nodes, z, h", ZELEV(FROM_NODE(i)),HLINKTMP(FROM_NODE(i))
1095 ! if(TO_NODE(i) .gt. 0) then
1096 ! write(6,*) "to nodes, z, h", ZELEV(TO_NODE(i)), HLINKTMP(TO_NODE(i))
1097 ! else
1098 ! write(6,*) "no to nodes "
1099 ! endif
1100 write(6,*) "CHANLEN(i), MannN(i), Bw(i), ChSSlp(i) ", CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)
1101 write(6,*) "*******end of diag ***************"
1102 #endif
1104 goto 999
1105 endif
1106 enddo
1108 999 continue
1109 #ifdef MPP_LAND
1110 call mpp_same_int1(flag)
1111 #endif
1114 if(flag < 0 .and. DTCT >0.1) then
1116 ! call smoth121(HLINK,nlinks,maxv_p,pnode,to_node)
1118 if(DTCT .gt. minDTCT) then !-- timestep in seconds
1119 DTCT = max(DTCT/2 , minDTCT) !-- 1/2 timestep
1120 KRT = 0 !-- restart counter
1121 HLINKTMP = HLINK !-- set head and vol to start value of timestep
1122 CVOLTMP = CVOL
1123 CYCLE crnt !-- start cycle over with smaller timestep
1124 else
1125 write(6,*) "Courant error with smallest routing timestep DTCT: ",DTCT
1126 ! call hydro_stop("drive_CHANNEL")
1127 DTCT = 0.1
1128 HLINKTMP = HLINK !-- set head and volume to start values of timestep
1129 CVOLTMP = CVOL
1130 goto 998
1131 end if
1132 endif
1134 998 continue
1137 #ifdef MPP_LAND
1138 do iyw = 1,yw_MPP_NLINKS
1139 i = nlinks_index(iyw)
1140 #else
1141 do i = 1, NLINKS !--- compute volume and depth at each node
1142 #endif
1144 if(TYPEL(i) .eq. 0) then !-- regular channel grid point, compute volume
1145 CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) )* DTCT
1146 if((CVOLTMP(i) .lt. 0) .and. (gwChanCondSw == 0)) then
1147 #ifdef HYDRO_D
1148 print *, "WARNING! channel volume less than 0:i,CVOL,QSUM,QLat", &
1149 i, CVOLTMP(i),QSUM(i),QLateral(i),HLINK(i)
1150 #endif
1151 CVOLTMP(i) =0
1152 endif
1154 elseif(TYPEL(i) .eq. 1) then !-- pour point, critical depth downstream
1156 if (QSUM(i)+QLateral(i) .lt. 0) then
1157 else
1159 !DJG remove to have const. flux b.c.... CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
1160 CD(i) = HLINKTMP(i) !This is a temp hardwire for flow depth for the pour point...
1161 endif
1163 ! change in volume is inflow, lateral flow, and outflow
1164 !yw DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*DXRT),HLINKTMP(i), &
1165 CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
1166 DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
1167 CD(i),CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
1168 elseif (TYPEL(i) .eq. 2) then !--- into a reservoir, assume critical depth
1169 if ((QSUM(i)+QLateral(i) .lt. 0) .and. (gwChanCondSw == 0)) then
1170 #ifdef HYDRO_D
1171 print *, i, 'CrtDpth Qsum+QLat into lake< 0',QSUM(i), QLateral(i)
1172 #endif
1173 else
1174 !DJG remove to have const. flux b.c.... CD(i) =CRITICALDEPTH(i,abs(QSUM(i)+QLateral(i)), Bw(i), 1./ChSSlp(i))
1175 CD(i) = HLINKTMP(i) !This is a temp hardwire for flow depth for the pour point...
1176 endif
1178 !-- compute volume in reach (m^3)
1179 CVOLTMP(i) = CVOLTMP(i) + (QSUM(i) + QLateral(i) - &
1180 DIFFUSION(i,ZELEV(i),ZELEV(i)-(So(i)*CHANLEN(i)),HLINKTMP(i), &
1181 CD(i) ,CHANLEN(i), MannN(i), Bw(i), ChSSlp(i)) ) * DTCT
1182 else
1183 print *, "FATAL ERROR: This node does not have a type.. error TYPEL =", TYPEL(i)
1184 call hydro_stop("In drive_CHANNEL() - error TYPEL")
1185 endif
1187 if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
1188 HLINKTMP(i) = HEAD(i, CVOLTMP(i)/CHANLEN(i),Bw(i),1/ChSSlp(i)) !--updated depth
1189 else
1190 HLINKTMP(i) = CD(i) !!! CRITICALDEPTH(i,QSUM(i)+QLateral(i), Bw(i), 1./ChSSlp(i)) !--critical depth is head
1191 endif
1193 if(TO_NODE(i) .gt. 0) then
1194 if(LAKENODE(TO_NODE(i)) .gt. 0) then
1195 QLAKEI8(LAKENODE(TO_NODE(i))) = QLAKEI8(LAKENODE(TO_NODE(i))) + QLINK(i,1)
1196 endif
1197 endif
1199 END DO !--- done processing all the links
1201 #ifdef MPP_LAND
1202 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CVOLTMP,NLINKS,99)
1203 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,CD,NLINKS,99)
1204 if(NLAKES .gt. 0) then
1205 ! call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
1206 call sum_real8(QLAKEI8,NLAKES)
1207 QLAKEI = QLAKEI8
1208 endif
1209 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,HLINKTMP,NLINKS,99)
1210 #endif
1211 ! call check_channel(83,CVOLTMP,1,nlinks)
1212 ! call check_channel(84,CD,1,nlinks)
1213 ! call check_channel(85,HLINKTMP,1,nlinks)
1214 ! call check_lake(86,QLAKEI,lake_index,nlakes)
1216 do i = 1, NLAKES !-- mass balances of lakes
1217 #ifdef MPP_LAND
1218 if(lake_index(i) .gt. 0) then
1219 #endif
1221 ! Calls run for a single reservoir depending on its type as
1222 ! in whether it uses level pool, machine learning, or persistence.
1223 ! Inflow, lateral inflow, water elevation, and the routing period
1224 ! are passed in. Updated outflow and water elevation are returned.
1225 call rt_domain(did)%reservoirs(i)%ptr%run(QLAKEIP(i), QLAKEI(i), &
1226 QLLAKE(i), RESHT(i), QLAKEO(i), DTCT, rt_domain(did)%final_reservoir_type(i), &
1227 rt_domain(did)%reservoir_assimilated_value(i), rt_domain(did)%reservoir_assimilated_source_file(i))
1229 ! TODO: Encapsulate the lake state variable for water elevation (RESHT(i))
1230 ! inside the reservoir module, but it requires a redesign of the lake
1231 ! MPI communication model.
1233 QLAKEIP(i) = QLAKEI(i) !-- store total lake inflow for this timestep
1236 #ifdef MPP_LAND
1237 endif
1238 #endif
1239 enddo
1240 #ifdef MPP_LAND
1241 if(NLAKES .gt. 0) then
1242 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLLAKE,NLAKES,99)
1243 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,RESHT,NLAKES,99)
1244 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEO,NLAKES,99)
1245 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEI,NLAKES,99)
1246 !yw call MPP_CHANNEL_COM_REAL(LAKE_MSKRT,ixrt,jxrt,QLAKEIP,NLAKES,99)
1248 ! TODO: Change updateLake_grid calls below to updated distributed reservoir
1249 ! objects and not arrays as currently implemented.
1250 call updateLake_grid(QLLAKE, nlakes,lake_index)
1251 call updateLake_grid(RESHT, nlakes,lake_index)
1252 call updateLake_grid(QLAKEO, nlakes,lake_index)
1253 call updateLake_grid(QLAKEI, nlakes,lake_index)
1254 call updateLake_grid(QLAKEIP,nlakes,lake_index)
1255 endif
1256 #endif
1259 #ifdef MPP_LAND
1260 do iyw = 1,yw_MPP_NLINKS
1261 i = nlinks_index(iyw)
1262 #else
1263 do i = 1, NLINKS !--- compute volume and depth at each node
1264 #endif
1265 if(TYPEL(i) == 0) then !-- regular channel node, finalize head and flow
1266 QLINK(i,1)=DIFFUSION(i,ZELEV(FROM_NODE(i)),ZELEV(TO_NODE(i)), &
1267 HLINKTMP(FROM_NODE(i)),HLINKTMP(TO_NODE(i)), &
1268 CHANLEN(i), MannN(i), Bw(i), ChSSlp(i))
1269 endif
1270 enddo
1272 #ifdef MPP_LAND
1273 call MPP_CHANNEL_COM_REAL(Link_location,ixrt,jxrt,QLINK(:,1),NLINKS,99)
1274 #endif
1276 KRT = KRT + 1 !-- iterate on the timestep
1277 if(KRT .eq. DT_STEPS) exit crnt !-- up to the maximum time in interval
1279 end do crnt !--- DTCT timestep of DT_STEPS
1281 HLINK = HLINKTMP !-- update head based on final solution in timestep
1282 CVOL = CVOLTMP !-- update volume
1283 else !-- no channel option apparently selected
1284 print *, "FATAL ERROR: no channel option selected"
1285 call hydro_stop("In drive_CHANNEL() - no channel option selected")
1286 endif
1288 #ifdef HYDRO_D
1289 write(6,*) "finished call drive_CHANNEL"
1290 #endif
1292 if (KT .eq. 1) KT = KT + 1
1295 end subroutine drive_CHANNEL
1296 ! ----------------------------------------------------------------
1298 !-=======================================
1299 REAL FUNCTION AREAf(AREA,Bw,h,z)
1300 REAL :: AREA, Bw, z, h
1301 AREAf = (Bw+z*h)*h-AREA !-- Flow area
1302 END FUNCTION AREAf
1304 !-====critical depth function ==========
1305 REAL FUNCTION CDf(Q,Bw,h,z)
1306 REAL :: Q, Bw, z, h
1307 if(h .le. 0) then
1308 print *, "FATAL ERROR: head is zero, will get division by zero error"
1309 call hydro_stop("In CDf() - head is zero")
1310 else
1311 CDf = (Q/((Bw+z*h)*h))/(sqrt(9.81*(((Bw+z*h)*h)/(Bw + 2.0*z*h)))) - 1.0 !--critical depth function
1312 endif
1313 END FUNCTION CDf
1315 !=======find flow depth in channel with bisection Chapra pg. 131
1316 REAL FUNCTION HEAD(idx,AREA,Bw,z) !-- find the water elevation given wetted area,
1317 !--bottom widith and side channel.. index was for debuggin
1318 REAL :: Bw,z,AREA,test
1319 REAL :: hl, hu, hr, hrold
1320 REAL :: fl, fr,error !-- function evaluation
1321 INTEGER :: maxiter, idx
1323 error = 1.0
1324 maxiter = 0
1325 hl = 0.00001 !-- minimum depth is small
1326 hu = 30. !-- assume maximum depth is 30 meters
1328 if (AREA .lt. 0.00001) then
1329 hr = 0.
1330 else
1331 ! JLM: .gt. "0" is somewhat alarming. We really should have a constant like zero_r4
1332 do while ((AREAf(AREA,BW,hl,z)*AREAf(AREA,BW,hu,z)) .gt. 0.0 .and. maxiter .lt. 100)
1333 !-- allows for larger , smaller heads
1334 if(AREA .lt. 1.) then
1335 hl=hl/2
1336 else
1337 hu = hu * 2
1338 endif
1339 maxiter = maxiter + 1
1341 end do
1343 maxiter =0
1344 hr = 0
1345 fl = AREAf(AREA,Bw,hl,z)
1346 do while (error .gt. 0.0001 .and. maxiter < 1000)
1347 hrold = hr
1348 hr = (hl+hu)/2
1349 fr = AREAf(AREA,Bw,hr,z)
1350 maxiter = maxiter + 1
1351 if (hr .ne. 0) then
1352 error = abs((hr - hrold)/hr)
1353 endif
1354 test = fl * fr
1355 if (test.lt.0) then
1356 hu = hr
1357 elseif (test.gt.0) then
1358 hl=hr
1359 fl = fr
1360 else
1361 error = 0.0
1362 endif
1363 end do
1364 endif
1365 HEAD = hr
1367 22 format("i,hl,hu,Area",i5,2x,f12.8,2x,f6.3,2x,f6.3,2x,f6.3,2x,f9.1,2x,i5)
1369 END FUNCTION HEAD
1370 !=================================
1371 REAL FUNCTION MANNING(h1,n,Bw,Cs)
1373 REAL :: Bw,h1,Cs,n
1374 REAL :: z, AREA,R,Kd
1376 z = 1.0/Cs
1377 R = ((Bw + z*h1)*h1)/(Bw + 2.0*h1*sqrt(1.0 + z*z)) !-- Hyd Radius
1378 AREA = (Bw + z*h1)*h1 !-- Flow area
1379 Kd = (1.0/n) * (R**(2.0/3.0))*AREA !-- convenyance
1380 #ifdef HYDRO_D
1381 print *,"head, kd", h1,Kd
1382 #endif
1383 MANNING = Kd
1385 END FUNCTION MANNING
1387 !=======find flow depth in channel with bisection Chapra pg. 131
1388 REAL FUNCTION CRITICALDEPTH(lnk, Q, Bw, z) !-- find the critical depth
1389 REAL :: Bw, z, Q, test
1390 REAL :: hl, hu, hr, hrold
1391 REAL :: fl, fr, error !-- function evaluation
1392 INTEGER :: maxiter
1393 INTEGER :: lnk
1395 error = 1.0
1396 maxiter = 0
1397 hl = 1.0e-5 !-- minimum depth is 0.00001 meters
1398 ! hu = 35. !-- assume maximum critical depth 25 m
1399 hu = 100. !-- assume maximum critical depth 25 m
1401 if(CDf(Q, BW, hl, z) * CDf(Q, BW, hu, z) .gt. 0.0) then
1402 if(Q .gt. 0.001) then
1403 #ifdef HYDRO_D
1404 print *, "interval won't work to find CD of lnk ", lnk
1405 print *, "Q, hl, hu", Q, hl, hu
1406 print *, "cd lwr, upr", CDf(Q, BW, hl, z), CDf(Q, BW, hu, z)
1407 ! call hydro_stop("In CRITICALDEPTH()")
1408 CRITICALDEPTH = -9999
1409 return
1410 #endif
1411 else
1412 Q = 0.0
1413 endif
1414 endif
1416 hr = 0.
1417 fl = CDf(Q, Bw, hl, z)
1419 if (Q .eq. 0.0) then
1420 hr = 0.
1421 else
1422 do while (error .gt. 0.0001 .and. maxiter < 1000)
1423 hrold = hr
1424 hr = (hl+hu)/2.0
1425 fr = CDf(Q, Bw, hr, z)
1426 maxiter = maxiter + 1
1427 if (hr .ne. 0.0) then
1428 error = abs((hr - hrold)/hr)
1429 endif
1430 test = fl * fr
1431 if (test .lt. 0) then
1432 hu = hr
1433 elseif (test .gt. 0) then
1434 hl=hr
1435 fl = fr
1436 else
1437 error = 0.0
1438 endif
1440 end do
1441 endif
1443 CRITICALDEPTH = hr
1444 END FUNCTION CRITICALDEPTH
1447 REAL FUNCTION SGNf(val) !-- function to return the sign of a number
1448 REAL:: val
1450 if (val .lt. 0) then
1451 SGNf= -1.
1452 elseif (val.gt.0) then
1453 SGNf= 1.
1454 else
1455 SGNf= 0.
1456 endif
1458 END FUNCTION SGNf
1459 !================================================
1461 REAL FUNCTION fnDX(qp,Tw,So,Ck,dx,dt) !-- find channel sub-length for MK method
1462 REAL :: qp,Tw,So,Ck,dx, dt,test
1463 REAL :: dxl, dxu, dxr, dxrold
1464 REAL :: fl, fr, error
1465 REAL :: X
1466 INTEGER :: maxiter
1468 error = 1.0
1469 maxiter =0
1470 dxl = dx*0.9 !-- how to choose dxl???
1471 dxu = dx
1472 dxr=0
1474 do while (fnDXCDT(qp,Tw,So,Ck,dxl,dt)*fnDXCDT(qp,Tw,So,Ck,dxu,dt) .gt. 0 &
1475 .and. dxl .gt. 10) !-- don't let dxl get too small
1476 dxl = dxl/1.1
1477 end do
1480 fl = fnDXCDT(qp,Tw,So,Ck,dxl,dt)
1481 do while (error .gt. 0.0001 .and. maxiter < 1000)
1482 dxrold = dxr
1483 dxr = (dxl+dxu)/2
1484 fr = fnDXCDT(qp,Tw,So,Ck,dxr,dt)
1485 maxiter = maxiter + 1
1486 if (dxr .ne. 0) then
1487 error = abs((dxr - dxrold)/dxr)
1488 endif
1489 test = fl * fr
1490 if (test.lt.0) then
1491 dxu = dxr
1492 elseif (test.gt.0) then
1493 dxl=dxr
1494 fl = fr
1495 else
1496 error = 0.0
1497 endif
1498 end do
1499 FnDX = dxr
1501 END FUNCTION fnDX
1502 !================================================
1503 REAL FUNCTION fnDXCDT(qp,Tw,So,Ck,dx,dt) !-- function to help find sub-length for MK method
1504 REAL :: qp,Tw,So,Ck,dx,dt,X
1505 REAL :: c,b !-- coefficients on dx/cdt log approximation function
1507 c = 0.2407
1508 b = 1.16065
1509 X = 0.5-(qp/(2.0*Tw*So*Ck*dx))
1510 if (X .le. 0.0) then
1511 fnDXCDT = -1.0 !0.115
1512 else
1513 fnDXCDT = (dx/(Ck*dt)) - (c*LOG(X)+b) !-- this function needs to converge to 0
1514 endif
1515 END FUNCTION fnDXCDT
1516 ! ----------------------------------------------------------------------
1518 subroutine check_lake(unit,cd,lake_index,nlakes)
1519 use module_RT_data, only: rt_domain
1520 implicit none
1521 integer :: unit,nlakes,i,lake_index(nlakes)
1522 real cd(nlakes)
1523 #ifdef MPP_LAND
1524 call write_lake_real(cd,lake_index,nlakes)
1525 #endif
1526 write(unit,*) cd
1527 call flush(unit)
1528 return
1529 end subroutine check_lake
1531 subroutine check_channel(unit,cd,did,nlinks)
1532 use module_RT_data, only: rt_domain
1533 #ifdef MPP_LAND
1534 USE module_mpp_land
1535 #endif
1536 implicit none
1537 integer :: unit,nlinks,i, did
1538 real cd(nlinks)
1539 #ifdef MPP_LAND
1540 real g_cd(rt_domain(did)%gnlinks)
1541 call write_chanel_real(cd,rt_domain(did)%map_l2g,rt_domain(did)%gnlinks,nlinks,g_cd)
1542 if(my_id .eq. IO_id) then
1543 write(unit,*) "rt_domain(did)%gnlinks = ",rt_domain(did)%gnlinks
1544 write(unit,*) g_cd
1545 endif
1546 #else
1547 write(unit,*) cd
1548 #endif
1549 call flush(unit)
1550 close(unit)
1551 return
1552 end subroutine check_channel
1553 subroutine smoth121(var,nlinks,maxv_p,from_node,to_node)
1554 implicit none
1555 integer,intent(in) :: nlinks, maxv_p
1556 integer, intent(in), dimension(nlinks):: to_node
1557 integer, intent(in), dimension(nlinks):: from_node(nlinks,maxv_p)
1558 real, intent(inout), dimension(nlinks) :: var
1559 real, dimension(nlinks) :: vartmp
1560 integer :: i,j , k, from,to
1561 integer :: plen
1562 vartmp = 0
1563 do i = 1, nlinks
1564 to = to_node(i)
1565 plen = from_node(i,1)
1566 if(plen .gt. 1) then
1567 do k = 1, plen-1
1568 from = from_node(i,k+1)
1569 if(to .gt. 0) then
1570 vartmp(i) = vartmp(i)+0.25*(var(from)+2.*var(i)+var(to))
1571 else
1572 vartmp(i) = vartmp(i)+(2.*var(i)+var(from))/3.0
1573 endif
1574 end do
1575 vartmp(i) = vartmp(i) /(plen-1)
1576 else
1577 if(to .gt. 0) then
1578 vartmp(i) = vartmp(i)+(2.*var(i)+var(to)/3.0)
1579 else
1580 vartmp(i) = var(i)
1581 endif
1582 endif
1583 end do
1584 var = vartmp
1585 return
1586 end subroutine smoth121
1588 ! SUBROUTINE drive_CHANNEL for NHDPLUS
1589 ! ------------------------------------------------
1591 subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, &
1592 LAKEINFLORT, QSTRMVOLRT, TO_NODE, FROM_NODE, &
1593 TYPEL, ORDER, MAXORDER, CH_LNKRT, &
1594 LAKE_MSKRT, DT, DTCT, DTRT_CH,MUSK, MUSX, QLINK, &
1595 CHANLEN, MannN, So, ChSSlp, Bw, Tw, Tw_CC, n_CC, &
1596 ChannK, RESHT, CVOL, QLAKEI, QLAKEO, LAKENODE, &
1597 QINFLOWBASE, CHANXI, CHANYJ, channel_option, &
1598 nlinks,NLINKSL, LINKID, node_area, qout_gwsubbas, &
1599 LAKEIDA, LAKEIDM, NLAKES, LAKEIDX, &
1600 #ifdef MPP_LAND
1601 nlinks_index,mpp_nlinks,yw_mpp_nlinks, &
1602 LNLINKSL, &
1603 gtoNode,toNodeInd,nToNodeInd, &
1604 #endif
1605 CH_LNKRT_SL, landRunOff &
1606 #ifdef WRF_HYDRO_NUDGING
1607 , nudge &
1608 #endif
1609 , accSfcLatRunoff, accBucket &
1610 , qSfcLatRunoff, qBucket &
1611 , QLateral, velocity, qloss &
1612 , HLINK &
1613 , nsize , OVRTSWCRT, SUBRTSWCRT, channel_only, channelBucket_only, &
1614 channel_bypass)
1616 use module_UDMAP, only: LNUMRSL, LUDRSL
1617 use config_base, only: nlst
1619 #ifdef WRF_HYDRO_NUDGING
1620 use module_RT_data, only: rt_domain
1621 use module_stream_nudging, only: setup_stream_nudging, &
1622 nudge_term_all, &
1623 nudgeWAdvance, &
1624 nudge_apply_upstream_muskingumCunge
1625 #endif
1628 implicit none
1630 ! -------- DECLARATIONS ------------------------
1631 integer, intent(IN) :: did, IXRT,JXRT,channel_option, OVRTSWCRT, SUBRTSWCRT
1632 integer, intent(IN) :: NLAKES, NLINKSL, nlinks
1633 integer, intent(INOUT) :: KT ! flag of cold start (1) or continue run.
1634 real, intent(IN), dimension(IXRT,JXRT) :: QSTRMVOLRT
1635 real, intent(IN), dimension(IXRT,JXRT) :: LAKEINFLORT
1636 real, intent(IN), dimension(IXRT,JXRT) :: QINFLOWBASE
1637 real, dimension(ixrt,jxrt) :: landRunOff
1639 integer(kind=int64), intent(IN), dimension(IXRT,JXRT) :: CH_LNKRT
1640 integer(kind=int64), intent(IN), dimension(IXRT,JXRT) :: CH_LNKRT_SL
1642 integer, intent(IN), dimension(IXRT,JXRT) :: LAKE_MSKRT
1643 integer, intent(IN), dimension(:) :: ORDER, TYPEL !--link
1644 integer(kind=int64), intent(in), dimension(:) :: TO_NODE, FROM_NODE
1645 integer, intent(IN), dimension(:) :: CHANXI, CHANYJ
1646 real, intent(IN), dimension(:) :: MUSK, MUSX
1647 real, intent(IN), dimension(:) :: CHANLEN
1648 real, intent(IN), dimension(:) :: So, MannN
1649 real, intent(IN), dimension(:) :: ChSSlp,Bw,Tw !--properties of nodes or links
1650 real, intent(IN), dimension(:) :: Tw_CC, n_CC ! properties of compound channel
1651 real, intent(IN), dimension(:) :: ChannK !--channel infiltration
1652 real :: Km, X
1653 real , intent(INOUT), dimension(:,:) :: QLINK
1654 real , intent(INOUT), dimension(:) :: HLINK
1656 logical, intent(in) :: channel_bypass
1658 #ifdef WRF_HYDRO_NUDGING
1659 !! inout for applying previous nudge to upstream components of flow at gages
1660 real, intent(inout), dimension(:) :: nudge
1661 #endif
1663 real, dimension(:), intent(inout) :: QLateral !--lateral flow
1664 real, dimension(:), intent(out) :: velocity, qloss
1665 real*8, dimension(:), intent(inout) :: accSfcLatRunoff, accBucket
1666 real , dimension(:), intent(out) :: qSfcLatRunoff , qBucket
1668 real , dimension(NLINKSL,2) :: tmpQLINK
1669 real, intent(IN) :: DT !-- model timestep
1670 real, intent(IN) :: DTRT_CH !-- routing timestep
1671 real, intent(INOUT) :: DTCT
1672 real :: minDTCT !BF minimum routing timestep
1673 integer, intent(IN) :: MAXORDER
1674 real , intent(IN), dimension(:) :: node_area
1676 !DJG GW-chan coupling variables...
1677 real, dimension(NLINKS) :: dzGwChanHead
1678 real, dimension(NLINKS) :: Q_GW_CHAN_FLUX !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state...
1679 real, dimension(IXRT,JXRT) :: ZWATTBLRT !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state...
1681 !-- lake params
1682 integer(kind=int64), intent(IN), dimension(:) :: LAKEIDM !-- NHDPLUS lakeid for lakes to be modeled
1684 real, intent(INOUT), dimension(:) :: RESHT !-- reservoir height (m)
1685 real, intent(INOUT), dimension(:) :: QLAKEI !-- lake inflow (cms)
1686 real, dimension(NLAKES) :: QLAKEIP !-- lake inflow previous timestep (cms)
1687 real, intent(INOUT), dimension(NLAKES) :: QLAKEO !-- outflow from lake used in diffusion scheme
1689 integer(kind=int64), intent(IN), dimension(:) :: LAKENODE !-- outflow from lake used in diffusion scheme
1690 integer(kind=int64), intent(IN), dimension(:) :: LINKID !-- id of channel elements for linked scheme
1691 integer(kind=int64), intent(IN), dimension(:) :: LAKEIDA !-- (don't need) NHDPLUS lakeid for all lakes in domain
1692 integer, intent(IN), dimension(:) :: LAKEIDX !-- the sequential index of the lakes id by com id
1694 real, dimension(NLINKS) :: QSUM !--mass bal of node
1695 real, dimension(NLAKES) :: QLLAKE !-- lateral inflow to lake in diffusion scheme
1696 integer :: nsize
1698 !-- Local Variables
1699 integer :: i,j,k,t,m,jj,ii,lakeid, kk,KRT,node, UDMP_OPT
1700 integer :: DT_STEPS !-- number of timestep in routing
1701 real :: Qup,Quc !--Q upstream Previous, Q Upstream Current, downstream Previous
1702 real :: bo !--critical depth, bnd outflow just for testing
1704 real ,dimension(NLINKS) :: CD !-- critical depth
1705 real, dimension(IXRT,JXRT) :: tmp
1706 real, dimension(nlinks) :: tmp2
1707 real, intent(INOUT), dimension(:) :: CVOL
1709 #ifdef MPP_LAND
1710 real*8, dimension(LNLINKSL) :: LQLateral
1711 real*8, dimension(LNLINKSL) :: tmpLQLateral
1712 real, dimension(NLINKSL) :: tmpQLateral
1714 integer nlinks_index(:)
1715 integer iyw, yw_mpp_nlinks, mpp_nlinks
1716 real ywtmp(ixrt,jxrt)
1717 integer LNLINKSL
1718 integer, dimension(:) :: toNodeInd
1719 integer(kind=int64), dimension(:,:) :: gtoNode
1720 integer :: nToNodeInd
1721 real, dimension(nToNodeInd,2) :: gQLINK
1722 #else
1723 real*8, dimension(NLINKS) :: tmpLQLateral
1724 real, dimension(NLINKSL) :: tmpQLateral
1725 real, dimension(NLINKSL) :: LQLateral
1726 #endif
1727 integer flag
1728 integer, intent(in) :: channel_only, channelBucket_only
1730 integer :: n, kk2, nt, nsteps ! tmp
1731 real, intent(in), dimension(:) :: qout_gwsubbas
1732 real, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI, tmpRESHT
1733 integer, allocatable, dimension(:) :: tmpFinalResType
1734 real, allocatable,dimension(:) :: tmpAssimilatedValue
1735 character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile
1737 #ifdef MPP_LAND
1738 if(my_id .eq. io_id) then
1739 #endif
1740 allocate(tmpQLAKEO(NLAKES))
1741 allocate(tmpQLAKEI(NLAKES))
1742 allocate(tmpRESHT(NLAKES))
1743 allocate(tmpFinalResType(nlakes))
1744 #ifdef MPP_LAND
1745 endif
1746 #endif
1748 QLAKEIP = 0
1749 CD = 0
1750 node = 1
1751 QSUM = 0
1752 QLLAKE = 0
1753 dzGwChanHead = 0.
1754 nsteps = (DT+0.5)/DTRT_CH
1756 #ifdef WRF_HYDRO_NUDGING
1757 !! Initialize nudging for the current timestep.
1758 !! This establishes the data structure used to solve the nudges.
1759 call setup_stream_nudging(0) !! always zero b/c at beginning of hydro timestep
1760 #endif /* WRF_HYDRO_NUDGING */
1762 !---------------------------------------------
1763 if(channel_only .eq. 1 .or. channelBucket_only .eq. 1) then
1764 #ifdef HYDRO_D
1765 write(6,*) "channel_only or channelBucket_only is not zero. Special flux calculations."
1766 call flush(6)
1767 #endif /* HYDRO_D */
1769 ! if(nlst_rt(1)%output_channelBucket_influx .eq. 1 .or. &
1770 ! nlst_rt(1)%output_channelBucket_influx .eq. 2 ) &
1771 ! !! qScfLatRunoff = qLateral - qBucket
1772 ! qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL) - qout_gwsubbas(1:NLINKSL)
1774 if(nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1775 nlst(1)%output_channelBucket_influx .eq. 2 ) then
1777 if(channel_only .eq. 1) &
1778 !! qScfLatRunoff = qLateral - qBucket
1779 qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL) - qout_gwsubbas(1:NLINKSL)
1781 if(channelBucket_only .eq. 1) &
1782 !! qScfLatRunoff = qLateral - qBucket
1783 qSfcLatRunoff(1:NLINKSL) = qLateral(1:NLINKSL)
1785 end if
1787 if(nlst(1)%output_channelBucket_influx .eq. 3) &
1788 accSfcLatRunoff(1:NLINKSL) = qSfcLatRunoff * DT
1790 else
1792 QLateral = 0 !! the variable solved in this section. Channel only knows this already.
1793 LQLateral = 0 !-- initial lateral flow to 0 for this reach
1795 tmpQLateral = 0 !! WHY DOES THIS tmp variable EXIST?? Only for accumulations??
1796 tmpLQLateral = 0
1798 ! NHDPLUS maping
1799 if(OVRTSWCRT .eq. 0) then
1800 do k = 1, LNUMRSL
1801 ! get from land grid runoff
1802 do m = 1, LUDRSL(k)%ncell
1803 ii = LUDRSL(k)%cell_i(m)
1804 jj = LUDRSL(k)%cell_j(m)
1805 LQLateral(k) = LQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 &
1806 *LUDRSL(k)%cellArea(m)/DT
1807 tmpLQLateral(k) = tmpLQLateral(k)+landRunOff(ii,jj)*LUDRSL(k)%cellweight(m)/1000 &
1808 *LUDRSL(k)%cellArea(m)/DT
1809 end do
1810 end do
1812 #ifdef MPP_LAND
1813 call updateLinkV(tmpLQLateral, tmpQLateral)
1814 #endif
1815 if(NLINKSL .gt. 0) then
1816 if (nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1817 nlst(1)%output_channelBucket_influx .eq. 2 ) &
1818 qSfcLatRunoff(1:NLINKSL) = tmpQLateral(1:NLINKSL)
1819 if (nlst(1)%output_channelBucket_influx .eq. 3) &
1820 accSfcLatRunoff(1:NLINKSL) = accSfcLatRunoff(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
1821 endif
1822 tmpLQLateral = 0 !! JLM:: These lines imply that squeege runoff does not count towards
1823 tmpQLateral = 0 !! JLM:: accumulated runoff to be output but it does for internal QLateral?
1824 !! JLM: But then the next accumulation is added to the amt before zeroing? result
1825 !! JLM: should be identical to LQLateral.... I'm totally mystified.
1826 endif
1828 !! JLM:: if ovrtswcrt=0 and subrtswcrt=1, then this accumulation is calculated twice for LQLateral???
1829 !! This impiles that if overland routing is off and subsurface routing is on, that
1830 !! qstrmvolrt represents only the subsurface contribution to the channel.
1831 if(OVRTSWCRT .ne. 0 .or. SUBRTSWCRT .ne. 0 ) then
1832 do k = 1, LNUMRSL
1833 ! get from channel grid
1834 do m = 1, LUDRSL(k)%ngrids
1835 ii = LUDRSL(k)%grid_i(m)
1836 jj = LUDRSL(k)%grid_j(m)
1837 LQLateral(k) = LQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 &
1838 *LUDRSL(k)%nodeArea(m)/DT
1839 tmpLQLateral(k) = tmpLQLateral(k) + QSTRMVOLRT(ii,jj)*LUDRSL(k)%weight(m)/1000 &
1840 *LUDRSL(k)%nodeArea(m)/DT
1841 end do
1842 end do
1844 #ifdef MPP_LAND
1845 call updateLinkV(tmpLQLateral, tmpQLateral)
1846 #endif
1848 !! JLM:: again why output in this conditional ?? why not just output QLateral
1849 !! after this section ????
1850 if(NLINKSL .gt. 0) then
1851 if(nlst(1)%output_channelBucket_influx .eq. 1 .OR. &
1852 nlst(1)%output_channelBucket_influx .eq. 2 ) &
1853 qSfcLatRunoff(1:NLINKSL) = tmpQLateral(1:NLINKSL)
1854 if(nlst(1)%output_channelBucket_influx .eq. 3) &
1855 accSfcLatRunoff(1:NLINKSL) = accSfcLatRunoff(1:NLINKSL) + tmpQLateral(1:NLINKSL) * DT
1856 end if
1858 endif
1860 #ifdef MPP_LAND
1861 call updateLinkV(LQLateral, QLateral(1:NLINKSL))
1862 #else
1863 call hydro_stop("fatal error: NHDPlus only works for parallel now.")
1864 QLateral = LQLateral
1865 #endif
1866 endif !! (channel_only .eq. 1 .or. channelBucket_only .eq. 1) then; else; endif
1869 !---------------------------------------------
1870 !! If not running channelOnly, here is where the bucket model is picked up
1871 if(channel_only .eq. 1) then
1872 #ifdef HYDRO_D
1873 write(6,*) "channel_only is not zero. No bucket."
1874 call flush(6)
1875 #endif /* HYDRO_D */
1876 else
1877 !! REQUIRE BUCKET MODEL ON HERE?
1878 if(NLINKSL .gt. 0) QLateral(1:NLINKSL) = QLateral(1:NLINKSL) + qout_gwsubbas(1:NLINKSL)
1879 endif !! if(channel_only .eq. 1) then; else; endif
1881 if(nlst(1)%output_channelBucket_influx .eq. 1 .or. &
1882 nlst(1)%output_channelBucket_influx .eq. 2 ) &
1883 qBucket(1:NLINKSL) = qout_gwsubbas(1:NLINKSL)
1885 if(nlst(1)%output_channelBucket_influx .eq. 3) &
1886 accBucket(1:NLINKSL) = accBucket(1:NLINKSL) + qout_gwsubbas(1:NLINKSL) * DT
1888 ! Skip this section if we are NOT running any actual channel routing
1889 if (.not. channel_bypass) then
1891 !---------------------------------------------
1892 ! QLateral = QLateral / nsteps
1893 do nt = 1, nsteps
1895 #ifdef MPP_LAND
1896 gQLINK = 0
1897 call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2))
1898 call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1))
1899 !---------- route other reaches, with upstream inflow
1900 #endif
1902 tmpQlink = 0
1903 #ifdef MPP_LAND
1904 if(my_id .eq. io_id) then
1905 #endif
1906 tmpQLAKEO = QLAKEO
1907 tmpQLAKEI = QLAKEI
1908 tmpRESHT = RESHT
1909 tmpFinalResType = rt_domain(did)%final_reservoir_type
1910 tmpAssimilatedValue = rt_domain(did)%reservoir_assimilated_value
1911 tmpAssimilatedSourceFile = rt_domain(did)%reservoir_assimilated_source_file
1913 #ifdef MPP_LAND
1914 endif
1915 #endif
1917 do k = 1,NLINKSL
1919 Quc = 0
1920 Qup = 0
1922 !process as standard link or a lake inflow link, or lake outflow link
1923 ! link flowing out of lake, accumulate all the inflows with the revised TO_NODEs
1924 ! TYPEL = -999 stnd; TYPEL=1 outflow from lake; TYPEL = 3 inflow to a lake
1926 if(TYPEL(k) .ne. 2) then ! don't process internal lake links only
1928 #ifdef MPP_LAND
1929 !using mapping index
1930 do n = 1, gtoNODE(k,1)
1931 m = gtoNODE(k,n+1)
1932 !! JLM - I think gQLINK(,2) is actually previous. Global array never sees current. Seeing
1933 !! current would require global communication at the end of each loop through k
1934 !! (=kth reach). Additionally, how do you synchronize to make sure the upstream are all
1935 !! done before doing the downstream?
1936 if(gQLINK(m,2) .gt. 0) Quc = Quc + gQLINK(m,2) !--accum of upstream inflow of current timestep (2)
1937 if(gQLINK(m,1) .gt. 0) Qup = Qup + gQLINK(m,1) !--accum of upstream inflow of previous timestep (1)
1938 end do ! do i
1939 #else
1940 do m = 1, NLINKSL
1941 if (LINKID(k) .eq. TO_NODE(m)) then
1942 Quc = Quc + QLINK(m,2) !--accum of upstream inflow of current timestep (2)
1943 Qup = Qup + QLINK(m,1) !--accum of upstream inflow of previous timestep (1)
1944 endif
1945 end do ! do m
1946 #endif
1947 endif !note that we won't process type 2 links, since they are internal to a lake
1950 !yw ### process each link k,
1951 ! There is a situation that different k point to the same LAKEIDX
1952 ! if(TYPEL(k) .eq. 1 .and. LAKEIDX(k) .gt. 0) then !--link is a reservoir
1953 if(TYPEL(k) .eq. 1 ) then !--link is a reservoir
1955 lakeid = LAKEIDX(k)
1956 if(lakeid .ge. 0) then
1958 ! Calls run for a single reservoir depending on its type as
1959 ! in whether it uses level pool, machine learning, or persistence.
1960 ! Inflow, lateral inflow, water elevation, and the routing period
1961 ! are passed in. Updated outflow and water elevation are returned.
1962 call rt_domain(did)%reservoirs(lakeid)%ptr%run(Qup, Quc, 0.0, &
1963 RESHT(lakeid), tmpQLINK(k,2), DTRT_CH, rt_domain(did)%final_reservoir_type(lakeid), &
1964 rt_domain(did)%reservoir_assimilated_value(lakeid), rt_domain(did)%reservoir_assimilated_source_file(lakeid))
1966 ! TODO: Encapsulate the lake state variable for water elevation (RESHT(lakeid))
1967 ! inside the reservoir module, but it requires a redesign of the lake
1968 ! MPI communication model.
1970 QLAKEO(lakeid) = tmpQLINK(k,2) !save outflow to lake
1971 QLAKEI(lakeid) = Quc !save inflow to lake
1972 endif
1973 105 continue
1976 elseif (channel_option .eq. 1) then !muskingum routing
1977 Km = MUSK(k)
1978 X = MUSX(k)
1979 tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow
1981 elseif (channel_option .eq. 2) then ! muskingum cunge, don't process internal lake nodes TYP=2
1982 ! tmpQLINK(k,2) = MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), &
1983 ! QLateral(k), DTRT_CH, So(k), CHANLEN(k), &
1984 ! MannN(k), ChSSlp(k), Bw(k), Tw(k) )
1986 #ifdef WRF_HYDRO_NUDGING
1987 call nudge_apply_upstream_muskingumCunge( Qup, Quc, nudge(k), k )
1988 #endif
1990 call SUBMUSKINGCUNGE(&
1991 tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k), Qup, Quc, QLINK(k,1), &
1992 QLateral(k), DTRT_CH, So(k), CHANLEN(k), &
1993 MannN(k), ChSSlp(k), Bw(k), Tw(k), Tw_CC(k), n_CC(k), HLINK(k), ChannK(k) )
1995 else
1996 #ifdef HYDRO_D
1997 print *, " no channel option selected"
1998 #endif
1999 call hydro_stop("drive_CHANNEL")
2000 endif
2002 end do !--k links
2005 #ifdef MPP_LAND
2006 call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO)
2007 call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI)
2008 call updateLake_seq(RESHT,nlakes,tmpRESHT)
2009 call updateLake_seqInt(rt_domain(did)%final_reservoir_type, nlakes, tmpFinalResType)
2010 call updateLake_seq(rt_domain(did)%reservoir_assimilated_value, nlakes, tmpAssimilatedValue)
2011 !call updateLake_seq_char(rt_domain(did)%reservoir_assimilated_source_file, nlakes, tmpAssimilatedSourceFile)
2012 #endif
2014 do k = 1, NLINKSL !tmpQLINK?
2015 if(TYPEL(k) .ne. 2) then !only the internal lake nodes don't have info.. but need to save QLINK of lake out too
2016 QLINK(k,2) = tmpQLINK(k,2)
2017 endif
2018 QLINK(k,1) = QLINK(k,2) !assigng link flow of current to be previous for next time step
2019 end do
2022 #ifdef WRF_HYDRO_NUDGING
2023 if(.not. nudgeWAdvance) call nudge_term_all(qlink, nudge, int(nt*dtrt_ch))
2024 #endif /* WRF_HYDRO_NUDGING */
2027 !#ifdef HYDRO_D
2028 ! print *, "END OF ALL REACHES...",KRT,DT_STEPS
2029 !#endif
2031 end do ! nsteps
2033 endif ! channel_bypass
2035 if (KT .eq. 1) KT = KT + 1
2037 #ifdef MPP_LAND
2038 if(my_id .eq. io_id) then
2039 if(allocated(tmpQLAKEO)) deallocate(tmpQLAKEO)
2040 if(allocated(tmpQLAKEI)) deallocate(tmpQLAKEI)
2041 if(allocated(tmpRESHT)) deallocate(tmpRESHT)
2042 if(allocated(tmpFinalResType)) deallocate(tmpFinalResType)
2043 endif
2044 #endif
2046 if (KT .eq. 1) KT = KT + 1 ! redundant?
2048 end subroutine drive_CHANNEL_RSL
2050 ! ----------------------------------------------------------------
2052 end module module_channel_routing
2054 !! Is this outside the module scope on purpose?
2055 #ifdef MPP_LAND
2056 subroutine checkReach(ii, inVar)
2057 use module_mpp_land
2058 use module_RT_data, only: rt_domain
2059 use MODULE_mpp_ReachLS, only : updatelinkv, &
2060 ReachLS_write_io, gbcastvalue, &
2061 gbcastreal2
2062 implicit none
2063 integer :: ii
2064 real,dimension(rt_domain(1)%nlinksl) :: inVar
2065 real:: g_var(rt_domain(1)%gnlinksl)
2066 call ReachLS_write_io(inVar, g_var)
2067 if(my_id .eq. io_id) then
2068 write(ii,*) g_var
2069 call flush(ii)
2070 endif
2071 end subroutine checkReach
2072 #endif