SUBROUTINE SWMM5_NR_ITERATIVE
C EXTRAN BLOCK test for SWMM 5 beta solution - 12/12/2002
INCLUDE 'TAPES.INC'
INCLUDE 'STIMER.INC'
INCLUDE 'BD.INC'
INCLUDE 'BND.INC'
INCLUDE 'HYFLOW.INC'
INCLUDE 'CONTR.INC'
INCLUDE 'JUNC.INC'
INCLUDE 'PIPE.INC'
INCLUDE 'TIDE.INC'
INCLUDE 'OUT.INC'
INCLUDE 'ORF.INC'
INCLUDE 'WEIR.INC'
INCLUDE 'FLODAT.INC'
DOUBLE PRECISION AKON,QNEW,DELQ1,DELQ2,DELQ3,DELQ4,DELQ5
DIMENSION AS1(NEE)
DOUBLE PRECISION df,f,n_omega
integer good_nodes
C=======================================================================
C STORE OLD TIME STEP FLOW VALUES
C=======================================================================
DO N = 1,NTL
QO(N) = Q(N)
AT(N) = A(N)
VT(N) = V(N)
enddo
C=======================================================================
C INITIALIZE CONTINUITY PARAMETERS
C=======================================================================
DO J = 1,NJ
if(othercom(63).eq.1) then
Y(J) = Y(J) + 0.5 * ( YEX2(J) -
+ YEX1(J) + YEX1(J) - YO(J))
IF(Y(J).LT.FUDGE) Y(J)=FUDGE
IF(Y(J).GT.SURELEV(J)) Y(J)=SURELEV(J)-Z(J)
yex2(j) = yex1(j)
yex1(j) = YO(j)
endif
cred beginning time step value of the node area
asold(j) = as(j)
YO(J) = Y(J)
GoodNode(j) = .FALSE.
enddo
good_nodes = 0
omega = input_omega
n_omega = node_omega
loop_count = 0
C=======================================================================
C HALF-STEP AREA, RADIUS : VELOCITY
C FULL-STEP FLOW
C=======================================================================
big_loop: DO while(good_nodes.lt.nj.and.loop_count.le.itmax)
loop_count = loop_count + 1
if(loop_count.ge.itmax-5) then
omega = 0.50 * omega
n_omega = 0.50 * n_omega
endif
DO J = 1,NJ
AS(J) = AMEN
AS1(J) = 0.0
SUMQ(J) = QIN(J)
SUMQS(J) = QIN(J)
SUMAL(J) = 0.0
enddo
CIM FIRST COMPUTE GATED ORIFICE PARAMETERS
CALL OGATES(DELT,Y,V)
c
flow_loop: DO N = 1,NTC
NL = NJUNC(N,1)
NH = NJUNC(N,2)
C=======================================================================
H(N,1) = AMAX1(Y(NL) + Z(NL),ZU(N))
H(N,2) = AMAX1(Y(NH) + Z(NH),ZD(N))
CALL nhead(N,NL,NH,H(N,1),H(N,2),Q(N),A(N),V(N),HRAD,
+ ANH,ANL,RNL,RNH,YNL,YNH,width,IDOIT,LINK(N),AS1)
cred bypass loop for nodes already converged
bypass_loop: if(loop_count.gt.2.and.
+ goodnode(nl).eq..TRUE..and.goodnode(nh).eq..TRUE.) then
bypass = bypass + 1.0
else
IF(HRAD.GT.HMAX(N)) HMAX(N) = HRAD
IF(A(N).GT.AMAX(N)) AMAX(N) = A(N)
cred save information for the culvert classification
HRLAST(N) = HRAD
vup(n) = anl
vdn(n) = anh
if(loop_count.eq.1) then
aup(n) = anl
rup(n) = rnl
rmd(n) = hrad
endif
positive_flow: IF(IDOIT.gt.0) THEN
c
c Q/ANH = velocity at downstream end used for exit loss
c Q/ANL = velocity at upstream end used for entrance loss
c The loss = 1/2*K*V^2/g is factored into momentum
c equation similarly to the friction slope
c The loss momentum term = g*A*[1/2*K*V^2/g]
c = g*A*[1/2*K*Q/A*Q/A/g]
C =g *[1/2*K*|Q*A| /g] * Q
C = [1/2*K*|Q*A| ] * Q
c DELQ5 is sum of losses
c
cred calculate the froude number for every conduit
c
DELH = H(N,1) - H(N,2)
DELZP = ZU(N) - ZD(N)
DIFF_mid = 0.5 * (H(N,1) - ZU(N) + H(N,2) - ZD(N))
IF(DIFF_mid.GT.0.0) THEN
FROUDE_MID = ABS(Q(N))/A(N)/SQRT(GRVT*(DIFF_mid))
ELSE
FROUDE_MID = 0.0
ENDIF
bfactor = 0.0
if(FROUDE_MID.ge.1.0) THEN
bfactor = 0.0
elseif(froude_mid.gt.0.9) then
bfactor = 1.0 - (10.0*FROUDE_MID-9.0)
endif
cred test for zero sloped conduits
if(delzp.eq.0.0) then
del_ratio = delh/0.001
else
del_ratio = delh/delzp
endif
cred swmm 4 definition for normal flow
if(del_ratio.le.1.0) then
delfactor = 0.0
cred swmm 5 transition definition
else if(del_ratio.gt.1.10) then
delfactor = 1.0
else
delfactor = 10.0 * (del_ratio-1.0)
endif
DELQ5 = 0.0
IF(ANH.NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/ANH)*ENTK(N)
IF(ANL.NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/ANL)*EXITK(N)
IF(A(N).NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/A(N))*OTHERK(N)
DELQ5 = DELQ5 * DELT/LEN(N)
c
DELQ4 = DELT*V(N)*(ANH-ANL)/A(N)/LEN(N)
cred the mean area travels from the midpoint to the upstream area
cred as the value of delh/delzp changes
area_mean = wt*anl + wd*aup(n) + ( wt*a(n) + wd*at(n) -
+ wt*anl - wd*aup(n)) * delfactor
DELQ2 = DELT*GRVT*area_mean*((H(N,2) - H(N,1)))/LEN(N)
DELQ3 = 2.0*(A(N)-AT(N))/A(N)
cred the mean hydraulic radius travels from the midpoint to the
cred upstream hydraulic radius as the value of delh/delzp changes
hrad_mean = wt*rnl + wd*rup(n) + ( wt*hrad + wd*rmd(n) -
+ wt*rnl - wd*rup(n)) * delfactor
DELQ1 = DELT*(ROUGH(N)/hrad_mean**1.33333)*ABS(V(N))
QNEW = QO(N) - delq2
AKON = DELQ1 + DELQ5 - delq4*bfactor - DELQ3*bfactor
cred Newton-Raphson iteration
F = Q(N) * ( 1.0 + akon) - qnew
DF = 1.0 + akon
Q(N) = Q(N) - omega*f/df
DQDH = 1.0/(1.0+AKON)*GRVT*DELT*A(N)/LEN2(N)
dqdh_old(n) = dqdh
C=======================================================================
C DO NOT ALLOW A FLOW REVERSAL IN ONE TIME STEP
C=======================================================================
DIRQT = SIGN(1.0,QO(N))
DIRQ = SIGN(1.0,Q(N))
IF(DIRQT/DIRQ.LT.0.0) Q(N) = 0.001*DIRQ
v(n) = q(n)/A(N)
C=======================================================================
C COMPUTE CONTINUITY PARAMETERS
C=======================================================================
SELECT CASE (INGATE(N))
CASE (1)
Q(N) = AMAX1(Q(N),0.0)
CASE (2)
Q(N) = AMIN1(Q(N),0.0)
END SELECT
IF (NKLASS(N).LE.21) THEN
Q(N) = AMAX1(STHETA(N),Q(N))
Q(N) = AMIN1(SPHI(N),Q(N))
ENDIF
endif positive_flow
endif bypass_loop
SUMQS(NL) = SUMQS(NL) - Q(N)*barrels(n)
SUMAL(NL) = SUMAL(NL) + dqdh_old(n)*barrels(n)
SUMQS(NH) = SUMQS(NH) + Q(N)*barrels(n)
SUMAL(NH) = SUMAL(NH) + dqdh_old(n)*barrels(n)
SUMQ(NL) = SUMQ(NL) - WT*Q(N)*barrels(n) - WD*QO(N)*barrels(n)
SUMQ(NH) = SUMQ(NH) + WT*Q(N)*barrels(n) + WD*QO(N)*barrels(n)
ENDDO flow_loop
C=======================================================================
C SET FULL STEP OUTFLOWS AND INTERNAL TRANSFERS
C=======================================================================
CALL BOUND(Y,Y,Q,TIME,DELT)
C=======================================================================
N1 = NTC+1
DO 370 N = N1,NTL
NL = NJUNC(N,1)
NH = NJUNC(N,2)
IF(ABS(Q(N)).LT.1E-10) Q(N) = 0.0
C=======================================================================
SUMQ(NL) = SUMQ(NL) - WT*Q(N)*barrels(n) - WD*QO(N)*barrels(n)
SUMQS(NL) = SUMQS(NL) - Q(N)*barrels(n)
IF(NH.NE.0) THEN
SUMQ(NH) = SUMQ(NH) + WT*Q(N)*barrels(n) + WD*QO(N)*barrels(n)
SUMQS(NH) = SUMQS(NH) + Q(N)*barrels(n)
ENDIF
370 CONTINUE
C=======================================================================
C CALCULATE THE FULL-STEP HEAD
C=======================================================================
DO J = 1, NJ
AS(J) = AS(J) + AS1(J)
ENDDO
DO J = 1,NJ
IF(JSKIP(J).le.0) then
C=======================================================================
cred time weighted area average
if(othercom(62).eq.1) then
average_area = WT*as(j) + WD*asold(j)
else
average_area = as(j)
endif
cred Newton-Raphson iteration
yt(j) = y(j)
if(y(j).le.ycrown(j)) then
F = Y(J) * average_area - YO(J) * average_area
+ - sumq(j)*DELT
DF = average_area
ASFULL(J) = AS(J)
else
DF = sumal(j)
IF(Y(J).LT.1.25*YCROWN(J)) DF = DF +
+ (ASFULL(J)/DELT-SUMAL(J))*EXP(-15.*(Y(J)-YCROWN(J))/YCROWN(J))
cred correction from SWMM 5 QA testing - 11/12/2004
cred if a large delt or if there is a large value of sumal(j)
cred - usually from a very large conduit connected to the current node -
cred the expression ASFULL(J)/DELT-SUMAL(J) may be negative
cred invalidating the whole concept of a "transition" slot
if(Y(J).LT.1.25*YCROWN(J).le.ASFULL(J)/DELT.le.SUMAL(J))
+ denom = sumal(j)
CORR = 1.00
IF(NCHAN(J,2).EQ.0) CORR = 0.60
F = Y(J) * DF - YO(J) * DF - CORR * sumqs(j)
cred WRITE(*,*) ajun(j),df,f,y(j),ycrown(j),sumal(j)
endif
Y(J) = Y(J) - n_omega*F/DF
IF(Y(J).LT.0.0) Y(J) = 0.0
IF((Y(J)+Z(J)).GT.SURELEV(J)) Y(J) = SURELEV(J)-Z(J)
endif
enddo
cred converge until all nodes meet the convergence criteria
good_nodes = 0
error_node = 0
i_was_bad = 1
DO j = 1,nj
if(jskip(j).le.0) then
if(abs(y(j)-yt(j)).le.node_toler) then
good_nodes = good_nodes + 1
GoodNode(j) = .TRUE.
endif
if(abs(y(j)-yt(j)).gt.error_node) then
error_node = abs(y(j)-yt(j))
i_was_bad = j
endif
else
good_nodes = good_nodes + 1
GoodNode(j) = .TRUE.
endif
enddo
c write(*,669) loop_count,nj-good_nodes,error_node,
c + delt,ajun(i_was_bad),y(i_was_bad),
c + sumq(i_was_bad),omega
669 format(2I6,F10.4,F8.2,1x,a12,3F9.3)
enddo big_loop
cred culvert information - 8/6/2002
cred culvert information for the culvert comparison file - 8/6/2002
cred culvert information - 8/6/2002
if(othercom(90).eq.1) then
do n = 1,nc
if(nklass(n).eq.1.or.nklass(n).eq.21) then
if(abs(q(n)).gt.culvert_loss(n,5).and.abs(v(n)).lt.20.0)then
HW = H(N,1) - ZU(N)
TW = H(N,2) - ZD(N)
vup(n) = v(n)
vdn(n) = v(n)
head_loss = 0.5*ENTK(N)*vup(N)*vup(N)/GRVT +
+ 0.5*EXITK(N)*vdn(N)*vdn(N)/GRVT +
+ 0.5*OTHERK(N)*v(N)*v(N)/GRVT
sfloss = ROUGH(N)/GRVT *
+ abs(v(N))*ABS(v(N))/hrlast(n)**1.33333
cred sfloss = ROUGH(N)/GRVT *
cred + Q(N)*ABS(Q(N))/(A(N)**2*HRLAST(N)**1.33333)
culvert_loss(n,1) = len(n)*sfloss
culvert_loss(n,2) = head_loss
culvert_loss(n,3) = h(n,1)
culvert_loss(n,4) = h(n,2)
culvert_loss(n,5) = q(n)
endif
CALL DEPTHX(N,NKLASS(N),q(n),YC,YNORM)
cred type 1
IF(HW.LT.1.5*DEEP(N).AND.YC.LT.YNORM.and.
+ TW.LE.YC) then
culverted(N,1) = culverted(N,1) + DELT
culverted(n,8) = 1
endif
cred type 2
IF(HW.LT.1.5*DEEP(N).AND.YC.LT.YNORM.AND.
+ TW.GT.YC.AND.TW.LE.DEEP(N)) then
culverted(N,2) = culverted(N,2) + DELT
culverted(n,8) = 2
endif
cred type 3
IF(YC.GE.YNORM.AND.TW.LT.ZU(N)-ZD(N)) then
culverted(N,3) = culverted(N,3) + DELT
culverted(n,8) = 3
endif
cred type 4
IF(YC.GE.YNORM.AND.TW.GE.ZU(N)-ZD(N)+YC.
+ AND.TW.LT.ZU(N)-ZD(N)+DEEP(N)) then
culverted(N,4) = culverted(N,4) + DELT
culverted(n,8) = 4
endif
cred type 5
IF(YC.LT.YNORM.AND.TW.GE.DEEP(N)) then
culverted(N,5) = culverted(N,5) + DELT
culverted(n,8) = 5
endif
cred type 6
IF(YC.GE.YNORM.AND.TW.GE.ZU(N)-ZD(N)+
+ DEEP(N)) then
culverted(N,6) = culverted(N,6) + DELT
culverted(n,8) = 6
endif
cred type 7
IF(HW.GE.1.5*DEEP(N).AND.TW.LT.DEEP(N)) then
culverted(N,7) = culverted(N,7) + DELT
culverted(n,8) = 7
endif
endif
enddo
endif
RETURN
END
C EXTRAN BLOCK test for SWMM 5 beta solution - 12/12/2002
INCLUDE 'TAPES.INC'
INCLUDE 'STIMER.INC'
INCLUDE 'BD.INC'
INCLUDE 'BND.INC'
INCLUDE 'HYFLOW.INC'
INCLUDE 'CONTR.INC'
INCLUDE 'JUNC.INC'
INCLUDE 'PIPE.INC'
INCLUDE 'TIDE.INC'
INCLUDE 'OUT.INC'
INCLUDE 'ORF.INC'
INCLUDE 'WEIR.INC'
INCLUDE 'FLODAT.INC'
DOUBLE PRECISION AKON,QNEW,DELQ1,DELQ2,DELQ3,DELQ4,DELQ5
DIMENSION AS1(NEE)
DOUBLE PRECISION df,f,n_omega
integer good_nodes
C=======================================================================
C STORE OLD TIME STEP FLOW VALUES
C=======================================================================
DO N = 1,NTL
QO(N) = Q(N)
AT(N) = A(N)
VT(N) = V(N)
enddo
C=======================================================================
C INITIALIZE CONTINUITY PARAMETERS
C=======================================================================
DO J = 1,NJ
if(othercom(63).eq.1) then
Y(J) = Y(J) + 0.5 * ( YEX2(J) -
+ YEX1(J) + YEX1(J) - YO(J))
IF(Y(J).LT.FUDGE) Y(J)=FUDGE
IF(Y(J).GT.SURELEV(J)) Y(J)=SURELEV(J)-Z(J)
yex2(j) = yex1(j)
yex1(j) = YO(j)
endif
cred beginning time step value of the node area
asold(j) = as(j)
YO(J) = Y(J)
GoodNode(j) = .FALSE.
enddo
good_nodes = 0
omega = input_omega
n_omega = node_omega
loop_count = 0
C=======================================================================
C HALF-STEP AREA, RADIUS : VELOCITY
C FULL-STEP FLOW
C=======================================================================
big_loop: DO while(good_nodes.lt.nj.and.loop_count.le.itmax)
loop_count = loop_count + 1
if(loop_count.ge.itmax-5) then
omega = 0.50 * omega
n_omega = 0.50 * n_omega
endif
DO J = 1,NJ
AS(J) = AMEN
AS1(J) = 0.0
SUMQ(J) = QIN(J)
SUMQS(J) = QIN(J)
SUMAL(J) = 0.0
enddo
CIM FIRST COMPUTE GATED ORIFICE PARAMETERS
CALL OGATES(DELT,Y,V)
c
flow_loop: DO N = 1,NTC
NL = NJUNC(N,1)
NH = NJUNC(N,2)
C=======================================================================
H(N,1) = AMAX1(Y(NL) + Z(NL),ZU(N))
H(N,2) = AMAX1(Y(NH) + Z(NH),ZD(N))
CALL nhead(N,NL,NH,H(N,1),H(N,2),Q(N),A(N),V(N),HRAD,
+ ANH,ANL,RNL,RNH,YNL,YNH,width,IDOIT,LINK(N),AS1)
cred bypass loop for nodes already converged
bypass_loop: if(loop_count.gt.2.and.
+ goodnode(nl).eq..TRUE..and.goodnode(nh).eq..TRUE.) then
bypass = bypass + 1.0
else
IF(HRAD.GT.HMAX(N)) HMAX(N) = HRAD
IF(A(N).GT.AMAX(N)) AMAX(N) = A(N)
cred save information for the culvert classification
HRLAST(N) = HRAD
vup(n) = anl
vdn(n) = anh
if(loop_count.eq.1) then
aup(n) = anl
rup(n) = rnl
rmd(n) = hrad
endif
positive_flow: IF(IDOIT.gt.0) THEN
c
c Q/ANH = velocity at downstream end used for exit loss
c Q/ANL = velocity at upstream end used for entrance loss
c The loss = 1/2*K*V^2/g is factored into momentum
c equation similarly to the friction slope
c The loss momentum term = g*A*[1/2*K*V^2/g]
c = g*A*[1/2*K*Q/A*Q/A/g]
C =g *[1/2*K*|Q*A| /g] * Q
C = [1/2*K*|Q*A| ] * Q
c DELQ5 is sum of losses
c
cred calculate the froude number for every conduit
c
DELH = H(N,1) - H(N,2)
DELZP = ZU(N) - ZD(N)
DIFF_mid = 0.5 * (H(N,1) - ZU(N) + H(N,2) - ZD(N))
IF(DIFF_mid.GT.0.0) THEN
FROUDE_MID = ABS(Q(N))/A(N)/SQRT(GRVT*(DIFF_mid))
ELSE
FROUDE_MID = 0.0
ENDIF
bfactor = 0.0
if(FROUDE_MID.ge.1.0) THEN
bfactor = 0.0
elseif(froude_mid.gt.0.9) then
bfactor = 1.0 - (10.0*FROUDE_MID-9.0)
endif
cred test for zero sloped conduits
if(delzp.eq.0.0) then
del_ratio = delh/0.001
else
del_ratio = delh/delzp
endif
cred swmm 4 definition for normal flow
if(del_ratio.le.1.0) then
delfactor = 0.0
cred swmm 5 transition definition
else if(del_ratio.gt.1.10) then
delfactor = 1.0
else
delfactor = 10.0 * (del_ratio-1.0)
endif
DELQ5 = 0.0
IF(ANH.NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/ANH)*ENTK(N)
IF(ANL.NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/ANL)*EXITK(N)
IF(A(N).NE.0.0) DELQ5 = DELQ5 + 0.5 * ABS(Q(N)/A(N))*OTHERK(N)
DELQ5 = DELQ5 * DELT/LEN(N)
c
DELQ4 = DELT*V(N)*(ANH-ANL)/A(N)/LEN(N)
cred the mean area travels from the midpoint to the upstream area
cred as the value of delh/delzp changes
area_mean = wt*anl + wd*aup(n) + ( wt*a(n) + wd*at(n) -
+ wt*anl - wd*aup(n)) * delfactor
DELQ2 = DELT*GRVT*area_mean*((H(N,2) - H(N,1)))/LEN(N)
DELQ3 = 2.0*(A(N)-AT(N))/A(N)
cred the mean hydraulic radius travels from the midpoint to the
cred upstream hydraulic radius as the value of delh/delzp changes
hrad_mean = wt*rnl + wd*rup(n) + ( wt*hrad + wd*rmd(n) -
+ wt*rnl - wd*rup(n)) * delfactor
DELQ1 = DELT*(ROUGH(N)/hrad_mean**1.33333)*ABS(V(N))
QNEW = QO(N) - delq2
AKON = DELQ1 + DELQ5 - delq4*bfactor - DELQ3*bfactor
cred Newton-Raphson iteration
F = Q(N) * ( 1.0 + akon) - qnew
DF = 1.0 + akon
Q(N) = Q(N) - omega*f/df
DQDH = 1.0/(1.0+AKON)*GRVT*DELT*A(N)/LEN2(N)
dqdh_old(n) = dqdh
C=======================================================================
C DO NOT ALLOW A FLOW REVERSAL IN ONE TIME STEP
C=======================================================================
DIRQT = SIGN(1.0,QO(N))
DIRQ = SIGN(1.0,Q(N))
IF(DIRQT/DIRQ.LT.0.0) Q(N) = 0.001*DIRQ
v(n) = q(n)/A(N)
C=======================================================================
C COMPUTE CONTINUITY PARAMETERS
C=======================================================================
SELECT CASE (INGATE(N))
CASE (1)
Q(N) = AMAX1(Q(N),0.0)
CASE (2)
Q(N) = AMIN1(Q(N),0.0)
END SELECT
IF (NKLASS(N).LE.21) THEN
Q(N) = AMAX1(STHETA(N),Q(N))
Q(N) = AMIN1(SPHI(N),Q(N))
ENDIF
endif positive_flow
endif bypass_loop
SUMQS(NL) = SUMQS(NL) - Q(N)*barrels(n)
SUMAL(NL) = SUMAL(NL) + dqdh_old(n)*barrels(n)
SUMQS(NH) = SUMQS(NH) + Q(N)*barrels(n)
SUMAL(NH) = SUMAL(NH) + dqdh_old(n)*barrels(n)
SUMQ(NL) = SUMQ(NL) - WT*Q(N)*barrels(n) - WD*QO(N)*barrels(n)
SUMQ(NH) = SUMQ(NH) + WT*Q(N)*barrels(n) + WD*QO(N)*barrels(n)
ENDDO flow_loop
C=======================================================================
C SET FULL STEP OUTFLOWS AND INTERNAL TRANSFERS
C=======================================================================
CALL BOUND(Y,Y,Q,TIME,DELT)
C=======================================================================
N1 = NTC+1
DO 370 N = N1,NTL
NL = NJUNC(N,1)
NH = NJUNC(N,2)
IF(ABS(Q(N)).LT.1E-10) Q(N) = 0.0
C=======================================================================
SUMQ(NL) = SUMQ(NL) - WT*Q(N)*barrels(n) - WD*QO(N)*barrels(n)
SUMQS(NL) = SUMQS(NL) - Q(N)*barrels(n)
IF(NH.NE.0) THEN
SUMQ(NH) = SUMQ(NH) + WT*Q(N)*barrels(n) + WD*QO(N)*barrels(n)
SUMQS(NH) = SUMQS(NH) + Q(N)*barrels(n)
ENDIF
370 CONTINUE
C=======================================================================
C CALCULATE THE FULL-STEP HEAD
C=======================================================================
DO J = 1, NJ
AS(J) = AS(J) + AS1(J)
ENDDO
DO J = 1,NJ
IF(JSKIP(J).le.0) then
C=======================================================================
cred time weighted area average
if(othercom(62).eq.1) then
average_area = WT*as(j) + WD*asold(j)
else
average_area = as(j)
endif
cred Newton-Raphson iteration
yt(j) = y(j)
if(y(j).le.ycrown(j)) then
F = Y(J) * average_area - YO(J) * average_area
+ - sumq(j)*DELT
DF = average_area
ASFULL(J) = AS(J)
else
DF = sumal(j)
IF(Y(J).LT.1.25*YCROWN(J)) DF = DF +
+ (ASFULL(J)/DELT-SUMAL(J))*EXP(-15.*(Y(J)-YCROWN(J))/YCROWN(J))
cred correction from SWMM 5 QA testing - 11/12/2004
cred if a large delt or if there is a large value of sumal(j)
cred - usually from a very large conduit connected to the current node -
cred the expression ASFULL(J)/DELT-SUMAL(J) may be negative
cred invalidating the whole concept of a "transition" slot
if(Y(J).LT.1.25*YCROWN(J).le.ASFULL(J)/DELT.le.SUMAL(J))
+ denom = sumal(j)
CORR = 1.00
IF(NCHAN(J,2).EQ.0) CORR = 0.60
F = Y(J) * DF - YO(J) * DF - CORR * sumqs(j)
cred WRITE(*,*) ajun(j),df,f,y(j),ycrown(j),sumal(j)
endif
Y(J) = Y(J) - n_omega*F/DF
IF(Y(J).LT.0.0) Y(J) = 0.0
IF((Y(J)+Z(J)).GT.SURELEV(J)) Y(J) = SURELEV(J)-Z(J)
endif
enddo
cred converge until all nodes meet the convergence criteria
good_nodes = 0
error_node = 0
i_was_bad = 1
DO j = 1,nj
if(jskip(j).le.0) then
if(abs(y(j)-yt(j)).le.node_toler) then
good_nodes = good_nodes + 1
GoodNode(j) = .TRUE.
endif
if(abs(y(j)-yt(j)).gt.error_node) then
error_node = abs(y(j)-yt(j))
i_was_bad = j
endif
else
good_nodes = good_nodes + 1
GoodNode(j) = .TRUE.
endif
enddo
c write(*,669) loop_count,nj-good_nodes,error_node,
c + delt,ajun(i_was_bad),y(i_was_bad),
c + sumq(i_was_bad),omega
669 format(2I6,F10.4,F8.2,1x,a12,3F9.3)
enddo big_loop
cred culvert information - 8/6/2002
cred culvert information for the culvert comparison file - 8/6/2002
cred culvert information - 8/6/2002
if(othercom(90).eq.1) then
do n = 1,nc
if(nklass(n).eq.1.or.nklass(n).eq.21) then
if(abs(q(n)).gt.culvert_loss(n,5).and.abs(v(n)).lt.20.0)then
HW = H(N,1) - ZU(N)
TW = H(N,2) - ZD(N)
vup(n) = v(n)
vdn(n) = v(n)
head_loss = 0.5*ENTK(N)*vup(N)*vup(N)/GRVT +
+ 0.5*EXITK(N)*vdn(N)*vdn(N)/GRVT +
+ 0.5*OTHERK(N)*v(N)*v(N)/GRVT
sfloss = ROUGH(N)/GRVT *
+ abs(v(N))*ABS(v(N))/hrlast(n)**1.33333
cred sfloss = ROUGH(N)/GRVT *
cred + Q(N)*ABS(Q(N))/(A(N)**2*HRLAST(N)**1.33333)
culvert_loss(n,1) = len(n)*sfloss
culvert_loss(n,2) = head_loss
culvert_loss(n,3) = h(n,1)
culvert_loss(n,4) = h(n,2)
culvert_loss(n,5) = q(n)
endif
CALL DEPTHX(N,NKLASS(N),q(n),YC,YNORM)
cred type 1
IF(HW.LT.1.5*DEEP(N).AND.YC.LT.YNORM.and.
+ TW.LE.YC) then
culverted(N,1) = culverted(N,1) + DELT
culverted(n,8) = 1
endif
cred type 2
IF(HW.LT.1.5*DEEP(N).AND.YC.LT.YNORM.AND.
+ TW.GT.YC.AND.TW.LE.DEEP(N)) then
culverted(N,2) = culverted(N,2) + DELT
culverted(n,8) = 2
endif
cred type 3
IF(YC.GE.YNORM.AND.TW.LT.ZU(N)-ZD(N)) then
culverted(N,3) = culverted(N,3) + DELT
culverted(n,8) = 3
endif
cred type 4
IF(YC.GE.YNORM.AND.TW.GE.ZU(N)-ZD(N)+YC.
+ AND.TW.LT.ZU(N)-ZD(N)+DEEP(N)) then
culverted(N,4) = culverted(N,4) + DELT
culverted(n,8) = 4
endif
cred type 5
IF(YC.LT.YNORM.AND.TW.GE.DEEP(N)) then
culverted(N,5) = culverted(N,5) + DELT
culverted(n,8) = 5
endif
cred type 6
IF(YC.GE.YNORM.AND.TW.GE.ZU(N)-ZD(N)+
+ DEEP(N)) then
culverted(N,6) = culverted(N,6) + DELT
culverted(n,8) = 6
endif
cred type 7
IF(HW.GE.1.5*DEEP(N).AND.TW.LT.DEEP(N)) then
culverted(N,7) = culverted(N,7) + DELT
culverted(n,8) = 7
endif
endif
enddo
endif
RETURN
END
