subroutine imppro c c + + + PURPOSE + + + c c Channel water surface profile analysis subroutine. c SR IMPPRO uses the standard step method to obtain c the water surface profile and energy slopes at each c node in the channel. c c Called from: IMPINT c Author(s): A. Fogel c Reference in User Guide: c c Version: c Date recoded: c Recoded by: c c + + + KEYWORDS + + + c c + + + PARAMETERS + + + c c + + + ARGUMENT DECLARATIONS + + + c c + + + ARGUMENT DEFINITIONS + + + c c + + + COMMON BLOCKS + + + c include 'cimrou.inc' c read: b,ss,n,pypos,z,mdwit,dwt,rfrac,yc,frt,ncross c write: depth c c + + + LOCAL VARIABLES + + + c real area, hrad, cntrdp, v, vhead, theadc, y1, y2, hr23, r23, s1, 1 chr, dz, yy1, yy2, dy, ddx, stage, thead1, s2, sfavg, floss, 1 thdif, hhead, thead2, g, b, frt, ss, z(30), n, ydif integer iup, icsec, numnp1, nod, i, iter c c + + + LOCAL DEFINITIONS + + + c c area - the flow area of the channel c b - channel bottom width c chr - channel roughness coefficient c cntrdp - depth at the control section c ddx - distance downstream c dy - incremental change in depth c dz - incremental change in downstream distance c floss - friction loss c frt - flow rate c g - gravitational constant c hhead - head variable c hr23 - hydraulic radius ** 2/3 c hrad - hydraulic radius c icsec - number of cross sections c iter - number of iterations c iup - upstream routing flag c n - manning's n c nod - normal depth counter c numnp1 - number of points in cross section c r23 - 2/3 in manning's n equation c s1 - friction slope c s2 - friction slope c sfavg - friction slope c ss - side slopes c stage - stage of the flow depth c thdif - convergence head difference c thead1 - total head c thead2 - total head c theadc - total head c v - average flow velocity c vhead - velocity head c y1 - flow depth c y2 - flow depth c ydif - flow depth difference c yy1 - control section position c yy2 - control section position + 1 c z(i) - horizontal location c c + + + SAVES + + + c c + + + SUBROUTINES CALLED + + + c c + + + DATA INITIALIZATIONS + + + c c + + + END SPECIFICATIONS + + + c c calculate the hydraulic radius for the control section c b = bwes ss = sses n = nes c do 10 i = 1, ncross z(i) = zes(i) 10 continue c frt = qesr iup = 0 cntrdp = yc g = 32.2 c if (iup.eq.0) then icsec = ncross depth(icsec) = cntrdp area = b * depth(icsec) + ss * depth(icsec) ** 2 hrad = area / (b+2.*depth(icsec)*(ss**2+1)**0.5) else icsec = 1 depth(icsec) = cntrdp area = b * depth(icsec) + ss * depth(icsec) ** 2 hrad = area / (b+2.*depth(icsec)*(ss**2+1)**0.5) end if c c calculate the velocity head and the total head for the control c section c v = frt / area vhead = v * v / (2.0*g) theadc = pypos(icsec) + depth(icsec) + vhead y1 = depth(icsec) y2 = y1 c c calculate the friction slope using manning's equation c r23 = 2.0 / 3.0 chr = n / 1.486 hr23 = hrad ** r23 s1 = (v*chr/hr23) ** 2.0 c c begin loop over the cross sections and set up point to change depth c toward normal depth c numnp1 = ncross - 1 c do 30 nod = 1, numnp1 c if (iup.eq.0) then icsec = ncross - nod else icsec = nod + 1 end if c c calculate the distance downstream (or upstream) c if (iup.eq.0) then dz = z(icsec+1) - z(icsec) yy1 = pypos(icsec) yy2 = pypos(icsec+1) dy = yy2 - yy1 ddx = sqrt(dy*dy+dz*dz) else dz = z(icsec) - z(icsec-1) yy1 = pypos(icsec-1) yy2 = pypos(icsec) dy = yy2 - yy1 ddx = sqrt(dy*dy+dz*dz) end if c c set first guess at the stage to that of the previous cross section c and set the depth equal to that of the previous cross section c iter = 0 c 20 iter = iter + 1 c if (iter.gt.mdwit) then write (6,1000) iter, mdwit, icsec stop end if c stage = pypos(icsec) + y2 y1 = y2 c c calculate the cross sectional area, hydraulic radius, and velocity c and total head at the cross section c area = b * y1 + ss * y1 ** 2. hrad = area / (b+2.*y1*(ss**2.+1)**0.5) c v = frt / area vhead = v * v / (2.0*g) thead1 = stage + vhead c c calculate the friction slope using manning's equation c hr23 = hrad ** r23 s2 = (v*chr/hr23) ** 2.0 c c calculate the average friction slope over the reach and c the friction loss c sfavg = 0.5 * (s1+s2) floss = sfavg * ddx c c calculate the new total head c if (iup.eq.0) then thead2 = theadc + floss else thead2 = theadc - floss end if c c check for convergence c thdif = thead2 - thead1 c if (abs(thdif).lt.dwt) then theadc = thead2 s1 = s2 depth(icsec) = thead2 - vhead - pypos(icsec) y2 = yc go to 30 c else c hhead = 0.5 * (thead1+thead2) c c y2 = y1*(thead2/thead1) c y2 = hhead - pypos(icsec) - vhead ydif = y2 - y1 c if (abs(ydif).gt.rfrac*y1) then if (ydif.ge.0.0) then y2 = y1 + rfrac * y1 else y2 = y1 - rfrac * y1 end if end if c go to 20 end if c 30 continue c return c 1000 format (//12x,'!!!!!!!!!!!!!! E R R O R !!!!!!!!!!!!!!',// 1 7x,'WATER SURFACE PROFILE ITERATIONS (',i3, 1 ') EXCEEDED MAXIMUM (',i3,')',/28x,'PROGRAM HALTED'/19x, 1 'AT CROSS SECTION ',i2//12x, 1 '!!!!!!!!!!!!!! E R R O R !!!!!!!!!!!!!!') end