Full Equations (FEQ) Model for the Solution of the Full, Dynamic Equations of Motion for One-Dimensional Unsteady Flow in Open Channels and Through Control Structures

U.S. GEOLOGICAL SURVEY WATER-RESOURCES INVESTIGATIONS REPORT 96-4240

13.1 Run Control Block--Run Control Table

Purpose: Parameters that are used to control the time span to be simulated, including specification of the convergence tolerances, the number of iterations allowed in the computations, and other run control information are supplied with this block.

Headings:Three lines of user-selected information to identify the input sequence.
LINE 1

Variable: NBRA

Format: 5X, I5

Example: NBRA=00005

Explanation: NBRA is the number of branches in the stream system. The maximum number of branches allowed in FEQ is specified in parameter MNBRA in the INCLUDE file ARSIZE.PRM (appendix 3). This number may be increased when necessary and FEQ recompiled.

LINE 2

Variable: NEX

Format: 4X, I5

Example: NEX=00010

Explanation: NEX gives the number of flow-path end nodes in the stream system. This number must be even because flow-path end nodes are always in pairs. The maximum number of flow-path end nodes allowed in FEQ is specified in the parameter MNEX in the INCLUDE file ARSIZE.PRM (appendix 3). This number may be increased when necessary and FEQ recompiled.

LINE 3

Variable: SOPER

Format: 6X, A4

Example: SOPER=YES

Explanation: SOPER=YES means that control structures are present, that they will be operated dynamically in FEQ simulation and, therefore, that the Operation of Control Structures Block (section 13.12) must be in the input sequence. SOPER=NO means that the Operation of Control Structures Block is not needed.

LINE 4

Variable: POINT

Format: 6X, A4

Example: POINT=NO

Explanation: POINT=NO means that no point flows are present and that the Point Flows Block (section 13.7) will not be used. POINT=YES means that point flows are present and the Point Flows Block will be used.

LINE 5 Release information available for HECDSS support, DIFFUS = DSS

Variables: DIFFUS, MINPRT, LAGTSF, DMYEAR, DMMN

Format: 7X, A4, 2I1, 2I5

Example: DIFFUS=YES 10 1902 1

Explanation:

DIFFUS=YES means that diffuse inflows from the land area tributary to the stream system are used in the simulation. A Tributary Area Block (section 13.3) must be present in the primary input, and a DTSF in the format required in FEQ (appendix 2) must be present in the disk space for computations. DIFFUS=NO means that diffuse inflows will not be included in the unsteady-flow simulations, and no Tributary Area Block should be included in the primary input file.

MINPRT is a variable used to control the output of the normal summary information given at the completion of each time step when DIFFUS=YES. This output is not controlled by the variable PRTINT specified on Line 18. If MINPRT=0, then the summaries are printed; the summary output is suppressed if MINPRT=1. Error and warning messages are never suppressed.

LAGTSF is a time-lag variable applied to data in a time-series file. If LAGTSF is nonzero, the data in the time-series file are lagged; the result is a greater delay in water reaching the stream channels. This option is a simple means for approximating the effect of less abrupt variation of the lateral inflows to the stream. The lagging is always by one time interval of the time-series file. Instead of each runoff value for each time interval being treated as a constant, the runoff is lagged so that, for an isolated interval with runoff, the runoff increases from zero to the peak value over the interval. The runoff then decreases to zero from the peak over the subsequent time interval. Thus, the volume is not changed but the average delay time (that is, the time measured from centroid of runoff before and after lagging), is one-half the time interval. Adjacent time intervals are treated in the same manner, and superposition is used in combining the overlapping flows that result. This adjustment may not be realistic, but it can applied as a quick test to see whether greater refinement of the distribution of lateral inflows will increase the robustness of the computations.

DMYEAR is the year for the dummy event. If omitted, the year for the dummy event defaults to 1925.

DMMN is the month for the dummy event. If omitted, the month for the dummy event defaults to 1.

The interpretation of the starting and ending times (given on Lines 9 and 10) is affected by the DTSF. The DTSF contains one or more disjoint time segments (called events) containing the unit-area runoff intensities for one or more land-cover types or land uses. These unit-area runoff intensities are multiplied by areas specified in the Tributary Area Block (section 13.3) to derive the inflow from these land uses into the stream system. The events in the DTSF are given in ascending order of time. The first event must always be simulated to establish the initial conditions to be used for all subsequent events. The subsequent events are independent of each other and need not be simulated in any particular sequence. Thus, one can skip events and simulate only the event of interest. To do this, a starting time equal to or less than the starting time of the event of interest must be specified. The first event that has a starting time equal to or greater than the specified start time will be simulated. The ending time also is checked and computations will end at the ending time specified even if that time is within one of the events stored in the DTSF. The computations stop at the end of an event in the DTSF or at the end time specified by the user. Thus, if the entire DTSF is to be simulated, an end time known to be beyond the end time of the last event stored in the DTSF should be specified.

The information in the DTSF can come from any source, including continuous rainfall-runoff models such as the many variations of the Stanford Watershed Model, HEC-1, or one of the procedures from the Natural Resources Conservation Service. Data must be in the proper format in FEQ. The details of the format are presented in appendix 2, "Unformatted Data-File Structures."

LINE 6

Variable: WIND

Format: 5X, A4

Example: WIND=NO

Explanation: WIND=NO means that the effect of wind-shear stress on the water surfaces of the streams will not be simulated. WIND=YES means that the wind-shear stress on the water surfaces of the streams will be simulated, and a Wind Information Block (section 13.8) must be included in the primary input file in this case. Wind-shear stress is computed for branches only and cannot be computed for level-pool reservoirs.

LINE 7

Variable: UNDERFLOW

Format: 10X, A4

Example: UNDERFLOW=NO

Explanation: This variable is used to suppress underflow messages and should always be UNDERFLOW=NO. This variable is not currently used in FEQ, but it was previously needed for some compilers and operating systems. Some compilers and operating systems, by default, issue messages for each computational underflow resulting during the execution of a program. Many thousands of underflow messages may be generated in FEQ simulation. Issuance of these messages only slows the computations and clutters the output file. These messages can be suppressed by calling a compiler- and operating-system-specific routine. Therefore, this line is retained for possible future use.

LINE 8

Variable: ZL

Format: 3X, F10.0

Example: ZL=1.200

Explanation: ZL is the cutoff depth for applying the zero-inertia approximation in a computational element to simplify computation of supercritical flows. Inertia terms in the governing equations are suppressed whenever the depth at either end of the computational element is less than ZL. Unless special problems arise, ZL should be set to zero.

LINE 9

Variables: SYR, SMN, SDY, SFRAC

Format: 6X, I4, 1X, I2, 1X, I2, 1X, F12.0

Example: STIME=1980/01/20:12.0

Explanation:

The starting time for the current simulation is specified with this line. SYR is year (all four digits),

SMN is number of the month (two digits),

SDY is day of month (two digits), and

SFRAC is hour of the day.

The example gives noon on the 20th day of January 1980 as the starting time. The year should be 1859 or later for internal date/time computations to function properly.

LINE 10

Variables: EYR, EMN, EDY, EFRAC

Format: 6X, I4, 1X, I2, 1X, I2, 1X, F12.0

Example: ETIME=1981/12/31:14.56

Explanation: The ending time for the current simulation is specified with this line. Values are defined analogously to the starting times on Line 9.

LINE 11 Release information available for metric units.

Variable: GRAV

Format: 5X, F10.0

Example: GRAV=32.2

Explanation: The acceleration due to gravity is specified in the proper system of units. Two systems of units are supported in FEQ. If GRAV=32.2, then English units are used. In this case, the internal units are lengths in feet, areas in square feet, volumes in cubic feet, and times in seconds. Volume units in the English system are in acre-feet when printed. SFAC (specified on Line 24) is used to convert cross-section stations from the unit used in the input to the internal model unit. The square of SFAC is used to convert tributary area from the input unit to the internal unit. Tables prepared and input by the user may be in any units the user prefers. A multiplying factor specified by the user is applied to convert the values input to the internal units for FEQ simulation. For example, storage values for level-pool reservoirs may be input in acre-feet with an accompanying multiplying factor of 43,560 to convert acre feet to cubic feet. The unit of the unit-area runoff intensity from a diffuse time-series file must be consistent with the internal units; that is, feet per second. Example tables may be obtained by electronic retrieval from the World Wide Web (WWW) at http://water.usgs.gov/software/feq.html and by anonymous File Transfer Protocol (FTP) from water.usgs.gov in the pub/software/surface_water/feq/ directory.

If GRAV= 9.8 the internal units are meters for lengths, square meters for areas, cubic meters for volumes, and time in seconds. The volume unit for printing volumes is 1,000 cubic meters. Again, multiplying factors are available for making conversions between the units in the input and the internal units.

LINE 12

Variable: NODEID

Format: 7X, A4

Example: NODEID=YES

Explanation: NODEID=YES means that an identification string of as many as eight characters may be given in the Branch Tables for each node on a branch. Free nodes described in the Free-Node Initial Conditions Block (section 13.14) also may have an identification string of the same length. The string may be left blank if no printed node identification is desired. The identification string is printed with the node when results are reported. NODEID=NO means that no identification string can be given with the nodes. NODEID=YES should be used because clearer output is produced if node names are carefully selected.

LINE 13

Variable: SSEPS

Format: 6X, F10.0

Example: SSEPS=0.10

Explanation: SSEPS is the convergence tolerance for the volume of water in ponding storage when sewers are simulated. The maximum relative change in surcharge storage volume allowed anywhere in the system during the iterative computation of the surcharge storage is specified with SSEPS. If SSEPS is made too small, convergence problems will result. If SSEPS is too large, inaccurate ponding volumes may result.

LINE 14

Variables: NAME, PAGESP

Format: A4, 1X, I5

Example: PAGE=00022

Explanation:

NAME is the identifying character string.

PAGESP is the number of lines per page for the special output file. The maximum number of lines per page in the special outfile allowed in FEQ is specified in parameter MSPKNT in the INCLUDE file ARSIZE.PRM (appendix 3). This value should probably not be changed.

LINE 15

Variable: EPSSYS, ABSTOL, FAC

Format: 7X, 3F5.0

Example: EPSSYS=0.0500.010

Explanation:

The convergence tolerances for the matrix solution are specified in this line.

EPSSYS is the primary relative tolerance specifying the maximum relative change in the value of an unknown (flow or water-surface height depending on context) between iterations. The maximum relative change as computed is based on the last correction from Newton's method. The last correction is made so that the relative error in the unknowns is probably much less than the stated tolerance. Therefore, the tolerance should not be made too small. Some simulation experiments have shown that a change from 0.05 to 0.005 for EPSSYS changed the computed maximum water-surface elevations by about 0.02 ft or less. The flows differed by a few percent at most. The execution time with EPSSYS=0.05 was about 30 percent less than the execution time with EPSSYS=0.005.

ABSTOL is a tolerance on the change in water-surface height at nodes on branches in the stream system. If the absolute value of the change is less than ABSTOL, then the relative change for that variable is treated as if it were zero. ABSTOL must be nonzero if small depths are to be simulated. ABSTOL provides a criterion for small depths similar to that of QSMALL for low flow rates. The analogous tolerance for level-pool reservoirs, RESEPS, is set to 0.001 ft.

FAC is a factor used to establish secondary relative tolerance. If input for FAC is omitted, then FAC is set at 2.0. The secondary tolerance is computed by multiplying the primary relative tolerance (EPSSYS) by FAC. This secondary relative tolerance is used when NUMGT (given in next input line) is positive. Convergence is obtained if all but NUMGT variables satisfy the primary relative tolerance and the NUMGT variables that exceed the primary relative tolerance are all less than the secondary relative tolerance. The secondary relative tolerance can sometimes be used to increase the robustness of the computations in FEQ.

Variables: MKNT, NUMGT

Format: 5X, 2I5

Example: MKNT=00005 1

Explanation:

MKNT is the maximum number of iterations for the matrix system solution. If MKNT is exceeded, then the time step is reduced by multiplication with the LFAC value given on Line 28 and convergence is attempted again. This process continues until either convergence is obtained or the time step becomes less than the minimum time step given on Line 28. The maximum number of iterations allowed in FEQ is specified in the parameter MNITER in the INCLUDE file ARSIZE.PRM (appendix 3). The value of MNITER should be changed only as a last resort to obtain computational convergence. If MNITER is changed, then FEQ must be recompiled.

NUMGT is the number of nodes at which the secondary relative convergence criterion may be applied. Use of NUMGT solves a frequent problem of robustness. Convergence of the computations will frequently be prevented by problems at only one or two nodes; all other nodal variables may meet the convergence criterion for relative change by a large margin, yet convergence cannot be declared for the system until all nodal variables meet the criterion. Test simulations in which statistics on the relative change for the nodal variables were calculated and printed indicated that, in many cases, problems in only one or a few nodes can greatly delay or prevent convergence. Convergence can be obtained if the number of variables not meeting the current primary convergence criteria does not exceed NUMGT and if these variables satisfy the secondary relative tolerance defined in Line 15. The number of variables not meeting the current convergence criteria is printed in the iteration log. The value NUMGT should be small. A value of 1 generally results in greatly increased robustness and only minor effects on the final computed flows and depths. In some cases, the only way to obtain convergence without greatly modifying the stream-system representation is to apply NUMGT. The iteration log should be checked to determine when, where, and how often NUMGT is used in simulation. In most cases, relative changes will not exceed the convergence criterion. Often, cases in which the criterion is exceeded at one node result during low-flow periods and will have little effect on the higher flows and depths of interest.

Variables: OUTPUT, PROUT

Format: 7X, I5, 1X, I5

Example: OUTPUT=00001

Explanation:

OUTPUT controls the level of detail in the output. Setting OUTPUT equal to zero results in the minimum output detail, and setting OUTPUT equal to 5 results in the maximum output detail for debugging of simulation errors. A value of zero is used in normal operation. A value of 1 results in a complete echo of the cross-section tables in addition to the headings as for OUTPUT=00000. A dump of some internal table values for checking purposes also results for OUTPUT=00001. For values of OUTPUT greater than 1, the actual computations are detailed in the output from FEQ. OUTPUT=00005 should be applied sparingly because large amounts of output are produced. If OUTPUT is greater than 1, then Lines 19 and 20 must be in the input sequence.

PROUT specifies the output level applied after the simulation time has passed the optional date/time given after PRTINT on Line 18.

Variable: PRTINT, YR, MN, DY, HR

Format: 7X, I5, 1X, I4, 1X, I2, 1X, I2, 1X, F10.0

Example: PRTINT=00005

Explanation:

PRTINT controls the frequency of output of the stream-system results. PRTINT=00001 means that output is printed every time step, PRTINT=00002 means that output is printed every second time step, PRTINT=00003 means that output is printed every third time step, and so forth. The value of PRTINT is set to 1 so that the complete results of every time step are output. (This setting should be used only for debugging of simulation errors.)

YR is the four-digit year,

MN is the month number,

DY is the day number, and

HR is the hour of the day with possible fractional hours.

Variable: STPRNT

Format: 7X, I5

Example: STPRNT=00010

Explanation: STPRNT is the number of the time step at which detailed output is to start when OUTPUT is greater than 1. A more useful option for debugging is application of the options on Lines 17 and 18.

LINE 20

Variable: EDPRNT

Format: 7X, I5

Example: EDPRNT=00100

Explanation: This variable specifies the number of the time step at which detailed output is to stop when OUTPUT (Line 17) is greater than 1. A more useful option for debugging is application of the options on Lines 17 and 18.

LINE 21

Variables: NAME7, NAMEGE

Format: A7, A5

Example: GEQOPT=STDX

Explanation:

NAME7 is the identifying character string "GEQOPT=." The string GEQOPT= must be used because both parts of the input are read; if the first six characters are not GEQOPT, then a warning message is issued and the option STDX is selected.

NAMEGE is the default Governing EQuation OPTion for the branches in the stream system. The options for the governing equations are: STDX, STDW, STDCX, and STDCW. These options are described in detail in section 6.2. If STDW or STDCW are selected, then the Near-Zero-Depth Block (section 13.5) must appear in the input immediately before the Network-Matrix Control Input.

Variable: EPSB

Format: 5X, F10.0

Example: EPSB=0.00001

Explanation: EPSB is the convergence tolerance for the steady-state flow analysis used to define the initial conditions in the stream system. This value should be small to produce small residuals in the steady-flow governing algebraic equations because this residual is, in effect, multiplied by the time step in the unsteady-flow solution. A large residual may require extra computational effort during the frozen-time stage of the FEQ computations.

LINE 23

Variable: MAXIT

Format: 6X, I5

Example: MAXIT=00020

Explanation: This variable specifies the maximum number of iterations permitted in the steady-state flow analysis to define the initial conditions in the stream system.

LINE 24

Variable: SFAC

Format: 5X, F10.0

Example: SFAC=5280.0

Explanation: SFAC is the multiplying factor for converting the stations for the nodes along the stream system given in the Branch Description Block (section 13.2) to the internal units. The combination of SFAC and GRAV (Line 11) determine the internal units. If GRAV=32.2, English units are used and internal lengths are in feet, areas are in square feet, and volumes are in cubic feet. Time is in seconds. For convenience, SFAC is used to allow different units for the stations along the stream and for tributary areas. SFAC is applied to multiply whatever values are used for the stations, and SFACSFAC is used to adjust the tributary areas. If GRAV=32.2 and SFAC=5,280, then stations along the stream are given in miles and tributary areas are given in square miles.

LINE 25

Variables: QSMALL, QCHOP

Format: 7X, 2F10.0

Example: QSMALL=50.00.01

Explanation:

QSMALL is the value of flow rate to add to a flow before computing a relative change in the flow. Use of QSMALL prevents a division by zero if the flow is zero. QSMALL should be small relative to the flows of interest, but it should not be made too small because unneeded computational effort is required for reducing the relative changes between iterations to a value smaller than needed. A value on the order of a few percent of the maximum flows expected is generally reasonable.

QCHOP is the truncation level for output of flows. Any flow less than or equal to QCHOP is set to zero. This option is used to eliminate the output of small flows that result from roundoff during the computations. The default value for QCHOP is 0.001.

Variable: IFRZ

Format: 5X, I5

Example: IFRZ=00005

Explanation: Losses from flow through most control structures are not currently included in the steady-flow analysis used to determine initial conditions for simulation of unsteady flow. Therefore, simulated time is "frozen" for the number of time steps given by IFRZ to yield initial conditions that do include these losses. Because time is frozen, all the inflows and outflows that are functions of time remain constant. The following line gives the time steps to be used during this process. The concept of frozen time is explained in detail in section 8.3. The maximum number of frozen time steps allowed in FEQ is specified in the parameter MNFRDT in the INCLUDE file ARSIZE.PRM (appendix 3). The number may be changed, but FEQ must then be recompiled.

LINE 27

Variables: DTVEC(*)

Format: 10F8.0

Example: 1800.1250.850.600.400.275.200.150.100.

Explanation: The IFRZ time-step values to apply while time is frozen are specified on this line. The first time step is given last, and the last time step is given first. Normally, the first time step in frozen time should be small, and the time step should become progressively larger for each frozen-time period simulated. This line is not present if IFRZ=0.

LINE 28

Variables: MAXDT, MINDT, AUTO, SITER, HIGH, LOW, HFAC, LFAC

Format: 8F5.0

Example: 18000.10.72.93.22.42.00.5

Explanation:

The parameters used to select time steps are specified on this line.

MAXDT is the maximum time step allowed in FEQ simulation of unsteady flow.

MINDT is the minimum time step allowed in FEQ simulation of unsteady flow.

AUTO is a weighting factor used to compute a weighted running sum of the number of iterations, which is applied to detect convergence problems. AUTO must be in the range of 0 to 1.

SITER is an initial value for the computation of the weighted running sum of the number of iterations.

HIGH is the upper bound on the weighted running sum of the number of iterations. If the running sum is greater than HIGH, then the time step is decreased.

LOW is the lower bound on the weighted running sum of the number of iterations. If the running sum is less than LOW, then the time step is increased.

HFAC is a multiplier used to increase the computational time step to improve computational convergence.

LFAC is a multiplier used to decrease the computational time step to improve computational convergence.

In order to detect convergence problems, a weighted running sum, SITER, is computed from the number of iterations to reach convergence or to exceed MKNT iterations. The running sum is updated by use of the relation SITER=AUTOSITER + (1 - AUTO)ITER, where ITER is the current iteration count. To keep SITER from growing without bound, it is required in FEQ that 0 AUTO < 1.0. As an example, assume that the number of iterations, ITER, remains constant for many time steps. In this case, the value of SITER will approach ITER, and the speed of approach will depend on the value of AUTO. If AUTO is 0.0, then SITER is the number of iterations to convergence for the previous time step. This value of SITER may result in time-step changes that are too rapid. In contrast, if AUTO is nearly 1.0, the value of SITER will change only slowly and convergence will be obtained slowly.

With SITER determined, the time step is unchanged if LOW SITER HIGH. If SITER > HIGH, then the time step is multiplied by LFAC to reduce it. If SITER < LOW, then the time step is multiplied by HFAC to increase it. The goal is to keep the number of iterations at about 2. Keeping time steps so large that the number of iterations increases is not useful because convergence of the iterations does not mean that an accurate solution has been obtained. Suggested values are AUTO=0.7, HIGH=3.2, LOW=2.4, HFAC=2.0, and LFAC=0.5. The value given for SITER should be midway between HIGH and LOW. The maximum time step allowed, MAXDT, should be selected on the basis of the expected rate of change of flow during the simulation. If MAXDT is too large, major events may be missed; and if MAXDT is too small, unneeded computations will result. MINDT should be set small so that short-duration computational problems will not result in premature termination of the simulation. Values of MINDT of 1 second or less are suggested from modeling experience.

LINE 29

Variables: MRE, FAC1

Format: 4X, 5F10.0

Example: MRE=0.1 1.0

Explanation:

MRE is the maximum relative change permitted in a variable during extrapolation.

FAC1 is a factor used to regulate the degree of extrapolation from previous values at each node to aid in the solution of the equations. If FAC1=1.0, then full linear extrapolation in time from the two previous values at each node is used to estimate the unknowns at the end of the current time step. If FAC1=0.5, one-half the extrapolation is used and so on. If diffuse flows are included in the

simulation and are dominant, then the best results are obtained when FAC1=0. If diffuse flows are not included in the simulation, then use of FAC1=1.0 often reduces the time required to complete the computations.

Variable: DWT

Format: 4X, F10.0

Example: DWT=0.1

Explanation: DWT is the change in the current value of W _T, the weight factor for approximating integrals with respect to time. DWT and BWT (described in Line 31) are used to vary the value of W_T depending on the conditions simulated in the computations. Whenever the time step must be reduced because convergence was not obtained within MKNT iterations, the current value of W_T is incremented by DWT. If the resulting value for W_T >1.0, then W_T =1.0. Whenever the time step is decreased because SITER < LOW, then the current value of W_T is decreased by DWT. If the resulting value of W_T < BWT, then W_T = BWT. Varying W_T during the computations balances the need for accuracy and the need to damp nonlinear oscillations. A fixed value of W_T can still be applied throughout the simulation by setting DWT=0.0.

LINE 31

Variable: BWT

Format: 4X, F10.0

Example: BWT=0.5

Explanation: BWT is the base value of the weight factor for approximating integrals with respect to time. The valid range for this factor is 0.5 BWT 1.0. Values in the range of 0.55 to 0.7 represent a balance between accuracy and damping of nonlinear oscillations.

LINE 32

Variables: BWFDSN, BWFNAM

Format: 7X, I2, A64

Example: BWFDSN=15

Explanation:

BWFDSN is the Fortran unit number (see appendix 3) for the file used to save the initial conditions when DIFFUS=YES. A value is read but not used in FEQ when DIFFUS=NO.

BWFNAM is the name for the file opened in FEQ simulation to store the initial conditions needed for events in the DTSF. If BWFNAM is not blank, then a file will be opened with a name given by the contents of BWFNAM (in most microcomputers). On IBM mainframes, the name given in BWFNAM will be the ddname for the DD statement defining the data set. If BWFNAM is blank, then an implicit open procedure will be done on an IBM mainframe if the proper DD statement defining the unit number given by BWFDSN is present. In some microcomputer systems, the user will be prompted for the file name in this case, whereas execution of FEQ will be aborted on other microcomputer systems. This file is written at the end of the dummy event in the DTSF. The dummy event is always the first event and usually has been set in 1925. If an event prior to 1925 is simulated, then the date for the dummy event must be given on Line 5.

Variable: CHKGEO

Format: 7X, A4

Example: CHKGEO=YES

Explanation: CHKGEO is used to select geometric checking. If geometric checking is selected, then execution of FEQ is terminated at the end of the checking process so that changes to the input can be made.

The unsteady-flow computations in FEQ are often sensitive to the variations in channel geometry. If the geometric variations among cross sections are too large, then a solution may not be obtained or an invalid solution may result, producing very large Froude numbers that have no basis in physical reality. The governing equations are approximated by use of simple rules of numerical integration in each computational element. If the changes in top width, area, or conveyance from the upstream node to the downstream node of the element are too large, then large errors in the approximations are introduced. The list of depths that appear in the cross-section tables at each end of each computational element in the system are checked in CHKGEO. For each of these depths, the ratio of the downstream hydraulic-characteristic value to the upstream hydraulic-characteristic value is computed, where the hydraulic characteristics are top width, area, and conveyance. If this ratio is too small or too large, then a

message is printed alerting the user to the value of the ratio.

The absolute change in the cross-sectional hydraulic characteristics from one end of the computational element to the other is important, but the rate of change of these elements with distance at a constant depth also is important. Natural channels commonly include rapidly expanding or contracting reaches. The Manning's roughness values used to define the conveyance of the channels normally relate to the losses resulting in straight reaches. Additional losses must be included in reaches that are rapidly expanding or contracting. This input is considered in FEQ; the rate of change of top width and the rate of change of area, normalized to be in the same units as the rate of change of top width, is computed in CHKGEO as an indicator of the need to compute additional losses. If the rates of change are too large or too small, then a message is printed alerting the user to the value of the ratio.

LINE 34

Variable: ISTYLE

Format: 7X, A4

Example: ISTYLE=OLD

Explanation: The style of labeling flow-path end nodes is selected with ISTYLE. If ISTYLE=OLD, then the flow-path end nodes are numbered and the Branch-Exterior Node Block (section 13.4) must appear in the input. This is applied in versions of FEQ prior to 7.0. If ISTYLE=NEW, then the flow-path end nodes are labeled with a prefix character and the branch number if on a branch. An upstream branch end node is labeled with the prefix letter "U" followed by the branch number. A downstream branch end node is labeled with the prefix letter "D" followed by the branch number. Labels for level-pool reservoir and dummy-branch end nodes are prefixed by the letter "F" followed by a number chosen by the user in the range 1 to 999. As an example, the upstream end node on branch 231 is U231, and the downstream node is D231. A level-pool reservoir node and the associated inflow node could be labeled F500 and F499, respectively. The even number denotes the downstream end of the flow path represented by the level-pool reservoir.

LINE 35

Variable: EXTTOL

Format: 7X, F10.0

Example: EXTTOL= 0.5

Explanation: EXTTOL is the distance above the top of section that can be extrapolated during the unsteady-flow computations before a warning message is issued at the completion of the simulation. The extension of each cross section by the EXT option is checked to determine whether the maximum elevation of water exceeds the unextended maximum level of the cross-section table by more than EXTTOL length units.

LINE 36

Variable: SQREPS

Format: 7X, F10.0

Example: SQREPS=1.E25

Explanation: "SQREPS=" is always checked; therefore, it must be spelled as given here. SQREPS is used to solve special problems with Newton's method. An application of Newton's method will sometimes oscillate about a solution at one or more nodes that are close to a point of inflection on a n- dimensional solution surface, where n is the number of unknowns. In this case, the correct direction for the correction has been determined in Newton's method, but the magnitude of the correction is too large. A sketch of an example of this problem with one unknown is shown in figure 18B. The sum of the squares of the residuals is computed from all of the residual functions in the nonlinear system of equations describing the flow in the model of the stream system. If the magnitude of the correction computed from Newton's method is small enough, the sum of squares of the residuals should be smaller from iteration to iteration. The tolerance on the change in the sum of squares is given in SQREPS. If SQREPS=00000, then the sum of squares must stay the same or decrease to avoid reducing the Newton correction in simulation to what is called a partial Newton correction. If SQREPS=10000, then the sum of squares must increase by more than 10,000 before a partial Newton correction is applied.

The recommended procedure is to suppress partial Newton corrections until necessary. Therefore, if the value for SQREPS is left blank, a default value of 1.E30 is assigned, effectively suppressing any partial Newton corrections. Partial Newton corrections may be helpful if the maximum relative correction is at the same node and if the solution is close to convergence for successive iterations at that node. The corrected values of the sum of squares also should oscillate, at least approximately, between two values, both of which are close to but still not within the convergence criteria. The time step will continue to be reduced until the computations are terminated because the minimum-time-step limit has been reached. In this case, reducing SQREPS to a smaller value might circumvent the convergence problem. The iteration log, printed for each time step, includes (under the heading SUMSQR) the sum of squares of the residuals for the equation system. Examining these values for several time steps for which the computations did converge can help determine what size of tolerance would likely result in convergence. Too small a tolerance can result in unneeded partial Newton corrections and can lead to failure of convergence. The sum of squares does not always get smaller from iteration to iteration, especially in the earlier iterations for a time step. Therefore, partial Newton corrections should be used only in those cases that clearly indicate oscillation about a root.

LINE 37

Variables: NAME2, GETDSN, GETNAM

Format: A6, I5, A64

Example: GETIC=00009START

Explanation:

The file used to define the initial conditions for starting the computations is specified in this line. This file will have been output from FEQ with the PUTFC option in a previous simulation (described in Line 38).

NAME2 is the identifying character string "GETIC=." The string "GETIC=" must appear and must be spelled as shown.

GETDSN is the Fortran unit number (see appendix 3) for the file containing a previous FEQ simulation to be used as the initial condition for the current FEQ simulation. When the value of GETDSN is greater than zero, the initial conditions will be given in the file named after the number.

GETNAM is the name of the file containing the previous FEQ simulation to be used as the initial condition for the current FEQ simulation. If no file name is given (useful primarily on IBM mainframes), then the file is implicitly opened in FEQ following the IBM Fortran conventions. If the name is given and an IBM mainframe is being used, then the DD name is used to access the file. If this feature is not used, then the input following GETIC= should be left blank.

Variables: NAME3, PUTDSN, PUTNAM

Format: A6, I5, A64

Example: PUTFC=00009START

Explanation:

Saving the final conditions computed in FEQ for the simulation in the named file is specified with this line.

NAME3 is the identifying character string "PUTCF=." The string "PUTFC=" must appear and must be spelled as shown.

PUTDSN is the Fortran unit number (see appendix 3) for a file containing the final conditions computed in FEQ simulation to be used as the initial condition for a subsequent FEQ simulation. When the value of PUTDSN is greater than zero, the final conditions will be output to the file named after the number.

PUTNAM is the name of file in which the final conditions computed in FEQ simulation will be placed to be used as the initial condition for a subsequent FEQ simulation. The interpretations of the file name on various computer systems are listed in the explanation for Line 37.

A further and more common use for the PUTFC and GETIC options is to use them to save an initial condition from which to start FEQ simulation. The computations are often difficult to start. Once the computations are started the boundary conditions are held fixed, by use of dummy input tables or files, until transients from the assumed conditions damp out. Applying the PUTFC option at the end of the run then saves these results for later access by use of GETIC. An initial condition can be shifted by starting from a known initial condition and specifying the input hydrographs needed to attain the final condition of interest. Saving this final condition then makes it available to start subsequent simulations. This applies only if the stream-system model does not change. Addition or deletion of flow paths invalidates the saved initial conditions, but other changes, such as to an inflow hydrograph, do not.