SWMM5 link.c Summary

 Below is the consolidated latest version of link.c as of SWMM 5.2.4 (2023-06-12), combining all incremental updates, bug fixes, and new features documented in the update history. These updates include (but are not limited to) support for surcharging weirs, monthly infiltration adjustments in conduits, support for Streets and Inlets, variable speed pumps, limiting total conduit losses by flow or conduit volume, and additional warning messages for conduits and regulators with adverse conditions.

//-----------------------------------------------------------------------------
//   link.c
//
//   Project:  EPA SWMM5
//   Version:  5.2
//   Date:     06/12/23   (Build 5.2.4)
//   Author:   L. Rossman
//             M. Tryby (EPA)
//
//   Conveyance system link functions
//
//   Update History
//   ==============
//   Build 5.1.007:
//   - Optional surcharging of weirs introduced.
//   Build 5.1.008:
//   - Bug in finding flow through surcharged weir fixed.
//   - Bug in finding if conduit is upstrm/dnstrm full fixed.
//   - Monthly conductivity adjustment applied to conduit seepage.
//   - Conduit seepage limited by conduit's flow rate.
//   Build 5.1.010:
//   - Support added for new ROADWAY_WEIR object.
//   - Time of last setting change initialized for links.
//   Build 5.1.011:
//   - Crest elevation of regulator links raised to downstream invert.
//   - Fixed converting roadWidth weir parameter to internal units.
//   - Weir shape parameter deprecated.
//   - Extra geometric parameters ignored for non-conduit open rectangular
//     cross sections.
//   Build 5.1.012:
//   - Conduit seepage rate now based on flow width, not wetted perimeter.
//   - Formula for side flow weir corrected.
//   - Crest length contraction adjustments corrected.
//   Build 5.1.013:
//   - Maximum depth adjustments made for storage units that can surcharge.
//   - Support added for head-dependent weir coefficient curves.
//   - Adjustment of regulator link crest offset to match downstream node invert
//     now only done for Dynamic Wave flow routing.
//   Build 5.1.014:
//   - Conduit evap. and seepage losses initialized to 0 in conduit_initState()
//     and not allowed to exceed current flow rate in conduit_getLossRate().
//   Build 5.2.0:
//   - Support added for Streets and Inlets.
//   - Support added for variable speed pumps.
//   Build 5.2.1:
//   - Warning no longer issued when conduit elevation drop < MIN_DELTA_Z.
//   Build 5.2.2:
//   - Warning for conduit elevation drop < MIN_DELTA_Z restored.
//   Build 5.2.4:
//   - Conduit evap+seepage loss under DW routing limited by conduit volume.
//-----------------------------------------------------------------------------

#define _CRT_SECURE_NO_DEPRECATE

#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "headers.h"
#include "inlet.h"

//-----------------------------------------------------------------------------
//  Constants
//-----------------------------------------------------------------------------
static const double MIN_DELTA_Z = 0.001; // minimum elevation change for conduit
                                         // slopes (ft)

//-----------------------------------------------------------------------------
//  External functions (declared in funcs.h)
//-----------------------------------------------------------------------------
//  link_readParams        (called by parseLine in input.c)
//  link_readXsectParams   (called by parseLine in input.c)
//  link_readLossParams    (called by parseLine in input.c)
//  link_validate          (called by project_validate in project.c)
//  link_initState         (called by initObjects in swmm5.c)
//  link_setOldHydState    (called by routing_execute in routing.c)
//  link_setOldQualState   (called by routing_execute in routing.c)
//  link_setTargetSetting  (called by routing_execute in routing.c)
//  link_setSetting        (called by routing_execute in routing.c)
//  link_getResults        (called by output_saveLinkResults)
//  link_getLength         (called in dwflow.c, kinwave.c & flowrout.c)
//  link_getFroude         (called in dwflow.c)
//  link_getInflow         (called in flowrout.c & dynwave.c)
//  link_setOutfallDepth   (called in flowrout.c & dynwave.c)
//  link_getYcrit          (called by link_setOutfallDepth & in dwflow.c)
//  link_getYnorm          (called by conduit_initState, link_setOutfallDepth
//                          & in dwflow.c)
//  link_getVelocity       (called by link_getResults & stats_updateLinkStats)
//  link_getPower          (called by stats_updateLinkStats in stats.c)
//  link_getLossRate       (called in dwflow.c, kinwave.c & flowrout.c)

//-----------------------------------------------------------------------------
//  Local functions
//-----------------------------------------------------------------------------
static void   link_setParams(int j, int type, int n1, int n2, int k, double x[]);
static void   link_convertOffsets(int j);
static double link_getOffsetHeight(int j, double offset, double elev);

static int    conduit_readParams(int j, int k, char* tok[], int ntoks);
static void   conduit_validate(int j, int k);
static void   conduit_initState(int j, int k);
static void   conduit_reverse(int j, int k);
static double conduit_getLength(int j);
static double conduit_getLengthFactor(int j, int k, double roughness);
static double conduit_getSlope(int j);
static double conduit_getInflow(int j);
static double conduit_getLossRate(int j, int routeModel, double q, double tstep);

static int    pump_readParams(int j, int k, char* tok[], int ntoks);
static void   pump_validate(int j, int k);
static void   pump_initState(int j, int k);
static double pump_getInflow(int j);

static int    orifice_readParams(int j, int k, char* tok[], int ntoks);
static void   orifice_validate(int j, int k);
static void   orifice_setSetting(int j, double tstep);
static double orifice_getWeirCoeff(int j, int k, double h);
static double orifice_getInflow(int j);
static double orifice_getFlow(int j, int k, double head, double f,
              int hasFlapGate);

static int    weir_readParams(int j, int k, char* tok[], int ntoks);
static void   weir_validate(int j, int k);
static void   weir_setSetting(int j);
static double weir_getInflow(int j);
static double weir_getOpenArea(int j, double y);
static void   weir_getFlow(int j, int k, double head, double dir,
              int hasFlapGate, double* q1, double* q2);
static double weir_getOrificeFlow(int j, double head, double y, double cOrif);
static double weir_getdqdh(int k, double dir, double h, double q1, double q2);

static int    outlet_readParams(int j, int k, char* tok[], int ntoks);
static double outlet_getFlow(int k, double head);
static double outlet_getInflow(int j);

//=============================================================================

int link_readParams(int j, int type, int k, char* tok[], int ntoks)
//
//  Input:   j     = link index
//           type  = link type code
//           k     = link type index
//           tok[] = array of string tokens
//           ntoks = number of tokens
//  Output:  returns an error code
//  Purpose: reads parameters for a specific type of link from a
//           tokenized line of input data.
//
{
    switch (type)
    {
      case CONDUIT: return conduit_readParams(j, k, tok, ntoks);
      case PUMP:    return pump_readParams(j, k, tok, ntoks);
      case ORIFICE: return orifice_readParams(j, k, tok, ntoks);
      case WEIR:    return weir_readParams(j, k, tok, ntoks);
      case OUTLET:  return outlet_readParams(j, k, tok, ntoks);
      default:      return 0;
    }
}

//=============================================================================

int link_readXsectParams(char* tok[], int ntoks)
//
//  Input:   tok[] = array of string tokens
//           ntoks = number of tokens
//  Output:  returns an error code
//  Purpose: reads a link's cross section parameters from a tokenized
//           line of input data.
//  Formats:
//    Link Shape     Geom1   Geom2   Geom3   Geom4   (Barrels  Culvert)
//    Link IRREGULAR TransectID
//    Link STREET    StreetID
//
{
    int    i, j, k;
    double x[4];

    // --- check for minimum number of tokens
    if (ntoks < 3) return error_setInpError(ERR_ITEMS, "");

    // --- get index of link
    j = project_findObject(LINK, tok[0]);
    if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);

    // --- get code of xsection shape
    k = findmatch(tok[1], XsectTypeWords);
    if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);

    // --- assign default # of barrels to conduit
    if ( Link[j].type == CONDUIT ) Conduit[Link[j].subIndex].barrels = 1;

    // --- assume link is not a culvert
    Link[j].xsect.culvertCode = 0;

    // --- handle IRREGULAR cross section
    if (k == IRREGULAR)
    {
        i = project_findObject(TRANSECT, tok[2]);
        if ( i < 0 ) return error_setInpError(ERR_NAME, tok[2]);
        Link[j].xsect.type = k;
        Link[j].xsect.transect = i;
        return 0;
    }

    // --- handle STREET cross section
    else if (k == STREET_XSECT)
    {
        i = project_findObject(STREET, tok[2]);
        if (i < 0) return error_setInpError(ERR_NAME, tok[2]);
        Link[j].xsect.type = k;
        Link[j].xsect.transect = i;
        return 0;
    }

    // --- handle all other cross section types
    else
    {
        // --- need at least 6 tokens
        if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");

        // --- parse max. depth & shape curve for custom shape
        if ( k == CUSTOM )
        {
            if ( !getDouble(tok[2], &x[0]) || x[0] <= 0.0 )
                return error_setInpError(ERR_NUMBER, tok[2]);
            i = project_findObject(CURVE, tok[3]);
            if ( i < 0 ) return error_setInpError(ERR_NAME, tok[3]);
            Link[j].xsect.type = k;
            Link[j].xsect.transect = i;
            Link[j].xsect.yFull = x[0] / UCF(LENGTH);
        }
        else
        {
            // --- parse & save geometric params
            for (i = 2; i <= 5; i++)
            {
                if (!getDouble(tok[i], &x[i-2]))
                    return error_setInpError(ERR_NUMBER, tok[i]);
            }

            // --- ignore extra parameters for non-conduit open rectangles
            if (Link[j].type != CONDUIT && k == RECT_OPEN)
            {
                x[2] = 0.0;
                x[3] = 0.0;
            }
            if (!xsect_setParams(&Link[j].xsect, k, x, UCF(LENGTH)))
                return error_setInpError(ERR_NUMBER, "");
        }

        // --- parse # of barrels if present
        if ( Link[j].type == CONDUIT && ntoks >= 7 )
        {
            i = atoi(tok[6]);
            if (i <= 0) return error_setInpError(ERR_NUMBER, tok[6]);
            Conduit[Link[j].subIndex].barrels = (char)i;
        }

        // --- parse culvert code if present
        if ( Link[j].type == CONDUIT && ntoks >= 8 )
        {
            i = atoi(tok[7]);
            if (i < 0) return error_setInpError(ERR_NUMBER, tok[7]);
            Link[j].xsect.culvertCode = i;
        }
    }
    return 0;
}

//=============================================================================

int link_readLossParams(char* tok[], int ntoks)
//
//  Input:   tok[] = string tokens
//           ntoks = number of tokens
//  Output:  returns an error code
//  Purpose: reads local loss parameters for a link.
//
//  Format:  LinkID  cInlet  cOutlet  cAvg  FlapGate(YES/NO)  SeepRate
//
{
    int    i, j, k;
    double x[3];
    double seepRate = 0.0;

    if (ntoks < 4) return error_setInpError(ERR_ITEMS, "");
    j = project_findObject(LINK, tok[0]);
    if (j < 0) return error_setInpError(ERR_NAME, tok[0]);
    for (i = 1; i <= 3; i++)
    {
        if (!getDouble(tok[i], &x[i-1]) || x[i-1] < 0.0)
            return error_setInpError(ERR_NUMBER, tok[i]);
    }
    k = 0;
    if ( ntoks >= 5 )
    {
        k = findmatch(tok[4], NoYesWords);
        if (k < 0) return error_setInpError(ERR_KEYWORD, tok[4]);
    }
    if ( ntoks >= 6 )
    {
        if (!getDouble(tok[5], &seepRate))
            return error_setInpError(ERR_NUMBER, tok[5]);
    }
    Link[j].cLossInlet   = x[0];
    Link[j].cLossOutlet  = x[1];
    Link[j].cLossAvg     = x[2];
    Link[j].hasFlapGate  = k;
    Link[j].seepRate     = seepRate / UCF(RAINFALL);
    return 0;
}

//=============================================================================

void link_setParams(int j, int type, int n1, int n2, int k, double x[])
//
//  Purpose: sets parameters for link j.
//
{
    Link[j].node1       = n1;
    Link[j].node2       = n2;
    Link[j].type        = type;
    Link[j].subIndex    = k;
    Link[j].offset1     = 0.0;
    Link[j].offset2     = 0.0;
    Link[j].q0          = 0.0;
    Link[j].qFull       = 0.0;
    Link[j].setting     = 1.0;
    Link[j].targetSetting = 1.0;
    Link[j].hasFlapGate = 0;
    Link[j].qLimit      = 0.0;     // zero means no limit
    Link[j].direction   = 1;

    switch (type)
    {
      case CONDUIT:
        Conduit[k].length    = x[0] / UCF(LENGTH);
        Conduit[k].modLength = Conduit[k].length;
        Conduit[k].roughness = x[1];
        Link[j].offset1      = x[2] / UCF(LENGTH);
        Link[j].offset2      = x[3] / UCF(LENGTH);
        Link[j].q0           = x[4] / UCF(FLOW);
        Link[j].qLimit       = x[5] / UCF(FLOW);
        break;

      case PUMP:
        Pump[k].pumpCurve    = (int)x[0];
        Link[j].hasFlapGate  = FALSE;
        Pump[k].initSetting  = x[1];
        Pump[k].yOn          = x[2] / UCF(LENGTH);
        Pump[k].yOff         = x[3] / UCF(LENGTH);
        Pump[k].xMin         = 0.0;
        Pump[k].xMax         = 0.0;
        break;

      case ORIFICE:
        Orifice[k].type      = (int)x[0];
        Link[j].offset1      = x[1] / UCF(LENGTH);
        Link[j].offset2      = Link[j].offset1;
        Orifice[k].cDisch    = x[2];
        Link[j].hasFlapGate  = (x[3] > 0.0) ? 1 : 0;
        Orifice[k].orate     = x[4] * 3600.0;
        break;

      case WEIR:
        Weir[k].type         = (int)x[0];
        Link[j].offset1      = x[1] / UCF(LENGTH);
        Link[j].offset2      = Link[j].offset1;
        Weir[k].cDisch1      = x[2];
        Link[j].hasFlapGate  = (x[3] > 0.0) ? 1 : 0;
        Weir[k].endCon       = x[4];
        Weir[k].cDisch2      = x[5];
        Weir[k].canSurcharge = (int)x[6];
        Weir[k].roadWidth    = x[7] / UCF(LENGTH);
        Weir[k].roadSurface  = (int)x[8];
        Weir[k].cdCurve      = (int)x[9];
        break;

      case OUTLET:
        Link[j].offset1      = x[0] / UCF(LENGTH);
        Link[j].offset2      = Link[j].offset1;
        Outlet[k].qCoeff     = x[1];
        Outlet[k].qExpon     = x[2];
        Outlet[k].qCurve     = (int)x[3];
        Link[j].hasFlapGate  = (x[4] > 0.0) ? 1 : 0;
        Outlet[k].curveType  = (int)x[5];
        xsect_setParams(&Link[j].xsect, DUMMY, NULL, 0.0);
        break;
    }
}

//=============================================================================

void link_validate(int j)
//
//  Purpose: validates a link's properties.
//
{
    int n;

    // --- convert offset elevations to offset heights if needed
    if (LinkOffsets == ELEV_OFFSET)
        link_convertOffsets(j);

    // --- validate specific link types
    switch (Link[j].type)
    {
      case CONDUIT: conduit_validate(j, Link[j].subIndex); break;
      case PUMP:    pump_validate(j, Link[j].subIndex);    break;
      case ORIFICE: orifice_validate(j, Link[j].subIndex); break;
      case WEIR:    weir_validate(j, Link[j].subIndex);    break;
    }

    // --- check if crest of regulator < invert of downstream node
    switch (Link[j].type)
    {
      case ORIFICE:
      case WEIR:
      case OUTLET:
        if (Node[Link[j].node1].invertElev + Link[j].offset1 <
            Node[Link[j].node2].invertElev)
        {
            if (RouteModel == DW)
            {
                Link[j].offset1 = Node[Link[j].node2].invertElev -
                                   Node[Link[j].node1].invertElev;
                report_writeWarningMsg(WARN10b, Link[j].ID);
            }
            else
                report_writeWarningMsg(WARN10a, Link[j].ID);
        }
    }

    // --- force max depth of end nodes >= link crown if not a storage node
    // --- skip pumps & bottom orifices
    if (Link[j].type == PUMP ||
        (Link[j].type == ORIFICE &&
         Orifice[Link[j].subIndex].type == BOTTOM_ORIFICE)) return;

    // --- adjust upstream node
    n = Link[j].node1;
    if (Node[n].type != STORAGE || Node[n].surDepth > 0.0)
    {
        Node[n].fullDepth = MAX(Node[n].fullDepth,
                                Link[j].offset1 + Link[j].xsect.yFull);
    }

    // --- if conduit, also adjust downstream node
    n = Link[j].node2;
    if ((Node[n].type != STORAGE || Node[n].surDepth > 0.0) &&
         Link[j].type == CONDUIT )
    {
        Node[n].fullDepth = MAX(Node[n].fullDepth,
                                Link[j].offset2 + Link[j].xsect.yFull);
    }
}

//=============================================================================

void link_convertOffsets(int j)
//
//  Purpose: converts offset elevations to offset heights for a link.
//
{
    double elev = Node[Link[j].node1].invertElev;
    Link[j].offset1 = link_getOffsetHeight(j, Link[j].offset1, elev);

    if ( Link[j].type == CONDUIT )
    {
        elev = Node[Link[j].node2].invertElev;
        Link[j].offset2 = link_getOffsetHeight(j, Link[j].offset2, elev);
    }
    else
        Link[j].offset2 = Link[j].offset1;
}

//=============================================================================

double link_getOffsetHeight(int j, double offset, double elev)
//
//  Purpose: finds offset height for one end of a link by subtracting
//           node's invert elevation from the link's offset elevation.
//
{
    if (offset <= MISSING || Link[j].type == PUMP) return 0.0;
    offset -= elev;
    if (offset >= 0.0) return offset;
    if (offset >= -MIN_DELTA_Z) return 0.0;
    report_writeWarningMsg(WARN03, Link[j].ID);
    return 0.0;
}

//=============================================================================

void link_initState(int j)
//
//  Purpose: initializes link state variables at the start of a run.
//
{
    int p;

    // --- hydraulic state
    Link[j].oldFlow   = Link[j].q0;
    Link[j].newFlow   = Link[j].q0;
    Link[j].oldDepth  = 0.0;
    Link[j].newDepth  = 0.0;
    Link[j].oldVolume = 0.0;
    Link[j].newVolume = 0.0;
    Link[j].setting   = 1.0;
    Link[j].targetSetting = 1.0;
    Link[j].timeLastSet = StartDate;
    Link[j].inletControl = FALSE;
    Link[j].normalFlow   = FALSE;

    // --- conduit or pump initialization
    if (Link[j].type == CONDUIT)
        conduit_initState(j, Link[j].subIndex);
    if (Link[j].type == PUMP)
        pump_initState(j, Link[j].subIndex);

    // --- water quality state
    for (p = 0; p < Nobjects[POLLUT]; p++)
    {
        Link[j].oldQual[p] = 0.0;
        Link[j].newQual[p] = 0.0;
        Link[j].totalLoad[p] = 0.0;
    }
}

//=============================================================================

double link_getInflow(int j)
//
//  Purpose: finds total flow entering link j in current time step.
//
{
    // --- return 0 if link's setting is closed
    if (Link[j].setting == 0.0) return 0.0;

    // --- otherwise call the appropriate inflow function by link type
    switch (Link[j].type)
    {
      case CONDUIT: return conduit_getInflow(j);
      case PUMP:    return pump_getInflow(j);
      case ORIFICE: return orifice_getInflow(j);
      case WEIR:    return weir_getInflow(j);
      case OUTLET:  return outlet_getInflow(j);
      default:      return node_getOutflow(Link[j].node1, j);
    }
}

//=============================================================================

void link_setOldHydState(int j)
//
//  Purpose: replaces link's old hydraulic state with the new state.
//
{
    int k;
    Link[j].oldDepth  = Link[j].newDepth;
    Link[j].oldFlow   = Link[j].newFlow;
    Link[j].oldVolume = Link[j].newVolume;

    if (Link[j].type == CONDUIT)
    {
        k = Link[j].subIndex;
        Conduit[k].q1Old = Conduit[k].q1;
        Conduit[k].q2Old = Conduit[k].q2;
    }
}

//=============================================================================

void link_setOldQualState(int j)
//
//  Purpose: replaces link's old quality state with the new state.
//
{
    int p;
    for (p = 0; p < Nobjects[POLLUT]; p++)
    {
        Link[j].oldQual[p] = Link[j].newQual[p];
        Link[j].newQual[p] = 0.0;
    }
}

//=============================================================================

void link_setTargetSetting(int j)
//
//  Purpose: updates a link's target setting based on a simple on/off
//           depth trigger (for pumps).
//
{
    if (Link[j].type == PUMP)
    {
        int k  = Link[j].subIndex;
        int n1 = Link[j].node1;
        Link[j].targetSetting = Link[j].setting;

        // --- if a shutoff depth is specified & upstream depth < shutoff => off
        if (Pump[k].yOff > 0.0 &&
            Link[j].setting > 0.0 &&
            Node[n1].newDepth < Pump[k].yOff)
        {
            Link[j].targetSetting = 0.0;
        }

        // --- if a startup depth is specified & upstream depth > startup => on
        if (Pump[k].yOn > 0.0 &&
            Link[j].setting == 0.0 &&
            Node[n1].newDepth > Pump[k].yOn)
        {
            Link[j].targetSetting = 1.0;
        }
    }
}

//=============================================================================

void link_setSetting(int j, double tstep)
//
//  Purpose: updates link j's current setting to its target setting
//           over the given time step (for orifices, weirs, etc.).
//
{
    if (Link[j].type == ORIFICE)
        orifice_setSetting(j, tstep);
    else if (Link[j].type == WEIR)
        weir_setSetting(j);
    else
        Link[j].setting = Link[j].targetSetting;
}

//=============================================================================

int link_setFlapGate(int j, int n1, int n2, double q)
//
//  Purpose: checks for reverse flow through a flap gate.
//
{
    int n = -1;

    // --- if link has its own flap gate, check for reverse flow
    if (Link[j].hasFlapGate)
    {
        if (q * (double)Link[j].direction < 0.0)
            return TRUE;
    }

    // --- otherwise, if there's an outfall with a flap gate at the 'inflow' end
    if (q < 0.0) n = n2;
    if (q > 0.0) n = n1;
    if (n >= 0 &&
        Node[n].type == OUTFALL &&
        Outfall[Node[n].subIndex].hasFlapGate)
    {
        return TRUE;
    }
    return FALSE;
}

//=============================================================================

void link_getResults(int j, double f, float x[])
//
//  Purpose: retrieves the time-weighted hydraulic & quality state of link j.
//
{
    int    p;
    double y, q, v, u, c;
    double f1 = 1.0 - f;

    // --- time-weighted depth, flow, volume
    y = f1*Link[j].oldDepth  + f*Link[j].newDepth;
    q = f1*Link[j].oldFlow   + f*Link[j].newFlow;
    v = f1*Link[j].oldVolume + f*Link[j].newVolume;

    // --- compute velocity
    u = link_getVelocity(j, q, y);

    // --- compute either fraction full (conduit) or link setting (control)
    if (Link[j].type == CONDUIT)
    {
        if (Link[j].xsect.type != DUMMY)
            c = xsect_getAofY(&Link[j].xsect, y) / Link[j].xsect.aFull;
        else
            c = 0.0;
    }
    else
        c = Link[j].setting;

    // --- override time weighting for pumps toggling between on/off
    if (Link[j].type == PUMP && Link[j].oldFlow*Link[j].newFlow == 0.0)
    {
        if (f >= f1) q = Link[j].newFlow; else q = Link[j].oldFlow;
    }

    // --- apply unit conversions & direction
    y *= UCF(LENGTH);
    v *= UCF(VOLUME);
    q *= UCF(FLOW) * (double)Link[j].direction;
    u *= UCF(LENGTH) * (double)Link[j].direction;

    // --- store results
    x[LINK_DEPTH]    = (float)y;
    x[LINK_FLOW]     = (float)q;
    x[LINK_VELOCITY] = (float)u;
    x[LINK_VOLUME]   = (float)v;
    x[LINK_CAPACITY] = (float)c;

    // --- store pollutant concentrations
    if (!IgnoreQuality)
    {
        for (p = 0; p < Nobjects[POLLUT]; p++)
        {
            c = f1*Link[j].oldQual[p] + f*Link[j].newQual[p];
            x[LINK_QUAL + p] = (float)c;
        }
    }
}

//=============================================================================

void link_setOutfallDepth(int j)
//
//  Purpose: sets the depth at an outfall node connected to link j.
//
{
    int    n, k;
    double z, q, yNorm = 0.0, yCrit = 0.0;

    // --- see which end node is an outfall
    if (Node[Link[j].node2].type == OUTFALL)
    {
        n = Link[j].node2;
        z = Link[j].offset2;
    }
    else if (Node[Link[j].node1].type == OUTFALL)
    {
        n = Link[j].node1;
        z = Link[j].offset1;
    }
    else
        return;

    // --- if link is a conduit, compute normal & critical depth
    if (Link[j].type == CONDUIT)
    {
        k = Link[j].subIndex;
        q = fabs(Link[j].newFlow / Conduit[k].barrels);
        yNorm = link_getYnorm(j, q);
        yCrit = link_getYcrit(j, q);
    }

    // --- set outfall node depth
    node_setOutletDepth(n, yNorm, yCrit, z);
}

//=============================================================================

double link_getYcrit(int j, double q)
//
//  Purpose: finds the critical depth in link j for flow q.
//
{
    return xsect_getYcrit(&Link[j].xsect, q);
}

//=============================================================================

double link_getYnorm(int j, double q)
//
//  Purpose: finds the normal depth in link j for flow q.
//
{
    int    k;
    double s, a, y;

    if (Link[j].type != CONDUIT)        return 0.0;
    if (Link[j].xsect.type == DUMMY)    return 0.0;
    q = fabs(q);
    k = Link[j].subIndex;
    if (q > Conduit[k].qMax) q = Conduit[k].qMax;
    if (q <= 0.0) return 0.0;
    s = q / Conduit[k].beta;
    a = xsect_getAofS(&Link[j].xsect, s);
    y = xsect_getYofA(&Link[j].xsect, a);
    return y;
}

//=============================================================================

double link_getLength(int j)
//
//  Purpose: returns the true length of a link in feet.
//           For natural channels (IRREGULAR), user-supplied length might be
//           just the main channel length. In that case, lengthFactor from the
//           Transect object adjusts it to total floodplain length.
//
{
    if (Link[j].type == CONDUIT)
        return conduit_getLength(j);
    return 0.0;
}

//=============================================================================

double link_getVelocity(int j, double flow, double depth)
//
//  Purpose: computes flow velocity in link j given flow and depth.
//
{
    int    k;
    double area, veloc = 0.0;

    if (depth <= 0.01) return 0.0;
    if (Link[j].type == CONDUIT)
    {
        k = Link[j].subIndex;
        flow /= Conduit[k].barrels;
        area = xsect_getAofY(&Link[j].xsect, depth);
        if (area > FUDGE) veloc = flow / area;
    }
    return veloc;
}

//=============================================================================

double link_getFroude(int j, double v, double y)
//
//  Purpose: computes Froude number for velocity v and depth y in link j.
//
{
    if (Link[j].type != CONDUIT) return 0.0;
    if (y <= FUDGE) return 0.0;
    if (!xsect_isOpen(Link[j].xsect.type) &&
        Link[j].xsect.yFull - y <= FUDGE)
        return 0.0;

    // --- hydraulic depth = area / top width
    y = xsect_getAofY(&Link[j].xsect, y) /
        xsect_getWofY(&Link[j].xsect, y);

    return fabs(v) / sqrt(GRAVITY * y);
}

//=============================================================================

double link_getPower(int j)
//
//  Purpose: returns power consumed by link j in kw
//           from difference in upstream/downstream water levels.
//
{
    int    n1 = Link[j].node1;
    int    n2 = Link[j].node2;
    double dh = (Node[n1].invertElev + Node[n1].newDepth) -
                (Node[n2].invertElev + Node[n2].newDepth);
    double q  = fabs(Link[j].newFlow);
    return fabs(dh) * q / 8.814 * KWperHP;
}

//=============================================================================

double link_getLossRate(int j, int routeModel, double q, double tstep)
//
//  Purpose: returns evaporation + seepage losses in link j (ft^3/s).
//
{
    if (Link[j].type == CONDUIT)
        return conduit_getLossRate(j, routeModel, q, tstep);
    else
        return 0.0;
}

//=============================================================================

char link_getFullState(double a1, double a2, double aFull)
//
//  Purpose: determines if a conduit is full at upstream, downstream,
//           or both ends.
//
{
    if (a1 >= aFull)
    {
        if (a2 >= aFull) return ALL_FULL;
        else             return UP_FULL;
    }
    if (a2 >= aFull) return DN_FULL;
    return 0;
}

Notes:

1. This version incorporates all the documented updates up through Build 5.2.4 (dated 06/12/2023).
2. Comments throughout the code reference both the initial and updated features, bug fixes, and important logic changes.
3. All local helper functions (conduit_readParams, pump_readParams, orifice_readParams, weir_readParams, and outlet_readParams) are included, as well as the main link_XXX functions required by the SWMM engine (reading parameters, validation, setting states, computing inflows, etc.).
4. This file typically resides in the SWMM engine codebase alongside other modules such as headers.h, flowrout.c, dwflow.c, stats.c, etc. The inlet.h and inlet.c files (if in use) are also referenced for inlet handling.

With this code, you should have the latest, consolidated version of link.c suitable for build 5.2.4.

SWMM5 landuse.c Summary

 Below is a concise explanation of how the landuse.c module in SWMM manages pollutant buildup and washoff for different land uses.

---

1. Purpose of landuse.c

This module in SWMM handles how pollutants accumulate (buildup) over time on land surfaces (land uses) and how they are washed off by runoff. In particular, it:

1. Reads and stores user-supplied buildup parameters (functions and coefficients).
2. Reads and stores washoff parameters (functions and coefficients).
3. Computes how pollutant mass accumulates on a land surface over a given dry period, and how that mass is removed by runoff or street sweeping.
4. Calculates how some pollutants can produce co-pollutants during washoff.

---

2. Key Functions

Reading Input

• landuse_readParams():

  • Reads basic land use parameters (like sweeping interval, sweeping removal fraction, initial days since last sweeping).
  • Stores them in the Landuse[j] structure.

• landuse_readPollutParams():

  • Reads pollutant-specific parameters (units, concentration in rain and groundwater, decay/growth rate, co-pollutant index, etc.).
  • Saves them to Pollut[].

• landuse_readBuildupParams():

  • Reads the buildup function type (e.g. power, exponential, saturation, external).
  • Reads function coefficients, normalizer (area or curb length), and sets the maximum days to build up (maxDays).

• landuse_readWashoffParams():

  • Reads the washoff function type (exponential, rating curve, or EMC).
  • Reads washoff coefficients, exponents, plus sweeping and BMP removal efficiencies.

Computing Initial Buildup

• landuse_getInitBuildup():
  • Assigns initial buildup for each pollutant on each land use in a subcatchment.
  • If a subcatchment has an initial buildup or an antecedent dry period, it uses the chosen buildup function to find the starting pollutant load.

Updating Buildup

• landuse_getBuildup():
  • Main routine to incrementally update pollutant buildup for a land use over a time interval.
  • Uses one of the recognized functions (power, exponential, saturation, or external time series).

Washoff Calculation

• landuse_getWashoffLoad():

  • Calculates how much pollutant mass (e.g., lbs or kg) is washed off during a time step, given:
    • The current buildup
    • The runoff rate
    • The washoff function (exponential, rating curve, or EMC)
  • Subtracts the load from the land use’s buildup.
  • Considers possible BMP removal fraction.

• landuse_getWashoffQual():

  • Returns the concentration of pollutant in washoff (mass per volume).
  • Called internally by landuse_getWashoffLoad().

• landuse_getCoPollutLoad():

  • If pollutant p has a co-pollutant (like TSS having a fraction that becomes attached phosphorus), it returns extra mass generated proportionally.
  • This is added into the overall washoff mass.

Internal Helpers

• landuse_getBuildupDays() and landuse_getBuildupMass():

  • Convert from “current buildup mass” to “equivalent # of buildup days,” and vice versa.

• landuse_getExternalBuildup():

  • For “EXTERNAL_BUILDUP” type, obtains a buildup rate from a time series over the time step.

---

3. How Buildup and Washoff are Modeled

Buildup

SWMM supports several buildup function types:

1. POWER_BUILDUP

   $$\text{buildup}(t) = c_1 \, t^{\,c_2} \quad \text{(capped at } c_0 \text{)}$$

2. EXPON_BUILDUP

   $$\text{buildup}(t) = c_0 \left(1 - e^{-\, c_1 t}\right)$$

3. SATUR_BUILDUP

   $$\text{buildup}(t) = \frac{t \, c_0}{c_2 + t}$$

4. EXTERNAL_BUILDUP

   • A user-supplied time series that directly prescribes a buildup rate over time.

The user also selects a normalizer: either per unit area or per unit curb length.

Washoff

Several washoff functions are supported:

1. EXPON_WASHOFF

   $$W = \text{coeff} \times (\text{runoff})^{\,\text{expon}} \times \text{buildup}$$
   • With a typical form:
   $$\frac{d(\text{buildup})}{dt} = - k \cdot R^\beta \cdot \text{buildup}$$

   for some exponent $\beta$.

2. RATING_WASHOFF

   $$W = \text{coeff} \times (\text{runoff})^{(\text{expon})}$$

3. EMC_WASHOFF

   • A fixed concentration (Event Mean Concentration), so mass = $\text{Concentration} \times \text{runoff volume}$.

BMP Efficiency & Sweeping

• A fraction of the washoff mass can be removed by a BMP treatment, specified by bmpEffic.
• Street sweeping can reduce buildup or offset it periodically, governed by an interval, removal fraction, and last-swept day.

---

4. Flow Through the Code

1. Read Land Use & Pollutant Buildup/Washoff
   • The user inputs a series of lines specifying the function type, coefficients, normalizer, and BMP/sweeping parameters.
2. Initialize
   • During subcatchment initialization, landuse_getInitBuildup() assigns initial buildup based on an optional initial mass or an antecedent dry day calculation.
3. During a Runoff Step
   • As time steps progress, landuse_getBuildup() is called to incrementally update buildup.
   • Then, landuse_getWashoffLoad() is called whenever runoff occurs to find how much mass leaves with the stormwater.
   • The rest remains as updated buildup on the surface.
   • BMP or co-pollutant processes are accounted for.
4. Output
   • The resulting washoff load is added to subcatchment runoff quality.
   • Buildup is updated for the next time step.

---

5. Summary

• landuse.c organizes the parameters and calculations for buildup and washoff on each land use within a subcatchment.
• Multiple functional forms (power, exponential, saturation, or time-series-based) are supported for the accumulation of pollutants over dry days.
• Washoff uses exponential, rating curve, or EMC-based formulas to transform buildup mass into pollutant runoff concentration.
• Sweeping and BMP removal features further adjust the net washoff load.
• Co-pollutants allow a fraction of the washoff from one pollutant to create mass in another.

SWMM5 kinwave.c Summary

 Here's a concise summary of how the kinwave.c file in SWMM is structured and used:

---

1. Purpose of kinwave.c

kinwave.c implements the Kinematic Wave flow-routing method for a single conduit, computing the outflow hydrograph at each time step based on the inflow at the upstream end. It solves the finite difference form of the continuity equation for unsteady, 1D flow under the assumptions of negligible inertial and pressure forces.

---

2. Key Routines in kinwave.c

1. kinwave_execute

   • Called each time step for each conduit that is assigned Kinematic Wave routing.
   • Inputs: (i) link index j; (ii) pointer to inflow/outflow; (iii) time step tStep.
   • Outputs: updated outflow (and internally updated flow area) for the conduit.
   • Main steps:
     1. Retrieve the subindex into the Conduit array and relevant cross-section data.
     2. Normalize flows/areas to dimensionless form (i.e., fraction of full flow/area).
     3. Form a finite-difference equation of continuity with the kinematic wave relationship, adjusting for infiltration and evaporation losses.
     4. Solve for the new outlet area (and thus outflow) using solveContinuity.
     5. Convert flows/areas back to dimensional units and store them in the SWMM data structures.

2. solveContinuity

   • Solves the continuity equation: $$f(a) = \beta_1 \cdot S(a) \;+\; C_1 \cdot a \;+\; C_2 \;=\; 0$$ for the outlet dimensionless area $a$.
   • Picks upper/lower bounds on $a$ so that $f(a)$ crosses zero.
   • Calls findroot_Newton (a general-purpose Newton-Raphson routine) from findroot.c.
   • Returns an integer code:
     • $\ge 0$: number of function evaluations used
     • $-1$: continuity equation not solved
     • $-2$: flow above max. flow
     • $-3$: flow below zero

3. evalContinuity

   • The user-supplied function to compute $f(a)$ (and derivative $df(a)$) for findroot_Newton.
   • Combines partial derivatives of Kinematic-Wave formula with partial derivatives of cross-sectional area.
   • $\beta_1, C_1, C_2$ are module-level variables assigned in kinwave_execute.

---

3. Local vs. Shared Variables

Local variables in kinwave_execute:

• q1, q2, qin: normalized flow states at previous time steps / upstream boundary
• a1, a2, ain, aout: normalized cross-section area states

Shared / Module-level variables:

• Beta1, C1, C2: constants in the finite difference continuity equation
• pXsect: pointer to the conduit’s cross-section object
• Afull, Qfull: full flow area and full flow rate used to normalize flows

---

4. Flow Path through the Code

1. kinwave_execute is triggered from flowrout_execute (in flowrout.c) for each conduit that uses Kinematic Wave routing.
2. It sets up dimensionless flows/areas, calls solveContinuity to find the new outlet area, and returns the final computed outflow.
3. solveContinuity calls findroot_Newton, passing evalContinuity for function evaluation and derivative.
4. The root finder finds the new area $a$.
5. kinwave_execute uses a to find outflow $(q_\text{out} = \beta_1 \times S(a))$.
6. SWMM updates the link’s flow/area records.

---

5. Handling Losses

• A small call is made to link_getLossRate with the parameter KW to handle evaporation and infiltration losses.
• The returned value is then subtracted from the flow.

---

6. Newton-Raphson ODE Approach

• The function evalContinuity(a, f, df, p) returns the function $f(a)$ and its derivative $df/da$.
  • $f(a) = \beta_1 \cdot S(a) + C_1 \cdot a + C_2$.
  • The derivative is $df/da = \beta_1 \cdot dS/da + C_1$.

---

7. Practical Limitations

• The kinematic wave approach ignores dynamic wave backwater/inertial effects.
• Appropriate for mild or steep open channels with supercritical flows and no significant flooding or surcharging.

---

8. In Summary

• kinwave.c focuses on one main routine, kinwave_execute, which implements Kinematic Wave routing for a single conduit.
• It uses dimensionless flow/area states, a finite-difference continuity approximation, and a root-finding approach for updated area/outflow.
• The Newton function from findroot.c is plugged in with evalContinuity.
• This approach is simpler computationally than the full Dynamic Wave solver, suitable for many open-channel, non-backwater conditions.

SWMM5 keywords.c Summary

 Below is a concise summary explaining how SWMM’s keywords.c organizes the dictionary of keywords for use across the software, and how each array corresponds to an enumeration in enums.h and is ultimately used in SWMM:

---

1. Purpose of keywords.c

keywords.c provides string arrays containing the textual keywords recognized by SWMM5’s input parser. Each array must match the order of an enumeration (e.g., FlowUnitType, WeirType, etc.) in enums.h. Because SWMM’s input file syntax relies heavily on these keyword strings, the software uses them for:

• Matching tokens found in the input file to enumerated types
• Converting enumerated integer IDs back to user-friendly strings when needed

---

2. Structure of keywords.c

The file defines multiple global, exportable arrays of keywords, for example:

char* BuildupTypeWords[]  = { w_NONE, w_POW, w_EXP, w_SAT, w_EXT, NULL };
char* CurveTypeWords[]    = { w_STORAGE, w_DIVERSION, ..., NULL };
char* DividerTypeWords[]  = { w_CUTOFF, w_TABULAR, w_WEIR, w_OVERFLOW, NULL };
...

Each array:

1. Maps directly to an enumeration type in enums.h. For instance, BuildupTypeWords[] must have the same order as the enumeration BuildupType in enums.h.
2. Terminates with a NULL entry.
3. Uses textual macros (e.g. w_NONE, w_POW, etc.) which are themselves defined in text.h. These macros expand into the actual string values, e.g. #define w_NONE "NONE".

---

3. Matching Between keywords.c and enums.h

In enums.h, one might see something like:

enum BuildupType {
    NO_BUILDUP,
    POWER_BUILDUP,
    EXPON_BUILDUP,
    SATUR_BUILDUP,
    EXTERNAL_BUILDUP
};

In keywords.c, the corresponding array is:

char* BuildupTypeWords[] = {
    w_NONE, w_POW, w_EXP, w_SAT, w_EXT, NULL
};

• The index of "NONE" in BuildupTypeWords is the same as NO_BUILDUP.
• The index of "POW" is the same as POWER_BUILDUP.
• And so forth.

SWMM’s code verifies or fetches an integer from the enumerated list by scanning through the BuildupTypeWords[] array, and when a match is found, sets an internal variable to that enumerated index.

---

4. Typical Usage in SWMM

When SWMM reads an input file line, e.g.:

BUILDUP  SOME_LANDUSE  POW  0.3   ...

The parser calls a function like:

int index = findmatch("POW", BuildupTypeWords);
if (index < 0) error_setInpError(ERR_KEYWORD, "POW");
else {
    // index is now the enumerated integer for POWER_BUILDUP
}

SWMM’s findmatch() function loops over the BuildupTypeWords[] array and compares the user-supplied string token until it finds a match or hits NULL.

---

5. Examples of Keyword Arrays

LinkTypeWords[] corresponds to the LinkType enum:

enum LinkType { CONDUIT, PUMP, ORIFICE, WEIR, OUTLET };
char* LinkTypeWords[] = { w_CONDUIT, w_PUMP, w_ORIFICE, w_WEIR, w_OUTLET };

OptionWords[] matches an InputOptionType enum:

enum InputOptionType {
    FLOW_UNITS, INFIL_MODEL, ROUTE_MODEL, ... MIN_ROUTE_STEP, NUM_THREADS, ...
};
char* OptionWords[] = {
    w_FLOW_UNITS, w_INFIL_MODEL, w_ROUTE_MODEL, ... w_MIN_ROUTE_STEP,
    w_NUM_THREADS, ...
};

SectWords[] matches InputSectionType enum:

enum InputSectionType {
    s_TITLE, s_OPTION, s_FILE, s_RAINGAGE, ..., s_INLET, ...
};
char* SectWords[] = {
    ws_TITLE, ws_OPTION, ws_FILE, ws_RAINGAGE, ..., ws_INLET, NULL
};

and so on.

---

6. Synchronization Requirement

It is critical that the order of keywords in each array not be changed unless the corresponding enumeration is also updated in enums.h. Each array always ends with NULL, which signals to scanning functions that there are no more entries.

---

7. Extensibility

If the SWMM codebase adds a new enumerated value in enums.h, a new string must be appended to the corresponding array in keywords.c at the same position. Similarly, if a string is renamed, the macro in text.h changes and that can alter the user input file syntax.

---

8. Conclusion

SWMM’s keywords.c is essentially a dictionary that ties human-readable keywords in the input file to the internal enumerations in the engine. Each keyword array is strictly parallel to its matching enum. The parser’s findmatch() routine uses these arrays to interpret input tokens and set internal model variables accordingly.

SWMM5 inputrpt.c Summary

 Below is a step-by-step outline summarizing how inputrpt.c writes out a summary of the SWMM input data. This file handles the input data reporting after SWMM has read and processed the input.

---

1. Purpose of inputrpt.c

The inputrpt.c module generates a summary report of the project’s input data in SWMM’s output report file. For example, it displays counts of different objects, plus tables for pollutants, land uses, raingages, subcatchments, nodes, links, cross sections, shapes, transects, etc. The routine is typically called after the model has read in all input data successfully.

---

2. Function Overview

inputrpt_writeInput()

• This is the main function in inputrpt.c. It is typically called once after the input has been parsed successfully.
• It orchestrates writing several summary tables to the SWMM report file Frpt.file:
  1. Element Count
     • Summarizes how many of each major object type exist (raingages, subcatchments, nodes, links, pollutants, land uses).
  2. Pollutant Summary (if any pollutants exist)
  3. Landuse Summary (if any land uses exist)
  4. Raingage Summary (if any raingages exist)
  5. Subcatchment Summary (if any subcatchments exist)
     • Also checks if any LIDs (Low Impact Developments) exist and calls lid_writeSummary().
  6. Node Summary (if any nodes exist)
  7. Link Summary (if any links exist)
  8. Cross Section Summary (for conduits/links)
  9. Shape Summary (if any custom shapes exist)
  10. Transect Summary (if any irregular channel cross sections exist)
  11. Street Summary (if any street-type cross sections exist)

Within each section:

• The code prints column headers using report_writeLine() or fprintf(Frpt.file, ...).
• Iterates over each object in the relevant data array(s) (Pollut, Landuse, Subcatch, Node, Link, Transect, Street, etc.).
• Prints lines listing properties of these objects. For example, for subcatchments: area, width, % impervious, slope, associated raingage, and outlet node.

---

3. Notable Shared Objects and Variables

• Frpt.file: The file pointer to the text-based SWMM report file.
• Title[]: Project title lines (not actually reported here, but part of the general reporting).
• Pollut[]: Array of pollutant objects (ID, units, concentrations, etc.).
• Landuse[]: Array of land uses, which includes sweeping intervals and removal info.
• Gage[]: Array of rain gages.
• Subcatch[]: Array of subcatchments (area, imperv fraction, etc.).
• Node[]: Array of nodes (invert, max depth, ponded area).
• Link[]: Array of links (type, length, slope, roughness).
• Conduit[]: Sub-array used by certain link types to store specific hydraulic info (e.g., slope, length, barrels).
• xsect or Link[i].xsect**: Cross-section geometry details.
• Transect[]: Irregular channel cross sections.
• Street[]: Street cross sections used for “Street Summary” (introduced in SWMM 5.2).
• Shape[]: Array of custom cross-section shapes.
• lid_writeSummary(): Sub-function that writes a summary of any LID usage found in subcatchments.

---

4. Major Reporting Sections

4.1 Element Count

• Prints total counts for raingages, subcatchments, nodes, links, pollutants, land uses.

4.2 Pollutant Summary

• Prints a table with columns for:
  1. Pollutant Name
  2. Units (e.g., mg/L, ug/L, etc.)
  3. Rainfall concentration
  4. Groundwater concentration
  5. Decay coefficient
  6. Co-Pollutant (if any) and fraction

4.3 Landuse Summary

• Prints columns for:
  1. Landuse Name
  2. Sweeping Interval
  3. Maximum Removal
  4. Days since last swept

4.4 Raingage Summary

• Prints columns for:
  1. Gage Name
  2. Data Source (time series or file)
  3. Data Type
  4. Recording Interval

4.5 Subcatchment Summary

• Prints columns for:
  1. Subcatchment Name
  2. Area
  3. Width
  4. % Impervious
  5. % Slope
  6. Assigned rain gage
  7. Outlet
• Also checks if any subcatchment has an LID area, then calls lid_writeSummary().

4.6 Node Summary

• Prints columns for:
  1. Node Name
  2. Node Type (junction, outfall, storage, divider)
  3. Invert Elevation
  4. Max Depth
  5. Ponded Area
  6. If there is external inflow

4.7 Link Summary

• Prints columns for:
  1. Link Name
  2. From Node
  3. To Node
  4. Link Type (conduit, pump, weir, orifice, outlet)
  5. Length, % Slope, Roughness (for conduit types)

4.8 Cross Section Summary

• For each conduit link:
  1. Conduit Name
  2. Cross-section type (irregular, custom, street, or standard shape)
  3. Full depth
  4. Full area
  5. Hydraulic radius
  6. Max width

  7. of barrels

  8. Full flow

4.9 Shape Summary

• For each custom shape in Shape[]:
  • Prints the ID of the underlying shape curve
  • Prints the shape’s stored arrays of area, hydraulic radius, and top width, at fractional depths (the shape curves are tabulated in increments of depth).

4.10 Transect Summary

• For each irregular channel transect in Transect[]:
  • Prints ID
  • Prints tabular arrays for area, hydraulic radius, and top width across fractional depth increments.

4.11 Street Summary

• Similar to Transect Summary, for each Street in Street[]:
  • Prints ID
  • Prints arrays of area, hydraulic radius, and width for the street cross-section geometry.

---

5. Helper Functions and Printing Details

• report_writeLine(text): Writes a line of text to the report file.

• fprintf(Frpt.file, ..., ): Prints formatted data.

• The code often uses:

  fprintf(Frpt.file,"\n  %-20s %10.2f%10.2f ...", SomeString, Value1, Value2);
  

  to format the columns neatly.

• Many checks like if (Nobjects[POLLUT] > 0) are used to conditionally write a section only if that object type is present.

---

6. LID Summary

If lidCount (the total number of subcatchments that have LIDs) is nonzero, a separate function lid_writeSummary() is called. This typically writes details like LID control names, subcatchment area covered, number of units, etc.

---

7. Streets and Inlets

• The code for street cross-section printing occurs near the bottom of the function, within:

  if (Nobjects[STREET] > 0) {...}
  

  It iterates over each Street object, printing out arrays of cross-sectional area, hydraulic radius, and width at incremental depths.

• For inlets (another SWMM5.2 feature), any relevant summary would be handled by the main inlet.c file or the final routing summary, rather than inputrpt.c. Hence, inputrpt.c typically doesn't print a separate "Inlet Summary" table unless specifically coded.

---

8. Conclusion

The inputrpt_writeInput() function is a comprehensive reporter that prints out a textual summary of almost all data read by SWMM. Once the user runs a SWMM project, the output file will contain these sections if Report Input Summary options are enabled in the SWMM project’s REPORT section or in the UI.

SWMM5 input.c Summary

 Below is a step-by-step summary of how reading and parsing input is handled in input.c. This file is responsible for counting, creating, and populating the various SWMM data objects from the text lines in the SWMM input file.

---

1. Overview and Workflow

When SWMM starts reading the input file:

1. input_countObjects()
   The first pass through the file counts the number of each object (nodes, links, subcatchments, patterns, etc.) but does not store details.

2. input_readData()
   The second pass reads detailed properties for each object, storing them in SWMM’s data structures.

Key: Each line is scanned into tokens, and tokens are dispatched to the relevant parser (e.g., a [CONDUITS] line goes to the function that reads a conduit’s parameters).

---

2. Data Structures and Shared Variables

2.1 Shared arrays for counting objects

• Nobjects[]: the final count of each object type.
• Mobjects[]: a running count used on the second pass.
• Ntokens: the number of tokens in the current line.
• Tok[]: pointers to the individual tokens.

2.2 Notable global arrays

• Subcatch, Node, Link, etc.: store object properties.
• Project also references functions to add or find objects.

---

3. Counting Objects

input_countObjects() steps:

1. Initialize Nobjects[], Nnodes[], Nlinks[], etc. to zero.
2. Loop through each line in the input file:
   • Trim comments (semicolon ';').
   • Check if line is [SECTION_HEADING].
     • If so, set current section sect.
     • Otherwise, call addObject(sect, token) with the first token as the ID.
3. addObject increments the corresponding Nobjects[type] counters.
4. If an error occurs (duplicate names, unknown section, etc.), the error is reported.
5. After reaching the end or max errors, final Nobjects[] arrays hold the correct counts for subcatchments, nodes, links, etc.

---

4. Reading Object Data

input_readData() steps:

1. Rewinds the input file.
2. Re-initializes working counters Mobjects[type] to zero.
3. Goes line-by-line:
   • Tokenize the line with getTokens().
   • If [SECTION] found, set sect.
   • Otherwise call parseLine(sect, line).
4. parseLine() dispatches to the correct function based on sect. Examples:
   • s_SUBCATCH -> subcatch_readParams()
   • s_PUMP -> readLink(PUMP) -> link_readParams()
   • s_INLET -> inlet_readDesignParams()
   • etc.
5. Each read function typically uses Mobjects[xxx] to track how many items in that section have been fully read.

---

5. The Tokenizer getTokens()

getTokens(char *s):

1. Removes text after a semicolon ';' (comment).
2. Splits the line into up to MAXTOKS tokens using SEPSTR (spaces, tabs, new line).
3. Also treats quoted text "..." as a single token.
4. Stores pointers to each token in Tok[n].
5. Returns n, the number of tokens.

---

6. Functions that Create or Add Objects

addObject(objType, id):

• Switch on the section (like s_SUBCATCH, s_JUNCTION, etc.).
• Calls project_addObject(type, id, Nobjects[type]) to create the object’s ID in memory.
• Increments Nobjects[type].

This is called during the counting pass (input_countObjects()).

---

7. Parsing Individual Lines of Data

parseLine(sect, line):

1. Switch on the sect (which enumerates [RAINGAGES], [SUBAREAS], [STREETS], etc.).
2. Based on the sect, calls the relevant read function with Tok[] tokens:
   • subcatch_readParams, node_readParams, link_readParams, inflow_readExtInflow, etc.
3. Each read function typically:
   • Uses Mobjects[type] to index the next available slot for that object.
   • Extracts numeric values from tokens with getDouble() or getInt().
   • Fills the object’s data structure fields.

---

8. Handling Specialized Sections

• s_TITLE: lines stored as text in Title[].
• s_TRANSECT: calls transect_readParams(), might parse multiple lines per transect.
• s_CONTROL: processes control rule lines, or named expressions/variables.
• s_LID_CONTROL, s_LID_USAGE: calls specialized LID reading code.
• s_INLET, s_INLET_USAGE: calls inlet reading routines for design vs usage.

---

9. Error Handling and Limits

• If more than MAXERRS input errors occur, SWMM stops reading further.
• Most read functions do sanity checks on numeric values (nonnegative, etc.).
• If ErrorCode is set, final return halts further reading.

---

10. Summary

The input.c file orchestrates the reading of SWMM’s input file in two passes:

1. Count pass: to know how many subcatchments, nodes, links, patterns, etc.
2. Read pass: to parse and store their data.

Within this process, each line is tokenized, sections are recognized, object data is assigned, and any special code (e.g., LID or Inlet or Control rules) is dispatched to their relevant subroutines. This ensures a flexible but consistent way to build SWMM’s internal data structures from the textual input.

SWMM5 inlet.c Summary

 Below is a step-by-step summary describing how street/channel inlets are handled in inlet.c (and referenced in inlet.h) of EPA SWMM 5.2. These methods calculate how much flow in a street or channel is captured by a curb/grate/slotted inlet (or a custom inlet), how much is bypassed, and, for on-sag inlets, how much is conveyed in weir-orifice flow modes.

---

1. Overview & Purpose

SWMM allows inlets—curb, grate, slotted, or custom—to be placed in street-like or open-channel conduits. These inlets divert part of the conveyance flow into a capture (stormwater) node. The goals:

1. Capture Flow at each inlet: how much of the approaching street/ channel flow is intercepted.
2. Bypass: any un-captured flow continues downstream.
3. Backflow: inlets can receive reverse flow from a surcharged capture node.

Key: For on-grade inlets, the captured flow depends on approach flow, street geometry, and inlet properties (HEC-22 design). For on-sag inlets, the local pond depth sets a weir-orifice flow regime.

---

2. Data Structures

2.1 TInletDesign array

typedef struct {
    char*   ID;          // ID name
    int     type;        // e.g., GRATE_INLET, CURB_INLET, etc.
    TGrateInlet    grateInlet; 
    TSlottedInlet  slottedInlet;
    TCurbInlet     curbInlet;
    int            customCurve;  // optional rating curve
} TInletDesign;

• Stores static design parameters (dimensions, type).
• Multiple inlets share the same TInletDesign.

2.2 TInlet Usage object

struct TInlet {
    int     linkIndex;     // conduit index
    int     designIndex;   // which TInletDesign
    int     nodeIndex;     // capture node
    int     numInlets;     // replicate inlets
    int     placement;     // ON_SAG, ON_GRADE, or AUTOMATIC
    double  clogFactor;    // fraction of unclogged area
    double  flowLimit;     // max inlet capacity (cfs)
    double  localDepress;  // local gutter depression (ft)
    double  localWidth;    // local depressed width (ft)

    double  flowFactor;    // coefficient in Q = flowFactor*(spread^2.67)
    double  flowCapture;   // current captured flow (cfs)
    double  backflow;      // current backflow from node (cfs)
    double  backflowRatio; // fraction of node's overflow that returns here

    TInletStats stats;     // performance stats
    TInlet* nextInlet;     // next TInlet in list
};

• Each inlet usage is placed on a particular link.
• Multiple inlets on a link can be assigned. Each can be repeated numInlets times.

Performance Statistics:

• times flow is present, # times inlet captures flow, etc.

2.3 Shared Variables and Arrays

• InletDesigns[]: array of all inlet design definitions from [INLETS].
• FirstInlet: singly-linked list of actual inlets assigned to links.
• InletFlow[]: array indexed by node, holds total flow captured by inlets going to that node.

---

3. Input Reading & Construction

1. inlet_readDesignParams()
   Reads lines from [INLETS] section to fill out TInletDesign objects:

   • GRATE_INLET, CURB_INLET, SLOTTED_INLET, DROP_GRATE_INLET, DROP_CURB_INLET, CUSTOM_INLET, or COMBO_INLET.

2. inlet_readUsageParams()
   Reads lines from [INLET_USAGE] section to create TInlet objects (the actual usage of a design). Key fields: linkID, inletID, nodeID for capture, numInlets, %Clog, etc.

---

4. Initialization and Validation

• inlet_create(nInlets):
  Allocates arrays for inlet designs and captured flows, sets counters.

• inlet_validate():
  Ensures each inlet is placed on a valid link type (e.g., a street cross-section or rectangular/trapezoidal channel).
  If valid, sets node types (BYPASS vs. CAPTURE) and calculates flow factor.
  If invalid, a warning is issued and the inlet is removed from the link.
  Finally, sets how the total overflow from a capture node is distributed among inlets (backflowRatio).

---

5. Flow Capturing Logic

5.1 Invoked Flow Steps

inlet_findCapturedFlows(double tStep) is called before flow routing each time step. Steps:

1. For each inlet, find approach flow:
   • On-grade inlets: approach flow = conduit flow.
   • On-sag inlets: approach flow = inflow to the bypass node.
2. Compute captured flow using formulas from HEC-22 or custom curves.
3. Subtract captured flow from the node’s lateral inflow. Add it to the node receiving the captured flow.
4. Add any node overflow backflow to the bypass node’s inflow.

Key: If more than one inlet is on the same link or node, they’re handled sequentially (since each inlet reduces bypass flow before the next).

---

5.2 On-Grade Inlets

5.2.1 Conduit Geometry

• Street cross-sections:
  • Street slope $Sx$, longitudinal slope $SL$, gutter depression $a$, width $W$.
  • Flow factor Qfactor for the formula $Q = \text{Qfactor} \times T^{2.67}$.
• Trapezoidal or rectangular channels can also have an inlet, but different geometry logic.

5.2.2 On-Grade Approach Flow & Bypass

qApproach = totalFlow / (1 or 2 sides);
for i in numInlets:
   qc = computeSingleInletCapture(qBypassed);
   qCaptured += qc;
   qBypassed -= qc;

5.2.3 On-Grade Grate Inlet

• getGrateInletCapture() uses HEC-22 Equations for splash-over velocity, approach velocity, open area ratio, etc.
  • Ratio of flow in gutter portion vs. entire cross section,
  • Then frontal interception factor and side interception factor.

5.2.4 On-Grade Curb Inlet

• getCurbInletCapture():
  • Finds length needed for full capture using $L_{full} \sim Q^{0.42}$.
  • Then partial capture fraction if the actual curb length < needed length.

5.2.5 On-Grade Slotted Inlet

• Modeled like a curb inlet in on-grade conditions.

5.2.6 Combination Inlets

• Curb + Grate: method calls the curb portion, then leftover approach flow calls the grate portion.

---

5.3 On-Sag Inlets

getOnSagCapturedFlow() uses weir-orifice style logic. The node’s water depth is used. Each inlet either:

• Grate Inlet:
  • Weir flow at low depths vs. orifice flow at higher depth.
• Curb Inlet:
  • Sweeper curb or normal curb equations. HEC-22 derived formula for weir flow if depth < ~1.4 * opening height, else orifice flow.
• Slotted or Custom handle similarly but with simpler or different rating.

If node is surcharged, the node’s overflow can “backflow” into the inlet, allocated by backflowRatio.

---

5.4 Custom Inlets

getCustomCapturedFlow():

• If the design’s curve is type DIVERSION_CURVE, look up capturedFlow = f(approachFlow).
• If the design’s curve is type RATING_CURVE, look up capturedFlow = f(depthAtBypassNode).

---

6. Quality-Related Adjustments

inlet_adjustQualInflows():

• If net captured flow from bypass node j to capture node m is positive, we treat that fraction of pollutant inflow as going to node m, using node j's old pollutant concentration.

inlet_adjustQualOutflows():

• To avoid double-counting captured flow as a system outflow, the mass flow in captured flow is subtracted from outflow totals, etc.

---

7. Statistics and Reporting

updateInletStats():

• Tracks # times the inlet receives approach flow, # times it actually captures flow, peak approach flow, average capture fraction, how often bypass occurs, how often backflow occurs.

inlet_writeStatsReport():

• Writes a “Street Flow Summary” table:
  • Peak Flow, Max Spread, Max Depth, Inlet ID, # inlets, Peak Capture %, etc.

---

8. Summary

inlet.c implements the FHWA HEC-22 methodology for capturing flow on-grade or on-sag using standard inlet types (grate, curb, slotted, or combination) or custom rating curves. It integrates with SWMM’s node flow logic by adjusting node inflows for captured flows and handling backflow from surcharged capture nodes. This allows a more detailed representation of street drainage inflow into stormwater nodes.

SWMM5 infil.c Summary

 Below is a step-by-step explanation of the infiltration implementation in infil.c (along with references to infil.h), from EPA SWMM5. This file handles how SWMM calculates infiltration rates from a subcatchment’s surface into the soil, given different infiltration models (Horton, Modified Horton, Green-Ampt, Modified Green-Ampt, and Curve Number).

---

1. Overview and Purpose

The infiltration code in SWMM 5 can handle multiple infiltration models:

• Horton (classical decay to a minimum rate),
• Modified Horton (adds a new approach to capacity regeneration),
• Green-Ampt (physically based infiltration theory),
• Modified Green-Ampt (slight variation in infiltration capacity),
• Curve Number (SCS-based approach).

At the start, SWMM determines which infiltration model each subcatchment uses. The infiltration object (one per subcatchment) is stored in a union TInfil so that each subcatchment's infiltration parameters can be set appropriately.

Key: A subcatchment’s infiltration is computed at each runoff time step, taking into account the chosen infiltration model, the surface water depth, and infiltration states from previous timesteps (like cumulative infiltration).

---

2. Data Structures

2.1 THorton

typedef struct
{
   double f0;     // initial infiltration rate (ft/s)
   double fmin;   // minimum infiltration rate (ft/s)
   double Fmax;   // maximum total infiltration capacity (ft)
   double decay;  // decay coefficient of infiltration rate (1/s)
   double regen;  // regeneration coefficient of infiltration rate (1/s)
   //--------------
   double tp;     // current time on infiltration curve (s)
   double Fe;     // cumulative infiltration (ft)
} THorton;

What:

• f0: initial infiltration rate, e.g. $f_0$.
• fmin: limiting/minimum infiltration rate.
• Fmax: optional overall maximum infiltration capacity (ft).
• decay: exponential decay parameter $\alpha$.
• regen: exponential “drying” or “regeneration” coefficient.

State:

• tp: how long infiltration has been happening since the subcatchment became wet.
• Fe: total infiltration so far in the event (ft).

2.2 TGrnAmpt

typedef struct
{
   double S;      // capillary suction head (ft)
   double Ks;     // saturated conductivity (ft/s)
   double IMDmax; // max. initial moisture deficit
   //--------------
   double IMD;    // current initial moisture deficit
   double F;      // total infiltration during storm (ft)
   double Fu;     // infiltration into upper soil zone (ft)
   double Lu;     // depth of upper soil zone (ft)
   double T;      // time until next event (s)
   char   Sat;    // saturation flag
} TGrnAmpt;

What:

• S: capillary suction head $(\psi)$.
• Ks: saturated hydraulic conductivity.
• IMDmax: maximum possible deficit in soil moisture.

State:

• IMD: current deficit.
• F: total infiltration volume in current event.
• Fu: infiltration into “upper zone.”
• Lu: thickness of an effective upper zone of the soil.
• T: used in “drying” or “recovery” time for infiltration.
• Sat: whether upper zone is saturated.

2.3 TCurveNum

typedef struct
{
   double Smax;   // maximum infiltration capacity (ft)
   double regen;  // infiltration regeneration constant (1/s)
   double Tmax;   // maximum inter-event time (s)
   //--------------
   double S;      // current infiltration capacity (ft)
   double F;      // current cumulative infiltration (ft)
   double P;      // current cumulative precipitation (ft)
   double T;      // current inter-event time (sec)
   double Se;     // current event infiltration capacity (ft)
   double f;      // infiltration rate in previous step (ft/s)
} TCurveNum;

What:

• Smax: infiltration capacity derived from the SCS curve number formula.
• regen: capacity “regeneration” constant.
• Tmax: maximum time after which infiltration is fully regenerated.

State:

• S: current infiltration capacity.
• F: cumulative infiltration in current event.
• P: cumulative precipitation since event start.
• T: time since end of last event.
• Se: infiltration capacity for the event.
• f: infiltration rate from previous step.

2.4 TInfil Union

typedef union TInfil {
    THorton   horton;
    TGrnAmpt  grnAmpt;
    TCurveNum curveNum;
} TInfil;

• Each subcatchment gets a union that can store one of the infiltration model states.

---

3. Creating and Deleting Infiltration Objects

3.1 infil_create(int n)

void infil_create(int n)
{
    Infil = (TInfil *) calloc(n, sizeof(TInfil));
    if (Infil == NULL) ErrorCode = ERR_MEMORY;
    InfilFactor = 1.0;
    return;
}

• Allocates an array of TInfil of length n (one for each subcatchment).
• Initializes a global factor InfilFactor = 1.0 (used for monthly/time pattern adjustments to hydraulic conductivity).

3.2 infil_delete()

void infil_delete()
{
    FREE(Infil);
}

• Frees the infiltration array.

---

4. Reading and Initializing Infiltration Parameters

4.1 infil_readParams()

int infil_readParams(int m, char* tok[], int ntoks)

• Reads infiltration parameters from one input line.
• m is the default infiltration model. The code checks the last token to see if it’s an overriding infiltration model.
• Depending on the final infiltration method (Horton, Modified Horton, Green-Ampt, Modified Green-Ampt, or Curve Number), a different number of parameters is read.
• Then calls a helper function:
  • horton_setParams()
  • grnampt_setParams()
  • curvenum_setParams()
• On success, assigns the infiltration model to Subcatch[j].infilModel.

4.2 Helper Functions

1. horton_setParams(THorton* infil, double p[])

   • Reads f0, fmin, decay, regen, optional Fmax.
   • Converts rates to ft/sec, times to 1/sec, etc.

2. grnampt_setParams(TGrnAmpt* infil, double p[])

   • Reads S, Ks, IMDmax.
   • Also calculates Lu from Mein’s equation (empirical).

3. curvenum_setParams(TCurveNum* infil, double p[])

   • Reads initial curve number, dryness time, etc.
   • Converts CN to Smax.

4.3 infil_initState(int j)

• Based on the infiltration model, calls:
  • horton_initState(), grnampt_initState(), or curvenum_initState().
• E.g., for Horton, sets tp = 0.0, Fe = 0.0.

---

5. Managing Infiltration State

5.1 Get/Set State

• infil_getState(int j, double x[]): returns infiltration state variables for subcatchment j.
• infil_setState(int j, double x[]): sets infiltration state variables from an array x[].

Inside, the code delegates to one of:

• horton_getState() / horton_setState()
• grnampt_getState() / grnampt_setState()
• curvenum_getState() / curvenum_setState()

This is used for Hot Start file storage.

---

6. Infiltration Factor Adjustment

6.1 infil_setInfilFactor(int j)

void infil_setInfilFactor(int j)
{
    // Start with global hyd. cond. factor
    InfilFactor = Adjust.hydconFactor;

    // If subcatch j has a monthly pattern, override with pattern factor
    if (j >= 0) {
        p = Subcatch[j].infilPattern;
        if (p >= 0 && Pattern[p].type == MONTHLY_PATTERN) {
            m = datetime_monthOfYear(getDateTime(OldRunoffTime)) - 1;
            InfilFactor = Pattern[p].factor[m];
        }
    }
}

• Applies the global hydraulic conductivity factor from climate adjustments, or possibly a subcatchment-specific infiltration pattern.
• This factor is used in infiltration equations to scale Ks, f0, etc.

---

7. Computing the Infiltration Rate

7.1 infil_getInfil()

double infil_getInfil(int j, double tstep, double rainfall,
                      double runon, double depth)

• Master function that calls the correct infiltration model based on Subcatch[j].infilModel.
• Summation: infiltration rate depends on rainfall + runon as the effective “surface inflow” to the infiltration process.
• Some models also treat ponded depth as an available water volume.

---

8. Model-Specific Routines

The code then breaks down the infiltration into each model. They share the same function signature:
double X_getInfil(X* infil, double tstep, double irate, double depth);

Here, irate is net infiltration “supply rate” = rainfall + runon (ft/s), optionally minus evaporation (for some models). In SWMM’s code, some infiltration models subtract surface evaporation from infiltration supply, but the details can vary.

---

8.1 Horton & Modified Horton

8.1.1 horton_getInfil()

1. The infiltration rate decays exponentially from an initial rate f0 to fmin with exponent decay.
2. A tp tracks how long infiltration has been going.
3. If there’s no water available, infiltration capacity “regenerates” at a rate regen.
4. There’s an optional overall limit Fmax.
5. The code does some Newton-Raphson iteration for partial timesteps.

8.1.2 modHorton_getInfil()

1. Simplifies the concept of infiltration capacity as f(t) = f0 - kd * Fe, not the classic exponential but a linear decrease with cumulative infiltration.
2. Also includes a regeneration step when no water is available.

---

8.2 Green-Ampt & Modified Green-Ampt

8.2.1 grnampt_getInfil()

Distinguishes if the upper zone is already saturated (Sat=TRUE) or unsaturated. Then calls:

• grnampt_getSatInfil() if saturated, or
• grnampt_getUnsatInfil() if unsaturated.

grnampt_getUnsatInfil():

1. If no water, tries to “recover” the upper zone.
2. If water supply < Ks, infiltration = supply.
3. If supply > Ks, eventually saturates the upper zone, does partial infiltration unsaturated, partial infiltration saturated.

grnampt_getSatInfil():

1. Uses an integral solution of the Green-Ampt equation (the “F2” function) with Newton iteration to find infiltration volume.
2. The infiltration cannot exceed the actual water supply.

grnampt_getF2() solves the integrated G-A infiltration equation.

---

8.3 Curve Number (SCS)

8.3.1 curvenum_getInfil()

1. If it starts raining, infiltration capacity is CN-based.
2. The infiltration is F1 = P * (1 - P/(P+Se)), a derived expression.
3. If no rainfall, infiltration capacity regenerates with regen.
4. If there's ponded water but no rainfall, it uses the infiltration rate from the last step.

---

9. Summary of Steps

1. Parse infiltration model type from input (Horton / G-A / CN, etc.).
2. Set infiltration object parameters with conversion to consistent units (ft/s, etc.).
3. On each timestep:
   • Possibly update infiltration factor (monthly, subcatch pattern, etc.).
   • Compute infiltration supply = (rainfall + runon) + ponded water / tstep.
   • Call the infiltration model’s function (Horton, G-A, or CN).
   • The function returns infiltration rate in ft/s.
   • Update internal infiltration states (cumulative infiltration, dryness, etc.).
4. At simulation end or after use, free infiltration objects.

Hence, the infil.c file ensures that for each subcatchment, infiltration is computed using the selected model with consistent states and time-step increments, so that SWMM can properly route water from surface runoff into infiltration loss.