ICM InfoWorks and ICM SWMM RED (Rainfall Event Data) file format

 This document describes the RED (Rainfall Event Data) file format,  used as an input for hydrological simulations, within the InfoWorks ICM software suite or other related tools.  

File Structure:

The RED file is a structured text file containing rainfall data and associated parameters. It consists of the following record types:

• Records 60-63 (Optional): These records, if present, must be at the beginning of the file. They define analysis parameters, ground condition parameters, and unit hydrograph parameters, that are used in a simulation or a batch set of simulations.
• Record 20 (Title): Contains a title or description of the rainfall data. There is only one Record 20 per file.
• Record 21 (Global Variables): Defines global parameters for the rainfall event, such as event number, start date/time, rainfall time step, number of rainfall profiles, and a global Urban Catchment Wetness Index (UCWI). There is only one Record 21 per file.
• Record 22 (Rainfall Event Data): Contains data for each rainfall profile, including a local UCWI (overriding the global value if not 0), an areal reduction factor, and a profile title. There is one Record 22 for each profile defined in the file.
• Record 23 (Rainfall Hyetographs): Contains the actual rainfall intensities (in mm/hr) for each time step and for each profile. There are multiple Record 23s, one for each time step, with each record containing data for all profiles.

Record Details:

• Record 20 (Title):

  • Column 1: Flag (must be 0)
  • Columns 2-80: Title (text)

• Record 21 (Global Variables):

  • Columns 1-5: Event Number (integer)
  • Columns 7-14: Start Date (dd/mm/yy or ddmmyyyy)
  • Columns 16-20: Start Time (hh:mm, 24-hour clock)
  • Columns 21-25: Rainfall Time Step (seconds, integer)
  • Columns 26-30: Number of Profiles (integer)
  • Columns 41-50: Urban Catchment Wetness Index (UCWI, decimal with 1 decimal place)

• Record 22 (Rainfall Event Data):

  • Columns 1-10: Local UCWI (mm, decimal with 1 decimal place. If 0, the Record 21 value is used).
  • Columns 21-30: Areal Reduction Factor (decimal with 1 decimal place)
  • Columns 41-80: Profile Title (text)

• Record 23 (Rainfall Hyetographs):

  • Columns 1-80: 10 columns of 8 characters each, representing rainfall intensities (mm/hr, decimal with 2 decimal places) for each profile. Each record represents a single time step.

• Record 60 (Analysis Information):

  • Columns 1-10: Analysis Interval (minutes, integer)
  • Columns 11-20: Analysis Duration (minutes, integer)
  • Columns 21-30: Antecedent Dry Weather Period (ADWP, hours, decimal with 3 decimal places)
  • Columns 31-40: Evaporation (mm/day, decimal with 3 decimal places)

• Record 61 (Ground Conditions):

  • Column 1: Profile Number (integer)
  • Columns 2-10: API 30 (Antecedent Precipitation Index for 30 days, mm, decimal with 3 decimal places)
  • Columns 11-19: depPaved (Depth of depression storage on paved areas, mm, decimal with 3 decimal places)
  • Columns 20-28: depPerv (Depth of depression storage on pervious areas, mm, decimal with 3 decimal places)
  • Columns 29-37: NAPI for Soil Type 1 (mm, decimal with 3 decimal places)
  • Columns 38-46: NAPI for Soil Type 2 (mm, decimal with 3 decimal places)
  • Columns 47-55: NAPI for Soil Type 3 (mm, decimal with 3 decimal places)
  • Columns 56-64: NAPI for Soil Type 4 (mm, decimal with 3 decimal places)
  • Columns 65-73: NAPI for Soil Type 5 (mm, decimal with 3 decimal places)

• Record 62 (NAPI Limiting Factors):

  • Column 1: Profile Number (integer)
  • Columns 2-10: Reserved.
  • Columns 11-80: Maximum cumulative NAPI values for Soil Types 1-5 (mm, decimal with 3 decimal places), one value per 9-column field.

• Record 63 (New Unit Hydrograph Parameters):

  • Column 1: Profile Number (integer)
  • Columns 2-10: Initial Soil Moisture (Cini, mm, decimal with 3 decimal places)
  • Columns 11-80: Reserved for future use.

Key Features and Implications:

• Multiple Profiles: The format supports multiple rainfall profiles within a single file, allowing for spatial variation of rainfall or different storm scenarios to be represented.
• Time-Series Data: Record 23 provides time-series rainfall data, defining the temporal distribution of rainfall intensity.
• Antecedent Conditions: Records 21, 60, 61, and 62 allow for the specification of antecedent conditions (UCWI, ADWP, API, NAPI, depression storage), which can significantly influence runoff generation.
• New Runoff and Unit Hydrograph Methods: Records 61, 62 and 63 specifically support newer methods for runoff calculation and unit hydrograph generation (likely ReFH, which stands for Revitalized Flood Hydrograph method).
• Batch Processing: Record 60 suggests that the format is designed to support batch processing of multiple simulations.

In essence, the RED file format provides a structured way to input rainfall data and associated parameters into a hydrological model, supporting multiple rainfall profiles, antecedent conditions, and newer hydrological methods. The format's flexibility makes it suitable for a range of hydrological modeling applications.

Summary of dwflow.c in SWMM5

Below is a detailed overview of dwflow.c, which is the portion of SWMM5 that handles dynamic wave hydraulic computations inside a single conduit. It specifically focuses on how SWMM uses the momentum equation (St. Venant) with various boundary conditions (subcritical, supercritical, surcharged) to update the conduit’s flow each time step. This file is tightly integrated with the rest of SWMM’s hydraulic engine (found in dynwave.c, flowrout.c, etc.), but here we focus on the finite‐difference solution for one conduit under dynamic wave routing.

---

1. Purpose of dwflow.c

• Dynamic wave flow is SWMM’s most detailed hydraulic mode.
• dwflow.c implements the specialized code to:
  1. Solve the momentum equation in a conduit (pipe or channel) under dynamic wave conditions.
  2. Perform one update of the conduit’s flow, depth, and flow classification at each iteration/time step.
  3. Account for:
     • Friction slope (Manning or Darcy-Weisbach for force mains).
     • Local (entrance/exit) losses.
     • Inertial (acceleration) terms.
     • Preissmann slot approach for surcharged pipes.
     • Evaporation/infiltration losses in the conduit barrel.

The dynamic wave engine calls these routines for each conduit, iterating until a system-wide flow solution converges.

---

2. Key Functions

2.1 dwflow_findConduitFlow(...)

void dwflow_findConduitFlow(int j, int steps, double omega, double dt)

• Inputs:
  • j: the link (conduit) index.
  • steps: how many internal iteration steps have already been taken for the current time step.
  • omega: an under‐relaxation factor (0–1).
  • dt: the current routing time step in seconds.
• Outputs:
  • It updates Link[j].newFlow, Link[j].newDepth, Link[j].newVolume, etc.
• Algorithm:
  1. Gathers upstream/downstream node heads $`h1`, `h2`$ and conduit geometry.
  2. Calculates “effective” flow depths (y1, y2) at each end, ensuring they are $\ge$ FUDGE to avoid zero area.
  3. Finds cross-sectional areas, hydraulic radii, midpoint area, checks if the conduit is full.
  4. Applies the finite difference form of the momentum equation: $$q = \frac{q_{\mathrm{old}} + [\text{inertial, slope, friction, local loss, seepage}]}{\text{denominator}}$$
     • Terms (dq1 through dq6) represent friction slope, gravity slope, inertial difference, local losses, etc.
     • Some inertial damping might apply based on the Froude number or user‐chosen inertial damping option.
  5. Possibly limit the flow to normal flow or inlet‐controlled culvert flow if that is smaller (depending on user’s normal flow or Froude criteria).
  6. Applies under‐relaxation: q_new = (1 - omega)*q_last + omega*q.
  7. Updates various variables: conduit’s mid-area, the new flow depth, volume, flow, etc.

In short, dwflow_findConduitFlow is the core routine that:

• Estimates the new flow in a conduit,
• Balances the momentum equation’s terms,
• Enforces flow or velocity constraints (e.g., normal flow limitation, infiltration, culvert, flap gates).

---

2.2 findSurfArea(...)

static void findSurfArea(int j, double q, double length,
    double* h1, double* h2, double* y1, double* y2)

• Purpose: Divides the conduit’s surface area between its upstream and downstream nodes (for stability in dynamic wave).
• How:
  1. Based on the conduit’s flow classification (subcritical, supercritical, critical, or partially dry), it calculates an “effective top width” and how much portion belongs to each node’s “surface area.”
  2. Adjusts h1, h2, y1, y2 if a node is at critical or normal depth.
  3. Helps the dynamic wave solver by distributing the conduit’s free surface to the node surcharges.

---

2.3 getFlowClass(...)

static int getFlowClass(
    int link, double q, double h1, double h2,
    double y1, double y2, double* yC, double* yN, double* fasnh)

• Classifies the conduit’s flow regime as:
  • SUBCRITICAL, SUPCRITICAL,
  • UP_CRITICAL, DN_CRITICAL,
  • UP_DRY, DN_DRY,
  • or DRY.
• It checks depths and heads at each end, compares them to normal and critical depths (found via link_getYnorm, link_getYcrit).
• fasnh is the fraction between normal & critical depth used to distribute free surface area.

---

2.4 findLocalLosses(...)

static double findLocalLosses(int link, double a1, double a2, double aMid, double q)

• If the user sets local minor loss coefficients (cLossInlet, cLossOutlet, cLossAvg), then for each end or midpoint, we add a term $\sim \text{coefficient} \times q/a$.
• Summed up as part of the momentum equation.

---

2.5 getSlotWidth(...), getWidth(...), getArea(...), getHydRad(...)

• Various helpers to derive:
  • The Preissmann slot top width if surcharged (SLOT method),
  • The top width of the cross‐section at a given depth,
  • The area of flow at a given depth,
  • The hydraulic radius.
• These rely heavily on cross‐section geometry (TXsect) and the presence (or not) of the Preissmann slot.

---

2.6 checkNormalFlow(...)

static double checkNormalFlow(int j, double q, double y1, double y2, double a1, double r1)

• Goal: Possibly override the dynamic wave flow with normal flow if:
  1. The user has set “normal flow limit” (SLOPE, FROUDE, or BOTH).
  2. A slope or Froude condition is met, meaning the flow is “supercritical or steep enough” that normal flow may apply.
• If normal flow is less than the dynamic wave flow, we reduce flow to that normal capacity.
• We do: $$q_{\mathrm{norm}} = \beta \times a_1 \times (r_1)^{2/3}$$ (where $\beta$ is basically $1/n$ or a force main friction factor).
• Then, if qNorm < q, we set q = qNorm.

---

3. Key Data Fields & Variables

• Conduit[k].q1 and Conduit[k].q2: store new flows at the start/end or iteration.
• Conduit[k].barrels: number of parallel barrels. Flow in SWMM is computed per barrel, then multiplied.
• Link[j].oldFlow, Link[j].newFlow: flows in cfs at the end of the last time step vs. current iteration result.
• Link[j].flowClass: the flow classification code (dry, subcritical, supercritical, etc.).
• Link[j].dqdh: partial derivative $\partial Q / \partial H$ used by the dynamic wave solver’s global matrix.
• Link[j].froude: Froude number to detect supercritical conditions.

---

4. Typical Flow of Computation

1. At each iteration in the dynamic wave solver (in dynwave_execute(...)), SWMM calls: \text{dwflow_findConduitFlow}(j, \dots) \quad \text{for each conduit } j.
2. The function:
   1. Reads previous iteration flow (qLast) and node depths (h1, h2).
   2. Determines partial or full flow, calculates cross‐section areas.
   3. Summons the momentum equation in finite difference form: $$q = \frac{q_\mathrm{old} - \text{(gravity slope term)} + \text{(inertial term)} + \dots}{1 + \text{(friction)} + \text{(local losses)}}$$
   4. Possibly modifies q for:
      • Culvert inlet control,
      • Normal flow limits,
      • Flap gates disallowing reverse flow,
      • User flow limit.
   5. Under‐relaxes with omega.
3. The updated flow is stored in Conduit[k].q1, and Link[j].newFlow is q * barrels.
4. Once each conduit is updated, the dynamic wave solver checks for convergence and possibly repeats more iteration passes.

---

5. Key Takeaways

• dwflow.c is specialized code for one conduit’s momentum equation solution under dynamic wave.
• The code ensures we capture:
  • Full or partial pipe, plus the possibility of Preissmann slot for surcharging.
  • Inertial terms that can be damped (none, partial, or full).
  • Local and friction losses.
  • Flow constraints (inlet control, normal flow, flap gates).
• Ultimately, the rest of SWMM orchestrates calls to dwflow_findConduitFlow() within a global loop, each iteration refining flows in all conduits until the entire network converges at a time step—thus giving a robust unsteady flow simulation.

Summary of consts.h in SWMM5

 Below is an overview of consts.h in EPA SWMM5, highlighting what each constant represents and why the file includes them. Essentially, consts.h provides global numeric definitions and common constants that other SWMM modules rely on—making the code more readable, maintainable, and consistent.

---

1. Purpose of consts.h

• Global numeric constants: Instead of scattering “magic numbers” throughout the code, SWMM consolidates them here. That way, a single header file holds units, tolerances, array limits, large/small values, etc.
• Compile-time definitions: By #define’ing these constants, we ensure every SWMM source file that includes consts.h sees the exact same definitions.

---

2. General Constants

#define VERSION 52004

• VERSION: 5.2.004.
  • This integer encodes SWMM version 5.2.4.

#define MAGICNUMBER 516114522

• MAGICNUMBER: A unique integer used in binary or hotstart file headers to confirm the file format.

#define EOFMARK 0x1A

• EOFMARK: Some OS’s or legacy systems used 0x1A (ASCII 26) or 0x04 (ASCII EOT) to mark end-of-file.

#define MAXTITLE 3
#define MAXMSG 1024
#define MAXLINE 1024
#define MAXFNAME 259
#define MAXTOKS 40
#define MAXSTATES 10
#define MAXODES 4

• MAXTITLE: SWMM allows up to 3 title lines in input.
• MAXMSG: Up to 1024 characters in a warning or error message.
• MAXLINE: Up to 1024 characters in any single input line.
• MAXFNAME: Max file name length, 259 characters (on Windows).
• MAXTOKS: Up to 40 tokens (space-delimited) can exist in one input line.
• MAXSTATES: Maximum number of computed hydraulic variables (e.g., flow, depth, velocity, etc.) in routing.
• MAXODES: Max. number of ODE’s SWMM might solve in infiltration or water quality sub-models.

#define NA -1
#define TRUE 1
#define FALSE 0

• NA: Typically “not applicable” or “not available.”
• TRUE, FALSE: C doesn’t have a built-in bool in older standards, so SWMM uses 1/0.

#define BIG 1.E10
#define TINY 1.E-6
#define ZERO 1.E-10
#define MISSING -1.E10

• BIG: A large number (10^10). Often used for initialization or bounding.
• TINY: A small number (10^-6). Used to avoid division by zero or floating underflow.
• ZERO: Even smaller (10^-10).
• MISSING: A placeholder if data is unavailable (-10^10).

#define PI 3.141592654

• PI: Approx. 3.14159…

#define GRAVITY 32.2
#define SI_GRAVITY 9.81

• GRAVITY: Acceleration of gravity (ft/s^2) in US Customary.
• SI_GRAVITY: 9.81 m/s^2 in SI units.

---

3. Manning Equation Factor

#define PHI 1.486

• PHI: The constant 1.486 appears in Manning’s formula for flow in US units: $$Q = 1.486 \cdot A \cdot R^{2/3} / n \;.$$

---

4. Minimum Flow/Depth Constants

#define MIN_RUNOFF_FLOW 0.001
#define MIN_EXCESS_DEPTH 0.0001
#define MIN_TOTAL_DEPTH 0.004167
#define MIN_RUNOFF 2.31481e-8

• MIN_RUNOFF_FLOW = 0.001 cfs: SWMM often zeroes out smaller flows.
• MIN_TOTAL_DEPTH = 0.004167 ft (which is 0.05 inches). Small threshold to differentiate between a dry vs. wet node or conduit.
• **MIN_RUNOFF = 2.31481e-8 ft/s = 0.001 in/hr. Another minimal precipitation-runoff threshold.

---

5. Steady State Tolerances

#define FLOW_TOL 0.00001
#define DEPTH_TOL 0.00001
#define VOLUME_TOL 0.01

• Used to check if changes in flow, depth, or volume are negligible, i.e., “steady state.”

---

6. Minimum Depth for Dynamic Wave

#define FUDGE 0.0001

• FUDGE: A small number (0.0001 ft).
  • If a computed link depth is below FUDGE, SWMM sets it to FUDGE to avoid numeric instability (like dividing by zero area).

---

7. Unit Conversion Factors

#define GPMperCFS 448.831
#define AFDperCFS 1.9837
#define MGDperCFS 0.64632
#define IMGDperCFS 0.5382
#define LPSperCFS 28.317
#define LPMperCFS 1699.0
#define CMHperCFS 101.94
#define CMDperCFS 2446.6
#define MLDperCFS 2.4466
#define M3perFT3 0.028317
#define LperFT3 28.317
#define MperFT 0.3048
#define PSIperFT 0.4333
#define KPAperPSI 6.895
#define KWperHP 0.7457
#define SECperDAY 86400
#define MSECperDAY 8.64e7
#define MMperINCH 25.40

• Provides conversions from cubic feet per second (cfs) to typical volumetric flow units: gallons per minute (GPM), million gallons per day (MGD), etc.
• Also length, pressure, power conversions:
  • MperFT = 0.3048 => 1 ft = 0.3048 m
  • PSIperFT = 0.4333 => psi per foot of water
  • SECperDAY = 86400 => # of seconds in a day
  • MMperINCH = 25.4 => mm per inch

---

8. Token Separators

#define SEPSTR " \t\n\r"

• Characters used as delimiters (space, tab, newline, carriage return) when parsing lines.

---

9. Summary

• consts.h centralizes constants used across SWMM:
  • Array size limits (lines, tokens).
  • Unit conversions (cfs to MGD, ft to m, etc.).
  • Minimum thresholds (flow, depth, volume).
  • Key “magic” numbers (version codes, missing values).

Thanks to these common definitions, the rest of SWMM’s code can cleanly reference ZERO instead of 1.e-10, or MperFT for a length conversion, making the software more modular and maintainable.

Summary of climate.c in SWMM5

 Below is a step‐by‐step overview of how climate.c manages climate‐related data in SWMM5—specifically, temperature, evaporation, wind speed, and monthly adjustment factors. This file handles:

1. Reading climate input parameters (temperature sources, evaporation sources, monthly pattern adjustments).
2. Opening and interpreting external climate data files (various formats).
3. Updating daily/hourly climate states (temperature, evaporation, wind) each simulation day/time step.
4. Applying monthly adjustments (e.g., temperature offset, rainfall multiplier, hydraulic conductivity factor, etc.).

---

1. Main Responsibilities

1. Temperature

   • Allows temperature to come from:
     • A time series in the SWMM input file, or
     • An external climate file, or
     • No temperature data (some default or not used).
   • For min/max daily temperature (from climate file), it models hourly temperature via a sinusoidal interpolation.

2. Evaporation

   • Supports various methods:
     • Constant daily evaporation.
     • Monthly average evaporation for each month of the year.
     • A time series of evaporation rates.
     • File with daily evaporation data.
     • Temperature-based evaporation (Hargreaves method) derived from daily temperature extremes.

3. Wind Speed

   • Can be:
     • Monthly fixed average speeds, or
     • Daily data from the climate file.

4. Adjustments

   • Monthly adjustments for:
     • Temperature, Evaporation, Rainfall, Hyd. conductivity.
   • Patterns for subcatchment infiltration parameters, roughness, etc.

5. File Format Handling

   • Recognizes several climate file formats:
     • User‐prepared (simple ASCII with year, month, day, TMAX, TMIN, EVAP, WIND),
     • NCDC GHCN Daily (GHCND),
     • NCDC TD3200,
     • Canadian DLY02 or DLY04.

---

2. Key Global Variables

Within climate.c, these are prominent:

1. Temp

   • Stores temperature data source, user file date range, current day’s temperature, average daily temp, etc.

2. Evap

   • Type of evaporation method (constant, monthly, timeseries, etc.).
   • Current evaporation rate Evap.rate.

3. Wind

   • Type of wind speed data (monthly or file).
   • Current wind speed Wind.ws.

4. Adjust

   • Holds arrays of 12 monthly adjustments for temp[i], evap[i], rain[i], hydcon[i].
   • Also stores subcatchment adjustment patterns (like Subcatch[i].nPervPattern).

5. File‐specific (like Fclimate)

   • Fclimate.name = external climate file name, Fclimate.file = file pointer.
   • FileYear, FileMonth, FileDay = track which day’s data we are currently on.
   • FileData[var][day] = a month’s worth of daily data for each climate variable.

---

3. Important Functions

3.1 Reading Input Lines

• climate_readParams(...)

  • Reads lines in the [TEMPERATURE] section.
  • Possible tokens:
    • TIMESERIES name
    • FILE name [startDate] [tempUnits]
    • WINDSPEED MONTHLY ...
    • WINDSPEED FILE
    • SNOWMELT v1 ... v6
    • ADC IMPERV/PERV ... (areal depletion curve data)
  • Configures the data source for temperature and wind speeds, plus any snowmelt parameters.

• climate_readEvapParams(...)

  • Reads lines in the [EVAPORATION] section.
  • Possible tokens:
    • CONSTANT value
    • MONTHLY v1 ... v12
    • TIMESERIES name
    • FILE (monthly pan coefficients)
    • TEMPERATURE
    • RECOVERY patternName
    • DRY_ONLY YES/NO

• climate_readAdjustments(...)

  • Reads lines in the [ADJUSTMENTS] section.
  • Could be:
    • TEMPERATURE 12vals
    • EVAPORATION 12vals
    • RAINFALL 12vals
    • CONDUCTIVITY 12vals
    • N-PERV subcatchID patternID
    • DSTORE subcatchID patternID
    • INFIL subcatchID patternID

---

3.2 Validation & Initialization

• climate_validate()

  • Ensures a valid climate file is provided if needed.
  • Opens it (climate_openFile()) if Wind.type==FILE_WIND, Evap.type==FILE_EVAP or Temp.dataSource==FILE_TEMP.
  • Checks snowmelt parameters, latitude bounds, adjusts monthly temperature/evap arrays if UnitSystem == SI.

• climate_openFile()

  • Actually fopen() the climate file and determines file format via getFileFormat().
  • Positions the file at the correct starting date (either the simulation start or Temp.fileStartDate).

• climate_initState()

  • Called at start of simulation.
  • Resets last day, sets up time series evaporation (if used) to track next date & next rate, etc.
  • For temperature-based evaporation, sets up a 7-day moving average structure Tma.

---

3.3 Setting Climate State Each Day/Time

• climate_setState(theDate)
  • Called at each new routing time step to update:
    1. Temperature,
    2. Evaporation,
    3. Wind Speed,
    4. Monthly rainfall/hyd. conductivity multipliers,
    5. NextEvapDate (the next day/time evap changes).
  • Internally calls:
    • updateFileValues(theDate) if using external climate file.
    • setTemp(theDate), setEvap(theDate), setWind(theDate).

---

3.4 Temperature Routines

1. setTemp(theDate)

   • If using a file for TMIN/TMAX, it checks if a new day has started, updates daily TMIN/TMAX, times of sunrise/sunset, etc.
   • For hourly temperature, does sinusoidal interpolation:
     • min temp at sunrise, max at early afternoon, then decreasing.
   • If using a temperature time series, just looks up the temperature in the time series.

2. updateTempTimes(day)

   • Calculates approximate sunrise/sunset hour angles, sets Hrsr, Hrss, etc. to determine when TMIN/TMAX occur.

3. updateTempMoveAve(tmin, tmax)

   • Maintains a 7-day moving average of average daily temperature (tAve) and daily range (tRng).
   • Used by Hargreaves evaporation approach (when Evap.type == TEMPERATURE_EVAP).

---

3.5 Evaporation Routines

1. setEvap(theDate)

   • Chooses evaporation rate based on method:
     • CONSTANT_EVAP: a single value.
     • MONTHLY_EVAP: uses the monthly value for mon-1.
     • TIMESERIES_EVAP: if current date >= NextEvapDate, loads NextEvapRate.
     • FILE_EVAP: uses FileValue[EVAP], possibly scaled by pan coefficients.
     • TEMPERATURE_EVAP: uses the Hargreaves formula with the 7‐day average T.
   • Applies monthly climate adjustments from Adjust.evap[mon-1].

2. getTempEvap(day, tAve, tRng)

   • Implements Hargreaves daily evap formula from 7‐day average T and daily range.
   • Converts from mm/day to in/day if in US units.

3. climate_getNextEvapDate()

   • Returns the current NextEvapDate (the next day/time we might update evaporation rate if from a time series or climate file).

---

3.6 Wind Speed Routines

• setWind(theDate)
  • If monthly, picks the monthly average from Wind.aws[mon-1].
  • If file-based, uses FileValue[WIND].

---

3.7 Climate File Reading Internals

1. climate_openFile() + getFileFormat()

   • Reads the first line.
   • Checks if format is:
     • TD3200 (first 3 chars = "DLY", line[23..26] = "9999"),
     • DLY0204 (Canadian format, line length >= 233, some coded fields),
     • USER_PREPARED (attempts to parse year, month, day),
     • GHCND (checks for "DATE" + possible TMIN/TMAX/Evap lines).
   • Positions the file at the right year/month.

2. readFileValues()

   • Reads a month’s worth of data into FileData[var][day] arrays.
   • For each line in that month, calls specialized parse functions (like parseUserFileLine(), parseTD3200FileLine(), etc.) to fill daily values.

3. updateFileValues(theDate)

   • Each day, we advance the day counters. If we moved to a new month, read a fresh chunk of data from the file.

---

4. Typical Flow (No pun intended)

1. User specifies climate sections in the SWMM input:
   • [TEMPERATURE] lines, [EVAPORATION] lines, [ADJUSTMENTS], etc.
2. During Setup:
   1. climate_validate() ensures it’s consistent, maybe opens the climate file if needed.
   2. climate_initState() initializes internal variables.
3. At Each Simulation Day/Time:
   1. climate_setState(theDate) is called from the runoff or routing engine.
   2. Possibly reads new file data if a new day or new month has begun.
   3. Updates daily TMIN/TMAX or hour interpolation for temperature, sets evaporation, wind speed, etc.
   4. Applies monthly adjustments.

---

5. Key Takeaways

• climate.c is the central hub for SWMM’s climate data. It orchestrates:
  • Reading external climate data in multiple formats,
  • Maintaining internal arrays of daily min/max temperature, evaporation, wind,
  • Interpolating hourly temperature from daily extremes (for rainfall/runoff processes),
  • Calculating Hargreaves evaporation if required,
  • Applying monthly adjustments across temperature, evaporation, rainfall, and infiltration parameters.
• The flexibility (user time series, external file, monthly/constant) allows SWMM to handle many real-world climate scenarios.

Summary of exfil.c in EPA SWMM5

 Below is a step‐by‐step explanation of how exfil.c in EPA SWMM5 processes exfiltration for storage units. The main idea is to handle infiltration (seepage) from the base and sloped banks of a storage node into the surrounding soil—using either a simple constant exfiltration rate or a Green‐Ampt approach.

---

1. Overview: Storage Node Exfiltration in SWMM

• Storage nodes can lose water through their bottom and sidewalls. This is called exfiltration.

• SWMM provides the option for:

  1. Constant infiltration (Ksat) – If you only specify a saturated hydraulic conductivity.
  2. Green–Ampt infiltration – If you provide suction head, saturated hydraulic conductivity, and initial moisture deficit (IMD).

• The Green-Ampt formulation is the same as used for subcatchments, except it’s applied to the bottom and bank areas of a storage node.

---

2. Key Data Structures

1. TExfil

   • Found in exfil.h, this struct holds two TGrnAmpt objects:
     • btmExfil: infiltration parameters for the bottom of the storage unit.
     • bankExfil: infiltration parameters for the banks (sloped side).
   • Also stores btmArea, bankMinDepth, bankMaxDepth, bankMaxArea—the geometry parameters describing the bottom area and side areas.

2. Storage[k].exfil

   • Each storage node k has an optional pointer exfil to a TExfil. If non-NULL, we can compute infiltration through it.

3. TGrnAmpt

   • Data structure for Green-Ampt infiltration parameters (suction head, Ks, IMDmax), plus state variables tracking how much infiltration has occurred so far.

---

3. Function-by-Function Explanation

3.1 exfil_readStorageParams(...)

int exfil_readStorageParams(int k, char* tok[], int ntoks, int n)

• Reads exfiltration parameters from an input line for storage node k.
• Tokens might contain:
  • Just one value (Ksat).
  • Or three values (suction head, Ksat, IMDmax) for Green-Ampt infiltration.
• Steps:
  1. Checks how many tokens remain after n.
  2. If only one token, it interprets it as the saturated conductivity Ksat.
  3. Else, it reads three tokens as suction head, Ksat, and IMDmax.
  4. If Ksat = 0.0, no exfiltration is created (returns 0).
  5. Otherwise, it calls createStorageExfil(...) to actually build a TExfil object and store these parameters.

Why: This parses the user’s input for infiltration/exfiltration parameters and decides if infiltration is constant or if we use a Green-Ampt approach.

---

3.2 exfil_initState(...)

void exfil_initState(int k)

• Called at the start of a simulation (or re-initialization) for a storage node k.

• What it does:

  1. If Storage[k].exfil is not NULL, it calls grnampt_initState on both btmExfil and bankExfil to reset infiltration states (like cumulative infiltration, wetting front depth, etc.).

  2. Determines the geometry of the bottom area and bank area depending on the storage shape:

     • TABULAR shape:

       • Uses the user-defined “Storage Curve” in Curve[i].
       • Looks up area at 0 ft depth to get the bottom area.
       • Then scans the table entries to find the “bank” min depth, max depth, and max area.
       • Converts these from user units to internal units.

     • FUNCTIONAL shape:

       • btmArea = Storage[k].a0 plus optionally a1 if a2 = 0.0.
       • No real max depth or max area, so sets them to BIG.

     • CYLINDRICAL, CONICAL, PYRAMIDAL shapes:

       • btmArea = a0, sets bank depth/area to BIG.

• Why: This sets up the baseline geometry—bottom area used for infiltration, plus the range of depths/areas for the “banks” infiltration.

---

3.3 exfil_getLoss(...)

double exfil_getLoss(TExfil* exfil, double tStep, double depth, double area)

• Main infiltration logic: how much water exfiltrates (cfs) in the current time step.
• Steps:
  1. Bottom infiltration:
     • If IMDmax == 0, we assume a constant infiltration rate = Ksat * Adjust.hydconFactor.
       • Adjust.hydconFactor is a monthly/temperature adjustment for hydraulic conductivity.
     • Else, calls grnampt_getInfil(...) with the bottom infiltration parameters to get the infiltration rate (ft/s).
     • Multiply that infiltration rate by exfil->btmArea to get a flow (cfs).
  2. Bank infiltration (for side slopes):
     • Only if water depth > exfil->bankMinDepth.
     • The bank area is min(area, exfil->bankMaxArea) - btmArea.
       • i.e., total water surface area minus the bottom area, but capped at bankMaxArea.
     • If IMDmax == 0, infiltration is just area × Ksat.
     • Else, we approximate the average depth above the bank using either:
       • If actual depth > bankMaxDepth, we do depth = depth - bankMaxDepth + (bankMaxDepth - bankMinDepth)/2.
       • If below top of bank, depth = (depth - bankMinDepth)/2.
     • Then calls grnampt_getInfil(...) again with bankExfil parameters and multiplies by the side bank area.
  3. Returns the total exfiltration rate (cfs).

Why: A storage node can have infiltration from both the bottom and the side banks, and the code uses Green-Ampt or constant infiltration to estimate these losses.

---

3.4 createStorageExfil(...)

int createStorageExfil(int k, double x[])

• Builds a new TExfil object for storage node k.
• Steps:
  1. Allocates memory for TExfil if not already created.
  2. Allocates two TGrnAmpt objects: one for bottom infiltration, one for bank infiltration.
  3. Calls grnampt_setParams(...) on both with the array x[] (suction head, Ksat, IMDmax).
  4. If memory allocation fails, sets error code.
• Why: This is how SWMM sets up the infiltration (exfiltration) mechanism for a particular storage node.

---

4. Putting It All Together

1. Reading Input: The user’s input line for a storage node might say “Green-Ampt parameters” or just a single Ksat. SWMM calls exfil_readStorageParams(...), which might create a TExfil if Ksat > 0.
2. Initialization: During exfil_initState(...), we determine the bottom area and potential bank area for infiltration, set initial infiltration states to zero, etc.
3. Every Time Step:
   • For each storage node with an exfil object, SWMM calls exfil_getLoss(...) to compute infiltration outflow.
   • That infiltration rate is subtracted from the water volume in the node, thus lowering the node’s depth over time.

The net effect is that storage nodes in SWMM can slowly lose water into the ground, either with a constant infiltration or a Green–Ampt infiltration approach that depends on soil suction, conductivity, and moisture deficit—applied separately to bottom and bank areas.

---

5. Key Highlights & Design Choices

1. Two Infiltration Regions

   • Bottom: Horizontal infiltration area.
   • Banks: Sloped infiltration area, which might vary with water depth.

2. Green‐Ampt vs. Constant

   • If IMDmax == 0.0, infiltration is simply constant Ksat × area.
   • Otherwise, the modified Green–Ampt approach is used.

3. Monthly/Temperature Adjustments

   • Adjust.hydconFactor is a multiplier that can be changed monthly in SWMM’s advanced settings (e.g., water temperature effect on viscosity/conductivity).

4. Shapes

   • For TABULAR storage shapes, the curve can show how cross-sectional area increases with depth—this code identifies the “bottom area” at depth = 0 and the portion of the surface area that could be considered “bank.”

5. Performance

   • Everything is done each time step for each storage node. This approach is straightforward but must be carefully handled if many storage nodes exist.

Overall, exfil.c is the part of SWMM that converts infiltration input parameters into actual water losses from storage nodes, integrating a Green‐Ampt infiltration model (or a simpler Ksat-based approach) into the dynamic routing of stored water.

SWMM5 flowrout.c Summary

 Below is a step‐by‐step explanation of how flowrout.c implements and coordinates flow routing in SWMM5. The functions here support Steady Flow, Kinematic Wave, and (indirectly) Dynamic Wave routing methods. They initialize link/node states, compute flow propagation each time step, and update depths and volumes at nodes and links.

---

1. Role of flowrout.c

• In SWMM5, flow routing is the process of moving water through the network of nodes and links based on continuity and momentum principles.
• flowrout.c is a high‐level module that:
  1. Initializes node/link depths and volumes when the simulation starts.
  2. Chooses the correct solver approach (steady, kinematic, or dynamic wave).
  3. Executes a routing step, updating flow rates, depths, and volumes at each node/link for the current time step.

Note that Dynamic Wave routing has much of its logic in dynwave.c, but flowrout.c still coordinates calls to it.

---

2. Key Externally Called Functions

1. flowrout_init(routingModel)

   • Called when the routing process is opened (one time at start).
   • If Dynamic Wave:
     1. Calls validateGeneralLayout() to ensure the network is topologically consistent for dynamic wave.
     2. Calls dynwave_init() to set up dynamic wave parameters.
     3. If no hotstart file is used, initializes node and link depths with initNodeDepths() and initLinkDepths().
   • If Kinematic Wave or Steady Flow:
     • Calls validateTreeLayout() which checks for a “tree-like” network structure: each node has at most one (or two) outlets, no adverse slopes unless dummy links, etc.
   • Then calls initNodes() and initLinks(routingModel) to set initial volumes/flows.

2. flowrout_close(routingModel)

   • Called at the end of the simulation.
   • If dynamic wave was used, calls dynwave_close() to free memory used by the dynamic wave engine.

3. flowrout_getRoutingStep(routingModel, fixedStep)

   • Returns the routing time step to use for the next iteration.
   • If Dynamic Wave, calls dynwave_getRoutingStep() to find a variable time step that meets Courant stability criteria. Otherwise, returns fixedStep.

4. flowrout_execute(links, routingModel, tStep)

   • Main function that routes flow through the network for a single time step (tStep).
   • Steps:
     1. Resets node overflow to 0.
     2. If routingModel == DW, calls dynwave_execute(tStep) and returns.
     3. Else (Steady or Kinematic wave), processes links in topological order:
        • If the upstream node is a storage, calls updateStorageState() to iteratively find the node’s outflow and final water level.
        • Gets inflow to the link (getLinkInflow()) and routes flow either with:
          • steadyflow_execute() or
          • kinwave_execute()
        • The result is assigned to Link[j].newFlow.
     4. After going through all links, updates each node’s and link’s final state (setNewNodeState() and setNewLinkState()).
   • Returns the number of computational sub-steps performed (e.g., for partial iterative updates in kinematic or steady flow).

---

3. Important Internal Functions

3.1 Initialization & Validation

1. validateTreeLayout()

   • Ensures the network is “tree-like” for Kinematic Wave or Steady Flow.
   • Checks each node’s outlet links (dividers, outfalls, storage) and ensures no node has multiple downstream links (unless it’s a valid exception).

2. validateGeneralLayout()

   • Used by Dynamic Wave. Checks for conditions like dummy links, ideal pumps, or outfalls that must have limited connectivity. Ensures at least one outlet node, etc.

3. initNodeDepths()

   • Assigns initial water depth to nodes if not specified.
   • E.g., for a non-storage node with no user-supplied depth, it might average the depths of its connecting links.

4. initLinkDepths()

   • For Dynamic Wave only—initial depth in each conduit is the average of the depths at its end nodes, unless a user-supplied q0 is given.

5. initNodes() & initLinks()

   • initNodes(): sets initial inflow/outflow, volume, overflow, etc.
   • initLinks(): sets initial link volumes (or zero flow for Steady Flow) and prepares internal states in conduits/pumps.

3.2 Steady/Kinematic Wave Execution

1. flowrout_execute(...) (as explained above).

2. updateStorageState(i, j, links, dt)

   • For a storage node i, we do an iterative approach to find the node’s new depth.
   • Because outflow from a storage node depends on the water depth in it, we do multiple iterations up to MAXITER (10) or until the depth change is under STOPTOL (0.005 ft).
   • We call getStorageOutflow(...) to sum flows leaving the node’s connected links.

3. getLinkInflow(link, dt)

   • Finds how much flow can enter a link from its upstream node.
   • For example, we might limit it if the node cannot supply more than node_getMaxOutflow(...).

4. steadyflow_execute(...)

   • A simplified approach for Steady Flow:
     1. If it’s a conduit, tries to see if the input flow exceeds full capacity; if so, it limits flow to full capacity.
     2. Adjusts for evaporation/infiltration losses via link_getLossRate(...).

5. kinwave_execute(...)

   • (Not shown in detail in this code snippet.) Similar approach but uses kinematic wave formulas.

3.3 Updating Final States

• After each link is routed, we do:
  • setNewNodeState(...):
    • For non‐storage nodes, updates volume = oldVolume + netInflow × dt.
    • If volume > fullVolume and ponding not allowed, set overflow and cap volume.
    • Depth is then node_getDepth(j, volume).
  • setNewLinkState(...):
    • For a conduit, we find average area from upstream and downstream cross-sectional areas, compute link volume, and set new depth based on y1 + offset, y2 + offset, etc.

---

4. Putting It All Together in a Typical Time Step

1. SWMM calls flowrout_getRoutingStep() to find tStep.
2. flowrout_execute(links, model, tStep) is invoked:
   1. If Dynamic Wave: goes directly to dynwave_execute(tStep), which solves St. Venant PDEs iteratively.
   2. Else (Steady or Kinematic):
      • Loops through each link in a topologically sorted order.
      • If the link’s upstream node is a storage, calls updateStorageState(...) to iteratively refine the storage node’s depth.
      • Finds link inflow, routes the flow with steady or kinematic formula.
      • Updates link’s newFlow, increments node outflow/inflow.
      • After all links are processed, calls setNewNodeState(...) for each node and setNewLinkState(...) for each link to finalize volumes and depths.
3. Return the count of computational steps used (1 step per link in steady/kinwave, or the iterative steps in dynamic wave).

---

5. Key Constants & Notes

• OMEGA = 0.55: Under‐relaxation factor for storage iteration—blends old depth with newly calculated depth to stabilize.
• MAXITER = 10: Max iteration passes for updating storage nodes.
• STOPTOL = 0.005 ft: Convergence tolerance on storage node depth.
• AllowPonding: If set, a node can have indefinite capacity beyond fullDepth by storing water in a “ponded area.”

---

6. Conclusions

• flowrout.c is the routing orchestrator. It sets up initial states (volumes/depths), picks the correct solver path, and after each step updates node/link states.
• For Steady and Kinematic wave, it uses simplified, topological progression: from upstream to downstream links, adjusting flows and volumes in each storage node.
• For Dynamic Wave, it delegates to dynwave_execute(), but still handles some setup/close routines.
• The iterative logic in updateStorageState() is a key mechanism for bridging the fact that outflow from a storage node depends on the node’s depth, and the node’s depth depends on how much outflow it has—hence the repeated approximation each time step.

SWMM5 Forcemain.c Summary

Below is an explanation of what each function in forcemain.c does, why it’s there, and how it fits into SWMM’s hydraulic calculations for pressurized force mains. This file handles Hazen-Williams and Darcy-Weisbach friction equations under force main conditions—i.e., when pipes are pressurized and do not follow standard Manning’s equation.

---

1. Background: Force Mains in SWMM

• In typical open-channel flow, SWMM uses Manning’s equation to model headloss.

• However, a force main is a pressurized conduit (often from a pump), where flow is not necessarily open-channel. In that scenario, many engineers use:

  • Hazen-Williams equation (common in North America for water distribution), or
  • Darcy-Weisbach equation (a more universal approach, used in many hydraulic contexts).

• SWMM’s user can choose either Hazen-Williams (H_W) or Darcy-Weisbach (D_W) for force mains. The code in forcemain.c provides the specialized calculations for friction losses, Reynolds number, friction factors, etc.

---

2. The Main Functions

2.1 forcemain_getEquivN(...)

double forcemain_getEquivN(int j, int k)

Purpose: Computes an “equivalent Manning’s n” for a force main flowing full under fully turbulent conditions, using the user’s chosen equation (Hazen-Williams or Darcy-Weisbach).

1. Inputs:
   • j = link index in SWMM’s data structures
   • k = conduit index
   • The conduit’s geometry is accessed via Link[j].xsect and Conduit[k].
2. Logic:
   • If Hazen-Williams:

     return 1.067 / xsect.rBot * pow(d/Conduit[k].slope, 0.04);
     

     • d is yFull (full diameter or full depth).
     • xsect.rBot here is repurposed to store the Hazen-Williams “C” factor.
     • The constant 1.067 and exponent 0.04 come from deriving an approximate equivalence between Hazen-Williams and Manning’s formula when the pipe is completely full.
   • If Darcy-Weisbach:

     f = forcemain_getFricFactor(xsect.rBot, d/4.0, 1.0e12);
     return sqrt(f/185.0) * pow(d, (1./6.));
     

     • xsect.rBot is used here to store the roughness height (e) for Darcy-Weisbach.
     • It calls forcemain_getFricFactor(...) with re = 1.0e12—effectively a giant Reynolds number—to ensure fully turbulent regime.
     • Then it converts that friction factor into an “equivalent Manning’s n.”
   • If neither (i.e., normal Manning’s is used), it just returns the conduit’s normal roughness.
3. Why do this?
   • Internally, SWMM might do normal flow checks or partial flow logic that sometimes calls on Manning’s equation. This function helps “fake” a Manning’s n that yields the same full-flow result as Hazen-Williams or Darcy-Weisbach.

---

2.2 forcemain_getRoughFactor(...)

double forcemain_getRoughFactor(int j, double lengthFactor)

Purpose: Computes an adjustment factor for the pipe’s roughness that compensates for any artificial lengthening the user or the model might impose. Sometimes, a pipe is artificially lengthened to maintain stability in the SWMM solver (i.e., if a pipe is too short, you might virtually lengthen it but keep the slope the same).

1. Inputs:
   • j = link index
   • lengthFactor = factor by which a pipe is “artificially lengthened.” (If the pipe length is doubled artificially, lengthFactor might be 2.0.)
2. Logic:
   • If Hazen-Williams (H_W):

     r = 1.318*xsect.rBot*pow(lengthFactor, 0.54);
     return GRAVITY / pow(r, 1.852);
     

     • xsect.rBot is the Hazen-Williams C-factor again.
     • The 1.318 and exponents 0.54, 1.852 come from the Hazen-Williams formula.
   • If Darcy-Weisbach (D_W):

     return 1.0/8.0/lengthFactor;
     

     • This modifies the friction factor for the artificial length.
3. Why do this?
   • If you artificially lengthen a pipe, you need to adjust the friction or roughness so that the final headloss is physically consistent with the original “real” pipe.

---

2.3 forcemain_getFricSlope(...)

double forcemain_getFricSlope(int j, double v, double hrad)

Purpose: Computes the headloss per unit length (i.e., friction slope $S_f$) in the dynamic wave flow routing for a pressurized force main. It uses either Hazen-Williams or Darcy-Weisbach, depending on ForceMainEqn.

1. Inputs:
   • j = link index
   • v = flow velocity (ft/s)
   • hrad = hydraulic radius (ft) (for a full pipe, hrad = diameter/4, but if partially pressurized, might differ).
2. Logic:
   • If Hazen-Williams:

     return xsect.sBot * pow(v, 0.852) / pow(hrad, 1.1667);
     

     • The code previously stored a “roughness factor” in xsect.sBot (during conduit validation), so that it can just multiply by v^0.852 / hrad^1.1667.
   • If Darcy-Weisbach:

     re = forcemain_getReynolds(v, hrad);
     f = forcemain_getFricFactor(xsect.rBot, hrad, re);
     return f * xsect.sBot * v / hrad;
     

     • First calculates Reynolds number.
     • Then obtains the friction factor f.
     • Then friction slope = f * sBot * (v / hrad).
3. Why do this?
   • SWMM’s dynamic wave solver needs friction slope (or friction headloss per length) at each time step. This is the crux of how SWMM calculates energy losses in pressure flow.

---

2.4 forcemain_getReynolds(...)

double forcemain_getReynolds(double v, double hrad)

Purpose: Computes the flow’s Reynolds number:

$$\text{Re} = \frac{4 \times \text{(hydraulic radius)} \times v}{\nu}$$

Here, $\nu = \text{VISCOS} = 1.1 \times 10^{-5} \,\text{ft}^2/\text{s}$ is the kinematic viscosity of water at 20°C.

1. Formula:

   return 4.0 * hrad * v / VISCOS;
   

2. Why: Reynolds number tells us if the flow is laminar, transitional, or turbulent, which the code then uses to pick the right friction factor formula.

---

2.5 forcemain_getFricFactor(...)

double forcemain_getFricFactor(double e, double hrad, double re)

Purpose: Computes the Darcy-Weisbach friction factor $f$ for a force main using the Swamee-Jain approximation to the Colebrook-White equation.

1. Inputs:
   • e: the roughness height (ft). In SWMM, this might be set by the user or derived from the “rBot” for Darcy-Weisbach.
   • hrad: hydraulic radius.
   • re: Reynolds number.
2. Logic:
   • Check for extremely low Re: if re < 10.0, set re = 10.0 artificially.
   • Laminar regime (Re <= 2000):

     f = 64.0 / re;
     

   • Transitional (2000 < Re < 4000):

     // recursively calls forcemain_getFricFactor with re=4000, then
     // linearly blend from f=0.032 up to that friction factor
     

     Actually, the code calls itself for Re=4000, gets that friction factor, and does a linear interpolation.
   • Turbulent (Re >= 4000):

     f = e/3.7/(4.0*hrad);
     if (re < 1.0e10) f += 5.74/pow(re, 0.9);
     f = log10(f);
     f = 0.25 / (f*f);
     

     • This is the Swamee-Jain formula: $$f = 0.25 \Big/ \Big[\log_{10}\big(\frac{e/(3.7 D)} + \frac{5.74}{Re^{0.9}}\big)\Big]^2$$ but here D is approximated as 4 * hrad if it’s a circular pipe.
3. Why:
   • The friction factor is crucial for Darcy-Weisbach calculations. The Swamee-Jain approximation is a direct (non-iterative) way to estimate friction factor from Colebrook-White.

---

3. Putting It All Together

1. Selecting Equation:

   • SWMM’s global setting ForceMainEqn is either H_W (Hazen-Williams) or D_W (Darcy-Weisbach).
   • Each force main (pressurized conduit) is flagged to use one or the other.

2. During Simulation:

   1. Initialization:
      • SWMM calls forcemain_getEquivN() at validation time to store a stand-in Manning’s n for certain normal-flow checks.
      • It also calls forcemain_getRoughFactor(...) to store the pipe’s friction coefficient in xsect.sBot.
   2. Every Time Step (Dynamic Wave):
      • SWMM calculates velocity v, hydraulic radius hrad, etc.
      • Then calls forcemain_getFricSlope(...) to get friction slope, which goes into the momentum equation.
      • Inside that function, if Darcy-Weisbach is chosen, the code calls:
        • forcemain_getReynolds(...) to get Re,
        • forcemain_getFricFactor(...) for Darcy-Weisbach’s friction factor.
      • This friction slope is used in the St. Venant equations for unsteady flow.

3. Result:

   • The code ensures force main headloss is modeled with either H-W or D-W.
   • The user sees more accurate pressurized-flow behavior than if SWMM were just forcing a Manning’s formula.

---

4. Key Takeaways

1. Hazen-Williams (H_W) vs. Darcy-Weisbach (D_W)

   • H_W: simpler, widely used in water distribution, but less general for different fluids/temperatures.
   • D_W: more universal, uses friction factor $f$ derived from roughness and Reynolds number (Colebrook-White/Swamee-Jain).

2. Equivalent Manning’s n

   • Because SWMM’s internal logic often uses Manning-based checks, the code calculates a “fake n” that yields the same full-flow rate as the chosen force main formula under ideal conditions.

3. Friction Factor

   • The code includes transitional laminar/turbulent handling.
   • Swamee-Jain is a closed-form approximation to Colebrook-White, avoiding iterative solutions for friction factor.

4. Artificial Pipe Lengthening

   • The user can artificially lengthen a short pipe for numerical stability. The code corrects friction to preserve real-world headloss.

Overall, forcemain.c is the specialized module that overrides normal Manning’s friction for pressurized pipes in SWMM, using either Hazen-Williams or Darcy-Weisbach—complete with Reynolds-number-based friction factors—so that force mains in your model behave more realistically than if they used open-channel Manning’s equations.

SWMM5 Controls.c Code

 Below is a walkthrough of controls.c from EPA SWMM5. This file handles rule-based controls, allowing a user to write logic rules (in the model input file) that dynamically change pump, gate, or orifice settings during a simulation. The code parses control rules, stores them, and executes them at each simulation time step.

---

1. Purpose of the File

• controls.c implements SWMM’s “Advanced Control Rules.”
• These rules are text-based, with an IF <condition> THEN <action> ELSE <action> format plus an optional PRIORITY.
• The code:
  1. Parses the rules from the input file into data structures.
  2. Evaluates those rules every time step to see which actions apply.
  3. Executes actions (e.g., change a pump from ON to OFF or set an orifice opening).

---

2. Overall Flow

1. During input parsing, each line that starts with RULE, IF, THEN, ELSE, PRIORITY, VARIABLE, or EXPRESSION is handled by functions in this file.
2. Variables (or expressions) are defined that can be referenced in rule premises.
3. Premises are stored in a linked list of conditions (e.g., “Node 123 Depth > 4”).
4. Actions (e.g., “Pump XYZ status = OFF”) are stored similarly, but separate from premises.
5. Each simulation time step:
   • controls_evaluate() goes through every rule, checks the premises, determines if THEN or ELSE actions should fire, and records which actions to apply.
   • The highest‐priority action for any single link overrides any lower‐priority actions for that link in that time step.
   • These actions then modify the link’s targetSetting (the new setting for a pump, orifice, etc.).

---

3. Key Data Structures

3.1 TRule (one per rule)

struct TRule
{
   char*    ID;             // rule's name/ID
   double   priority;       // rule's priority
   TPremise* firstPremise;  // pointer to the first premise (condition)
   TPremise* lastPremise;   // pointer to the last premise
   TAction*  thenActions;   // linked list of actions if premises are true
   TAction*  elseActions;   // linked list of actions if premises are false
};

• A single rule can have multiple “premise” lines (joined by AND/OR) and multiple “actions” in the THEN or ELSE block.

3.2 TPremise (linked list of conditions within each rule)

struct TPremise
{
    int     type;       // r_IF, r_AND, or r_OR
    int     exprIndex;  // if premise references a math expression
    TVariable lhsVar;   // left-hand side variable
    TVariable rhsVar;   // right-hand side variable
    int     relation;   // relational operator (>, <, =, etc.)
    double  value;      // numerical value on the right-hand side (if not a variable)
    TPremise *next;     // pointer to next premise in this rule
};

• Example premise: Node 10 Depth > 5.0
  • lhsVar = (object = NODE, index=10, attribute=DEPTH)
  • relation = GT (>)
  • value = 5.0

3.3 TAction (linked list of actions in the THEN/ELSE blocks)

struct TAction
{
   int     rule;       // index of the rule
   int     link;       // link index being controlled
   int     attribute;  // e.g. STATUS, SETTING, etc.
   int     curve;      // reference to a curve if modulated by a curve
   int     tseries;    // reference to a time series if modulated
   double  value;      // direct numeric setting
   double  kp, ki, kd; // PID coefficients if using a PID controller
   double  e1, e2;     // stored PID errors
   TAction* next;      // next action in THEN or ELSE
};

• Example action: Pump P1 status = OFF
  • link = index of link P1
  • attribute = STATUS
  • value = 0 (OFF)

3.4 TNamedVariable and TExpression

struct TNamedVariable
{
    TVariable variable;         
    char      name[MAXVARNAME]; // the variable's name used in rules
};

struct TExpression
{
    MathExpr* expression;       
    char       name[MAXVARNAME];
};

• A named variable might be declared in the input file as VARIABLE MyVar = Node 123 Depth.
• An expression might be declared as EXPRESSION MyExpr = MyVar + 1.5.

3.5 Global Arrays

• Rules: an array of size RuleCount.
• NamedVariable: array of named variables.
• Expression: array of parsed math expressions.

---

4. Main Functions

4.1 Initialization & Cleanup

1. controls_init()

   • Sets up defaults (Rule array pointers = NULL, counters = 0, etc.).

2. controls_addToCount(char* s)

   • If a line in the input is VARIABLE ..., increment the count of variables to be stored.
   • If a line is EXPRESSION ..., increment the count of expressions.

3. controls_create(int n)

   • Allocates memory for the global Rules array of size n.
   • Also allocates arrays for any named variables or expressions.

4. controls_delete()

   • Frees memory used by rules, premises, actions, expressions, etc.

4.2 Adding Variables & Expressions

• controls_addVariable(...)

  • Input format: VARIABLE varName = Node ID Depth (or SIMULATION TIME, etc.)
  • Stores a new TNamedVariable in NamedVariable[] so it can be referenced in rule premises or expressions.

• controls_addExpression(...)

  • Input format: EXPRESSION exprName = <math expression>
  • Uses mathexpr_create(...) to parse the expression into tokens.
  • This expression can be used later in premises, for example: IF MyExpr > 6.0 THEN ...

4.3 Creating & Parsing Rules

• controls_addRuleClause(r, keyword, tok[], nToks)
  • Called once a rule is encountered.
  • Dispatches to addPremise(...) if the clause is IF, AND, OR, or to addAction(...) if THEN, ELSE.
  • Also handles PRIORITY.

4.3.1 Premises

• addPremise(r, type, tok, nToks)
  • Creates a new TPremise (e.g., Node 10 Depth > 5.0).
  • If the premise’s right‐hand side is another variable, it sets that in rhsVar; otherwise, it stores a numeric value.

4.3.2 Actions

• addAction(r, tok, nToks)
  • Creates a TAction, e.g., “Pump P1 status = OFF.”
  • Identifies if the action is a direct numeric setting, or references a “Curve,” “TimeSeries,” or “PID.”

---

5. Rule Evaluation

5.1 controls_evaluate(currentTime, elapsedTime, tStep)

1. Sets some global time info (CurrentDate, CurrentTime, etc.).

2. Resets an ActionList (a temporary list of actions triggered this time step).

3. Loops over every rule r:

   1. Evaluate premises (TRUE or FALSE).
      • Each premise can be AND or OR type:
        • For AND, all premises must be TRUE to satisfy.
        • For OR, if any premise is true, the entire premise chain is satisfied.
   2. If premises are TRUE, pick the THEN actions; otherwise the ELSE actions.
   3. For each action:
      • Possibly compute a modulated value from a curve/time series or apply a PID loop.
      • Add the action to the ActionList with the rule’s priority.

4. After all rules have been processed, executeActionList(...) is called:

   • Goes through the ActionList.
   • If multiple rules propose conflicting actions for the same link, only the highest‐priority action is used.
   • The winning action’s value is written to the link’s targetSetting.

5.2 Premise Evaluation

• evaluatePremise(TPremise* p, double tStep)
  • Retrieves the lhsValue (left-hand side): either from a variable or an expression.
  • Retrieves the rhsValue: from a numeric constant or from a variable.
  • If it’s a time-based attribute (TIME, CLOCKTIME, etc.), uses compareTimes(...) with a half time-step tolerance so that an event of TIME = 10:00 is true if the time is within ± (tStep/2).
  • Otherwise uses compareValues(lhsValue, relation, rhsValue) (e.g., lhsValue > rhsValue).
  • The function returns TRUE or FALSE.

---

6. Special Features

1. Named Variables

   • Declared as VARIABLE MyVar = .... Allows referencing a node depth or flow by a simpler name.

2. Expressions

   • Declared as EXPRESSION MyExpr = MyVar + 3.0 (or any math formula).
   • The code uses mathexpr_create(...) and mathexpr_eval(...) to parse and evaluate these on the fly during premise evaluation.

3. PID Control

   • A rule action can be Pump P1 setting = PID Kp Ki Kd, providing continuous feedback control.
   • The code in getPIDSetting(...) implements the standard discrete‐time PID formula, adjusting the link’s setting each time step based on the error between a “SetPoint” and the measured “ControlValue.”

4. Time/Date Attributes

   • A premise can check SIMULATION TIME, DATE, CLOCKTIME, DAYOFYEAR, etc.
   • The code compares times with a half-step offset to avoid missing an exact match.

---

7. Action Execution & Priority

7.1 ActionList

• A temporary singly linked list of type TActionList that collects all actions triggered in the current time step.
• Each node points to a TAction.

7.2 updateActionList(a)

• If the same link is already on the list, compare rule priorities.
• The highest priority action wins for that link.

7.3 executeActionList(currentTime)

• Actually sets Link[a->link].targetSetting = a->value for each action.
• Counts how many new settings changed from the previous time step.
• Optionally writes a log of the action to the report file if RptFlags.controls is enabled.

---

8. Cleanup Routines

• deleteActionList(): frees the ActionList linked list.
• deleteRules(): for each rule, frees all premise nodes (TPremise) and action nodes (TAction). Finally frees the Rules array.

---

9. Key Takeaways

1. Rule Grammar

   • RULE name
   • IF/AND/OR <premise> lines
   • THEN/ELSE <action> lines
   • PRIORITY <p>

2. Premise

   • Compares attributes (depth, inflow, date, time, etc.) of nodes/links/gages to either numeric values or other attributes.
   • You can also use named variables or math expressions.

3. Action

   • Sets the status or setting of a link object (pump, orifice, weir, or conduit).
   • Could be a constant, or modulated by a curve/time series, or governed by PID.

4. Execution

   • Each time step, all rules are evaluated.
   • The highest‐priority action for each link takes effect.
   • The link’s targetSetting is updated accordingly, which the hydraulic solver then uses.

5. Integration with SWMM

   • Rules do not directly override the actual computed flow in a link; they only set the link’s target status or setting.
   • The hydraulic solver then uses that status/setting as a boundary condition.

In short, controls.c is the logic engine for SWMM’s Rule-Based Controls, handling everything from parsing the rule statements to actually applying them to the hydraulic model at each time step.