Updates and bug fixes for release 5.1.12

src/dwflow.c

@@ -5,6 +5,7 @@
 //   Version:  5.1
 //   Date:     03/20/14   (Build 5.1.001)
 //             03/19/15   (Build 5.1.008)
+//             03/14/17   (Build 5.1.012)
 //   Author:   L. Rossman (EPA)
 //             M. Tryby (EPA)
 //             R. Dickinson (CDM)
@@ -15,6 +16,8 @@
 //   Build 5.1.008:
 //   - Bug in finding if conduit was upstrm/dnstrm full was fixed.
 //
+//   Build 5.1.012:
+//   - Modified uniform loss rate term of conduit momentum equation.
 //-----------------------------------------------------------------------------
 #define _CRT_SECURE_NO_DEPRECATE
 
@@ -76,6 +79,8 @@ void  dwflow_findConduitFlow(int j, int steps, double omega, double dt)
     char   isFull = FALSE;             // TRUE if conduit flowing full
     char   isClosed = FALSE;           // TRUE if conduit closed
 
+
+
     // --- adjust isClosed status by any control action
     if ( Link[j].setting == 0 ) isClosed = TRUE;
 
@@ -208,7 +213,7 @@ void  dwflow_findConduitFlow(int j, int steps, double omega, double dt)
     }
 
     // --- 6. term for evap and seepage losses per unit length
-    dq6 = link_getLossRate(j, qOld, dt) * 1.5 * dt * v / link_getLength(j);    //(5.1.008)
+    dq6 = link_getLossRate(j, qOld, dt) * 2.5 * dt * v / link_getLength(j);    //(5.1.012)
 
     // --- combine terms to find new conduit flow
     denom = 1.0 + dq1 + dq5;

src/flowrout.c

@@ -6,6 +6,7 @@
 //   Date:     03/19/14  (Build 5.1.001)
 //             09/15/14  (Build 5.1.007)
 //             03/19/15  (Build 5.1.008)
+//             03/14/11  (Build 5.1.012)
 //   Author:   L. Rossman (EPA)
 //             M. Tryby (EPA)
 //
@@ -20,6 +21,10 @@
 //   - Determination of node crown elevations moved to dynwave.c.
 //   - Support added for new way of recording conduit's fullness state.
 //
+//   Build 5.1.012:
+//   - Overflow computed in updateStorageState() must be non-negative.
+//   - Terminal storage nodes now updated corectly.
+//
 //-----------------------------------------------------------------------------
 #define _CRT_SECURE_NO_DEPRECATE
 
@@ -570,6 +575,7 @@ void updateStorageState(int i, int j, int links[], double dt)
         {
             Node[i].overflow = (v2 - MAX(Node[i].oldVolume,
                                          Node[i].fullVolume)) / dt;
+            if ( Node[i].overflow < FUDGE ) Node[i].overflow = 0.0;            //(5.1.012)
             if ( !AllowPonding || Node[i].pondedArea == 0.0 )
                 v2 = Node[i].fullVolume;
         }
@@ -632,8 +638,15 @@ void setNewNodeState(int j, double dt)
     int   canPond;                     // TRUE if ponding can occur at node  
     double newNetInflow;               // inflow - outflow at node (cfs)
 
-    // --- storage nodes have already been updated
-    if ( Node[j].type == STORAGE ) return; 
+////  Following section revised for release 5.1.012.  ////                     //(5.1.012)
+    // --- update terminal storage nodes
+    if ( Node[j].type == STORAGE )
+    {	
+	if ( Node[j].updated == FALSE )
+	    updateStorageState(j, Nobjects[LINK], NULL, dt);
+        return; 
+    }
+//////////////////////////////////////////////////////////
 
     // --- update stored volume using mid-point integration
     newNetInflow = Node[j].inflow - Node[j].outflow - Node[j].losses;          //(5.1.007)
@@ -778,6 +791,11 @@ int steadyflow_execute(int j, double* qin, double* qout, double tStep)
             }
         }
         Conduit[k].a2 = Conduit[k].a1;
+
+        Conduit[k].q1Old = Conduit[k].q1;
+        Conduit[k].q2Old = Conduit[k].q2;
+
+
         Conduit[k].q1 = q;
         Conduit[k].q2 = q;
         (*qout) = q * Conduit[k].barrels;

src/globals.h

@@ -8,6 +8,7 @@
 //            09/15/14  (Build 5.1.007)
 //            03/19/15  (Build 5.1.008)
 //            08/01/16  (Build 5.1.011)
+//            03/14/17  (Build 5.1.012)
 //   Author:  L. Rossman
 //
 //   Global Variables
@@ -27,6 +28,9 @@
 //   - Added error message text as a variable. 
 //   - Added elapsed simulation time (in decimal days) variable.
 //   - Added variables associated with detailed routing events.
+//
+//   Build 5.1.012:
+//   - InSteadyState variable made local to routing_execute in routing.c.
 //-----------------------------------------------------------------------------
 
 EXTERN TFile
@@ -87,8 +91,8 @@ EXTERN int
                   SweepEnd,                 // Day of year when sweeping ends
                   MaxTrials,                // Max. trials for DW routing
                   NumThreads,               // Number of parallel threads used //(5.1.008)
-                  NumEvents,                // Number of detailed events       //(5.1.011)
-                  InSteadyState;            // System flows remain constant    //(5.1.011)
+                  NumEvents;                // Number of detailed events       //(5.1.011)
+                //InSteadyState;            // System flows remain constant    //(5.1.012)
 
 EXTERN double
                   RouteStep,                // Routing time step (sec)

src/lid.c

@@ -10,6 +10,7 @@
 //             04/30/15   (Build 5.1.009)
 //             08/05/15   (Build 5.1.010)
 //             08/01/16   (Build 5.1.011)
+//             03/14/17   (Build 5.1.012)
 //   Author:   L. Rossman (US EPA)
 //
 //   This module handles all data processing involving LID (Low Impact
@@ -60,6 +61,9 @@
 //   - The top of the storage layer is no longer used as a limit for an
 //     underdrain offset thus allowing upturned drains to be modeled.
 //   - Column headings for the detailed LID report file were modified.
+//
+//   Build 5.1.012:
+//   - Redefined initialization of wasDry for LID reporting.
 //-----------------------------------------------------------------------------
 #define _CRT_SECURE_NO_DEPRECATE
 
@@ -1898,5 +1902,6 @@ void initLidRptFile(char* title, char* lidID, char* subcatchID, TLidUnit* lidUni
     for ( i = 1; i < colCount; i++) fprintf(f, "\t%s", line9);
 
     //... initialize LID dryness state
-    lidUnit->rptFile->wasDry = TRUE;
+    lidUnit->rptFile->wasDry = 1;                                              //(5.1.012)
+    strcpy(lidUnit->rptFile->results, "");                                     //(5.1.012)
 }

src/lid.h

@@ -6,6 +6,7 @@
 //   Date:    03/20/14   (Build 5.1.001)
 //            03/19/15   (Build 5.1.008)
 //            08/01/16   (Build 5.1.011)
+//            03/14/17   (Build 5.1.012)
 //   Author:  L. Rossman (US EPA)
 //
 //   Public interface for LID functions.
@@ -19,6 +20,9 @@
 //   - Water depth replaces moisture content for LID's pavement layer. 
 //   - Arguments for lidproc_saveResults() modified.
 //
+//   Build 5.1.012:
+//   - Redefined meaning of wasDry in TLidRptFile structure.
+//
 //-----------------------------------------------------------------------------
 
 #ifndef LID_H
@@ -141,7 +145,7 @@ typedef struct
 typedef struct
 {
     FILE*     file;               // file pointer
-    int       wasDry;             // true if LID was dry                       //(5.1.008)
+    int       wasDry;             // number of successive dry periods          //(5.1.012)
     char      results[256];       // results for current time period           //(5.1.008)
 }   TLidRptFile;
 

src/lidproc.c

@@ -10,6 +10,7 @@
 //             04/30/15   (Build 5.1.009)
 //             08/05/15   (Build 5.1.010)
 //             08/01/16   (Build 5.1.011)
+//             03/14/17   (Build 5.1.012)
 //   Author:   L. Rossman (US EPA)
 //
 //   This module computes the hydrologic performance of an LID (Low Impact
@@ -46,6 +47,12 @@
 //     physically meaningful results.
 //   - Reporting of detailed results re-written.
 //
+//   Build 5.1.012:
+//   - Modified upper limit for soil layer percolation.
+//   - Modified upper limit on surface infiltration into rain gardens.
+//   - Modified upper limit on drain flow for LIDs with storage layers.
+//   - Used re-defined wasDry variable for LID reports to fix duplicate lines.
+//
 //-----------------------------------------------------------------------------
 #define _CRT_SECURE_NO_DEPRECATE
 
@@ -378,14 +385,15 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe
     if ( SurfaceInflow  < MINFLOW &&
          SurfaceOutflow < MINFLOW &&
          StorageDrain   < MINFLOW &&
-         StorageExfil   < MINFLOW && 
-         totalEvap      < MINFLOW ) isDry = TRUE;
+         StorageExfil   < MINFLOW &&
+		 totalEvap      < MINFLOW
+       ) isDry = TRUE;
 
     //... update status of HasWetLids
     if ( !isDry ) HasWetLids = TRUE;
 
-    //... write results to LID report file if at a reporting time              //(5.1.011)
-    if ( lidUnit->rptFile && fabs(NewRunoffTime - ReportTime) < 10. )          //(5.1.011)
+    //... write results to LID report file                                     //(5.1.012)
+    if ( lidUnit->rptFile )                                                    //(5.1.012)
     {
         //... convert rate results to original units (in/hr or mm/hr)
         ucf = ucfRainfall;
@@ -406,11 +414,13 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe
         rptVars[STOR_DEPTH] = theLidUnit->storageDepth*ucf;
 
         //... if the current LID state is wet but the previous state was dry
-        //    then write the saved previous results to the report file thus
-        //    marking the end of a dry period
-        if ( !isDry && theLidUnit->rptFile->wasDry )
+        //    for more than one period then write the saved previous results   //(5.1.012)
+        //    to the report file thus marking the end of a dry period          //(5.10012)
+        if ( !isDry && theLidUnit->rptFile->wasDry > 1)                        //(5.1.012)
+        {
             fprintf(theLidUnit->rptFile->file, "%s",
-            theLidUnit->rptFile->results);
+				  theLidUnit->rptFile->results);
+        }
 
         //... write the current results to a string which is saved between
         //    reporting periods
@@ -428,20 +438,25 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe
         {
             //... if the previous state was wet then write the current
             //    results to file marking the start of a dry period
-            if ( !theLidUnit->rptFile->wasDry )
+            if ( theLidUnit->rptFile->wasDry == 0 )                            //(5.1.012)
             {
-                fprintf(theLidUnit->rptFile->file, "%s", theLidUnit->rptFile->results);
+                fprintf(theLidUnit->rptFile->file, "%s",
+					theLidUnit->rptFile->results);
             }
-            theLidUnit->rptFile->wasDry = TRUE;
+
+            //... increment the number of successive dry periods               //(5.1.012)
+            theLidUnit->rptFile->wasDry++;                                     //(5.1.012)
         }
 
         //... if the current LID state is wet
         else
-
         {
             //... write the current results to the report file
-            fprintf(theLidUnit->rptFile->file, "%s", theLidUnit->rptFile->results);
-            theLidUnit->rptFile->wasDry = FALSE;
+			fprintf(theLidUnit->rptFile->file, "%s",
+			    theLidUnit->rptFile->results);
+
+            //... re-set the number of successive dry periods to 0             //(5.1.012)
+            theLidUnit->rptFile->wasDry = 0;                                   //(5.1.012)
         }
     }
 }
@@ -518,7 +533,7 @@ void greenRoofFluxRates(double x[], double f[])
 
     //... limit perc rate by available water
     availVolume = (soilTheta - soilFieldCap) * soilThickness;
-    maxRate = SurfaceInfil - SoilEvap + MAX(availVolume, 0.0) / Tstep;
+    maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;                        //(5.1.012)
     SoilPerc = MIN(SoilPerc, maxRate);
     SoilPerc = MAX(SoilPerc, 0.0);
 
@@ -542,8 +557,8 @@ void greenRoofFluxRates(double x[], double f[])
     else
     {
         //... limit drainmat outflow by available storage volume
-        maxRate = storageDepth * storageVoidFrac / Tstep +
-                  SoilPerc - StorageEvap;
+        maxRate = storageDepth * storageVoidFrac / Tstep - StorageEvap;        //(5.1.012)
+        if ( storageDepth >= storageThickness ) maxRate += SoilPerc;           //(5.1.012)
         maxRate = MAX(maxRate, 0.0);
         StorageDrain = MIN(StorageDrain, maxRate);
 
@@ -618,7 +633,7 @@ void biocellFluxRates(double x[], double f[])
 
     //... limit perc rate by available water
     availVolume =  (soilTheta - soilFieldCap) * soilThickness;
-    maxRate = SurfaceInfil - SoilEvap + MAX(availVolume, 0.0) / Tstep;
+    maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;                        //(5.1.012)
     SoilPerc = MIN(SoilPerc, maxRate);
     SoilPerc = MAX(SoilPerc, 0.0);
 
@@ -640,7 +655,15 @@ void biocellFluxRates(double x[], double f[])
         maxRate = MIN(SoilPerc, StorageExfil);
         SoilPerc = maxRate;
         StorageExfil = maxRate;
-    }
+
+////  Following code segment added to release 5.1.012  ////                    //(5.1.012)
+        //... limit surface infil. by unused soil volume
+        maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
+                  SoilPerc + SoilEvap;
+        SurfaceInfil = MIN(SurfaceInfil, maxRate);
+//////////////////////////////////////////////////////////
+
+	}
 
     //... storage & soil layers are full
     else if ( soilTheta >= soilPorosity && storageDepth >= storageThickness )
@@ -668,13 +691,14 @@ void biocellFluxRates(double x[], double f[])
     {
         //... limit storage exfiltration by available storage volume
         maxRate = SoilPerc - StorageEvap + storageDepth*storageVoidFrac/Tstep;
-        maxRate = MAX(maxRate, 0.0);
         StorageExfil = MIN(StorageExfil, maxRate);
+        StorageExfil = MAX(StorageExfil, 0.0);
 
         //... limit underdrain flow by volume above drain offset
         if ( StorageDrain > 0.0 )
         {
-            maxRate = SoilPerc - StorageExfil - StorageEvap;
+            maxRate = -StorageExfil - StorageEvap;                              //(5.1.012)
+            if ( storageDepth >= storageThickness) maxRate += SoilPerc;         //(5.1.012)
             if ( theLidProc->drain.offset <= storageDepth )
             {
                 maxRate += (storageDepth - theLidProc->drain.offset) *
@@ -766,14 +790,15 @@ void trenchFluxRates(double x[], double f[])
     StorageExfil = MIN(StorageExfil, maxRate);
     StorageExfil = MAX(StorageExfil, 0.0);
 
-    //... limit underdrain flow by available storage volume
+    //... limit underdrain flow by volume above drain offset
     if ( StorageDrain > 0.0 )
     {
-        maxRate = StorageInflow - StorageExfil - StorageEvap;
+        maxRate = -StorageExfil - StorageEvap;                                 //(5.1.012)
+        if (storageDepth >= storageThickness ) maxRate += StorageInflow;       //(5.1.012)
         if ( theLidProc->drain.offset <= storageDepth )
         {
             maxRate += (storageDepth - theLidProc->drain.offset) *
-                           storageVoidFrac/Tstep;
+                       storageVoidFrac/Tstep;
         }
         maxRate = MAX(maxRate, 0.0);
         StorageDrain = MIN(StorageDrain, maxRate);
@@ -869,7 +894,7 @@ void pavementFluxRates(double x[], double f[])
     {
         SoilPerc = getSoilPercRate(soilTheta);
         availVolume = (soilTheta - soilFieldCap) * soilThickness;
-        maxRate = PavePerc - SoilEvap + MAX(availVolume, 0.0) / Tstep;
+        maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;                    //(5.1.012)
         SoilPerc = MIN(SoilPerc, maxRate);
         SoilPerc = MAX(SoilPerc, 0.0);
     }
@@ -982,8 +1007,9 @@ void pavementFluxRates(double x[], double f[])
         //... limit underdrain flow by volume above drain offset
         if ( StorageDrain > 0.0 )
         {
-            maxRate = SoilPerc - StorageExfil - StorageEvap;
-            if ( theLidProc->drain.offset <= storageDepth )
+            maxRate = -StorageExfil - StorageEvap;                             //(5.1.012)
+            if (storageDepth >= storageThickness ) maxRate += SoilPerc;        //(5.1.012)
+            if ( theLidProc->drain.offset <= storageDepth ) 
             {
                 maxRate += (storageDepth - theLidProc->drain.offset) *
                            storageVoidFrac/Tstep;