subroutine impflo(itr,qi,htw) c c + + + PURPOSE + + + c c SR IMPFLO is an impoundment element subroutine responsible c for determining the change in stage and outflow over a c time step. It also calculates the duration of the next time c step. c c Called from: IMPMAI c Author(s): Mark Lindley c Reference in User Guide: c c Version: c Date recoded: 3/01/94 c Recoded by: Jim Ascough II c c + + + KEYWORDS + + + c c + + + PARAMETERS + + + c include 'pmximp.inc' c c + + + ARGUMENT DECLARATIONS + + + c real qi, htw integer itr c c + + + ARGUMENT DEFINITIONS + + + c c itr - outflow regime number c htw - depth of tailwater flow in downstream channel c qi - inflow rate c c + + + COMMON BLOCKS + + + c include 'cimitf.inc' c read: a,b,c,d,e,ha,ht,hlm,a0,a1,a2,l0,l1,l2,qinf,isize c include 'cimpnd.inc' c include 'cimqot.inc' c read: qo,h,t,dttry,deltat c write: qon,hn,tn,dtdid,dtnext c c + + + LOCAL VARIABLES + + + c real qo1, qo2, qo3, qo4, qo5, qo6, qo7, qo8, qo9, qo10, qo11, 1 qo12, qo13, qo14, qo15, tsav, hsav, dt, tt, htemp, errmax integer flg, k c c + + + LOCAL DEFINITIONS + + + c c dt - intermediate time step c errmax - error in making one step compared to two half steps c flg - flag for an attempted step that is too large c hsav - intermediate stage (at the beginning of the time step) c htemp - intermediate stage (at the end of the time step) c k - counter for number of attempts made c qo1,15 - outflow contribution for each stage-discharge function c tsav - intermediate time c tt - intermediate time step c c + + + SAVES + + + c c + + + SUBROUTINES CALLED + + + c c imphnw c c + + + DATA INITIALIZATIONS + + + c c + + + END SPECIFICATIONS + + + c c c initialize variables c k = 0 c tsav = t(ipond) hsav = h(ipond) dt = dttry(ipond) c 10 tt = 0.5 * dt c flg = 0 c c attempting integration with two consecutive half time steps c call imphnw(hsav,flg,tt,htemp,qi,itr,htw) c if (flg.eq.1) then dt = 0.5 * dt go to 10 end if c t(ipond) = tsav + tt c call imphnw(htemp,flg,tt,h(ipond),qi,itr,htw) c if (flg.eq.1) then dt = 0.5 * dt go to 10 end if c t(ipond) = tsav + dt c c attempting integration with one full time step c call imphnw(hsav,flg,dt,htemp,qi,itr,htw) c if (flg.eq.1) then dt = 0.5 * dt go to 10 end if c c checking error between two half steps and one full time step c if error is too large the time step is decreased and c attempted again, if the error is acceptable, the next c suggested time step is computed c htemp = h(ipond) - htemp errmax = max(0.0,abs(htemp)) / 1.0e-04 c if (errmax.gt.1) then dt = 0.9 * dt * errmax ** (-0.25) go to 10 else dtdid(ipond) = dt c if (errmax.gt.6.0e-04) then dtnext(ipond) = 0.9 * dt * (errmax**(-0.2)) else dtnext(ipond) = 4.0 * dt end if c end if c hn(ipond) = h(ipond) + htemp / 15.0 c c checking to see that the new stage is within the same flow c regime, if it is not the time step is attempted again c using the initial time step c k = k + 1 c if (k.ne.2) then c if (hn(ipond).gt.ht(itr+1,ipond)) then dt = deltat(ipond) go to 10 else c if (itr.gt.1) then c if (hn(ipond).lt.ht(itr,ipond)) then dt = deltat(ipond) go to 10 end if c end if end if end if c h(ipond) = hsav tn(ipond) = tsav + dtdid(ipond) t(ipond) = tsav c c computing the new outflow c if (hn(ipond).le.hlm(itr,ipond)) hn(ipond) = hlm(itr,ipond) c qo1 = (b(1,itr,ipond)*(hn(ipond)-ha(1,itr,ipond))** 1 c(1,itr,ipond)) c qo2 = (b(2,itr,ipond)*(hn(ipond)-ha(2,itr,ipond))** 1 c(2,itr,ipond)) c if (htw.gt.a(3,itr,ipond)) then qo3 = (b(3,itr,ipond)*(hn(ipond)-(ha(3,itr,ipond)+htw- 1 a(3,itr,ipond)))**c(3,itr,ipond)) else qo3 = (b(3,itr,ipond)*(hn(ipond)-ha(3,itr,ipond))** 1 c(3,itr,ipond)) end if c qo4 = a(4,itr,ipond) * ((hn(ipond)-ha(4,itr,ipond))/ 1 b(4,itr,ipond)) ** c(4,itr,ipond) c qo5 = a(5,itr,ipond) * (((hn(ipond)-ha(5,itr,ipond))/ 1 b(5,itr,ipond)+c(5,itr,ipond))/d(5,itr,ipond)) ** 0.5 c if (htw.gt.a(6,itr,ipond)) then qo6 = b(6,itr,ipond) * (hn(ipond)-(ha(6,itr,ipond)+htw- 1 a(6,itr,ipond))) ** c(6,itr,ipond) else qo6 = b(6,itr,ipond) * (hn(ipond)-ha(6,itr,ipond)) ** 1 c(6,itr,ipond) end if c qo7 = a(7,itr,ipond) * ((hn(ipond)-ha(7,itr,ipond))/ 1 b(7,itr,ipond)) ** c(7,itr,ipond) c qo8 = a(8,itr,ipond) * (((hn(ipond)-ha(8,itr,ipond))/ 1 b(8,itr,ipond)+c(8,itr,ipond))/d(8,itr,ipond)) ** 0.5 c if (htw.gt.a(9,itr,ipond)) then qo9 = b(9,itr,ipond) * (hn(ipond)-(ha(9,itr,ipond)+htw- 1 a(9,itr,ipond))) ** c(9,itr,ipond) else qo9 = b(9,itr,ipond) * (hn(ipond)-ha(9,itr,ipond)) ** 1 c(9,itr,ipond) end if c qo10 = a(10,itr,ipond) * ((hn(ipond)-ha(10,itr,ipond))/ 1 b(10,itr,ipond)) ** c(10,itr,ipond) + d(10,itr,ipond) * ( 1 hn(ipond)-e(10,itr,ipond)) ** 1.5 c qo11 = a(11,itr,ipond) + b(11,itr,ipond) * (hn(ipond)- 1 ha(11,itr,ipond)) + c(11,itr,ipond) * (hn(ipond)- 1 ha(11,itr,ipond)) ** 2. + d(11,itr,ipond) * (hn(ipond)- 1 ha(11,itr,ipond)) ** 3. + e(11,itr,ipond) * (hn(ipond)- 1 ha(11,itr,ipond)) ** 4.0 c qo12 = a(12,itr,ipond) * (hn(ipond)-ha(12,itr,ipond)) + ( 1 b(12,itr,ipond)+c(12,itr,ipond)*(hn(ipond)-d(12,itr,ipond))) * 1 (hn(ipond)-d(12,itr,ipond)) ** 1.5 c qo13 = a(13,itr,ipond) / (b(13,itr,ipond)+c(13,itr,ipond)/( 1 hn(ipond)-ha(13,itr,ipond))**1.5) c qo14 = a(14,itr,ipond) * (hn(ipond)-ha(14,itr,ipond)) ** 0.5 c qo15 = b(15,itr,ipond) * (hn(ipond)-ha(15,itr,ipond)) ** 1 c(15,itr,ipond) c qon(ipond) = min(qo1,qo2,qo3) + min(qo4,qo5,qo6) + 1 min(qo7,qo8,qo9) + qo10 + qo11 + qo12 + min(qo13,qo14,qo15) c return end