Sunday, August 12, 2018

Moment Calculations in the SWMM4 Stats Block

   Subject:  This is the code from the Stats Routine of SWMM3 and SWMM4

SUBROUTINE MOMENT(J,NLABEL)
C=======================================================================
C     This subroutine is called by the STATS Block.  For a given series
C     of event data, it computes and print the mean, variance, standard
C     deviation, coefficient of variation, and coefficient of skewness.
C
C     Updated September, 1990 by Laura B. Terrell (CDM).
C     Updated December,  1990 by Bob Dickinson.
C     WCH, 8/1/95.  Change rainfall station ID (LOCRN) to character.
C=======================================================================
      INCLUDE 'TAPES.INC'
      INCLUDE 'INTER.INC'
      INCLUDE 'STCOM.INC'
      DOUBLE PRECISION XEVNTS,SUMX,SUMX2,SUMX3,AMED,ACV,AXBAR,
     1                 ASDX,NUMER,FACTOR,DENOM,DPARM
C=======================================================================
C     Compute sum, sum of squares, and sum of cubes.
C=======================================================================
      JTEMP    = 0
      SUMX     = 0.0
      SUMX2    = 0.0
      SUMX3    = 0.0
      DO 100 I = 1,NEVNTS
      DPARM    = PARAM(I,J)
      IF(INLOG.GT.0.AND.DPARM.NE.0.0) DPARM = DLOG(DPARM)
      SUMX     = SUMX  + DPARM
      SUMX2    = SUMX2 + DPARM*DPARM
      SUMX3    = SUMX3 + DPARM*DPARM*DPARM
100   CONTINUE
C=======================================================================
C     Calculate moments.
C=======================================================================
      XEVNTS = DFLOAT(NEVNTS)
      XBAR   = SUMX/XEVNTS
      IF(NEVNTS.GT.1) THEN
                      VARX   = (SUMX2 - XEVNTS*(XBAR)**2)/(XEVNTS-1.0)
                      SDX    = SQRT(ABS(VARX))
                                      CV  = 0.0
                      IF(XBAR.NE.0.0) CV  = SDX/XBAR
                      ELSE
                      CV     = 0.0
                      SDX    = 0.0
                      VARX   = 0.0
                      ENDIF
      IF(NEVNTS.GT.2.AND.XBAR.NE.0.0) THEN
                   NUMER  = SUMX3/XEVNTS - 3.0*SUMX2*XBAR/XEVNTS +
     1                      2.0*(XBAR)**3
                   FACTOR = DSQRT(XEVNTS*(XEVNTS-1.0)) / (XEVNTS-2.0)
                   DENOM  = SUMX2/XEVNTS - XBAR*XBAR
                   IF(DENOM.LE.0.0) DENOM  = 1.0
                   DENOM  = DENOM**1.5
                   COSKEW = NUMER * FACTOR / DENOM
                   ELSE
                   COSKEW = 0.0
                   ENDIF
C=======================================================================
C     If analysis is for a pollutant, go to print statement for pollutants
C=======================================================================
      IF(NLABEL.EQ.0) THEN
                      J1    = J
                      JTEMP = 1
C#### WCH, 8/1/95.  LOCRN NOW CHARACTER.
C####                      IF(LOCRN.GT.0) J1 = J + 10
                      IF(LOCRN.NE.' ') J1 = J + 10
                      ELSE
                      J1    = J + 5
                      ENDIF
C=======================================================================
C     Print parameter identifier.
C=======================================================================
C     Perform log transform calculations on log emc data if desired.
C     This will override other flow and pollutant output.
C
C     From previously calculated log data (XBAR & SDX), arithmetic
C     statistics may be calculated.
C
C     XBAR  = Logarithmetic mean
C     SDX   = Logarithmetic standard deviation
C     AXBAR = Arithmetic mean
C     ACV   = Arithmetic coefficient of variation
C     ASDX  = Arithmetic standard deviation
C     AMED  = Arithmetic median
C=======================================================================
      IF(INLOG.EQ.1) THEN
           AXBAR = DEXP(XBAR + 0.5*SDX*SDX)
           ACV   = DSQRT(EXP(SDX*SDX) - 1.0)
           ASDX  = AXBAR * ACV
           AMED  = DEXP(XBAR)
           IF(JTEMP.EQ.1) THEN
              IF(LOCRQ.GT.0) WRITE(N6,1000) RAINF(2),PARLAB(J1),
     +                               XBAR,SDX,AXBAR,ASDX,ACV,AMED
C#### WCH, 8/1/95.  LOCRN NOW CHARACTER.
C####              IF(LOCRN.GT.0) WRITE(N6,1000) RAINF(1),PARLAB(J1),
              IF(LOCRN.NE.' ') WRITE(N6,1000) RAINF(1),PARLAB(J1),
     +                               XBAR,SDX,AXBAR,ASDX,ACV,AMED
              ELSE
              WRITE(N6,1000) PNAME(NLABEL),PARLAB(J1),
     +                               XBAR,SDX,AXBAR,ASDX,ACV,AMED
              ENDIF
           ELSE
           IF(JTEMP.EQ.1) THEN
              IF(LOCRQ.GT.0) WRITE(N6,1000) RAINF(2),PARLAB(J1),XBAR,
     +                                            VARX,SDX,CV,COSKEW
C#### WCH, 8/1/95.  LOCRN NOW CHARACTER.
C####            IF(LOCRN.GT.0) WRITE(N6,1000) RAINF(1),PARLAB(J1),XBAR,
              IF(LOCRN.NE.' ') WRITE(N6,1000) RAINF(1),PARLAB(J1),XBAR,
     +                                            VARX,SDX,CV,COSKEW
              ELSE
              WRITE(N6,1000) PNAME(NLABEL),PARLAB(J1),XBAR,VARX,
     +                                            SDX,CV,COSKEW
              ENDIF
           ENDIF
C=======================================================================
C     Print note regarding magnitude units for NDIM = 2.
C=======================================================================
      IF(NLABEL.EQ.0)       RETURN
      IF(NDIM(NLABEL).NE.2) RETURN
      IF(J.EQ.1) THEN
                 IF(METRIC.EQ.1) THEN
                                 WRITE(N6,1022) PUNIT(NLABEL)
                                 ELSE
                                 WRITE(N6,1026) PUNIT(NLABEL)
                                 ENDIF
                 ENDIF
      IF(J.GE.2.AND.J.LE.3) THEN
                 IF(METRIC.EQ.1) THEN
                                 WRITE(N6,1110) PUNIT(NLABEL)
                                 ELSE
                                 WRITE(N6,1130) PUNIT(NLABEL)
                                 ENDIF
                 ENDIF
      IF(J.GE.4) WRITE(N6,1160) PUNIT(NLABEL)
      RETURN
C=======================================================================
1000  FORMAT(1X,A8,A20,6(1X,1PG11.4))
1022  FORMAT(/,1X,'Note:  Magnitude has units of (cubic feet) * ',
     +   A8,'.  See user''s manual for explanation.')
1026  FORMAT(/,1X,'Note:  Magnitude has units of (liters) * ',
     +   A8,'.  See user''s manual for explanation.')
1110  FORMAT(/,1X,'Note:  Magnitude has units of (cfs) * ',
     +   A8,'  See user''s manual for explanation.')
1130  FORMAT(/,1X,'Note:  Magnitude has units of (liters/s) * ',
     +   A8,'.  See user''s manual for explanation.')
1160  FORMAT(/,1X,'Note:  Magnitude has units of ',
     +   A8,'.  See user''s manual for explanation.')
C=======================================================================
      END

Privileged and Confidential Communication: This electronic mail communication and any documents included hereto may contain confidential and privileged material for the sole use of the intended recipient(s) named above. If you are not the intended recipient (or authorized to receive for the recipient) of this message, any review, use, distribution or disclosure by you or others is strictly prohibited. Please contact the sender by reply email and delete and/or destroy the accompanying message.

No comments:

GitHub code and Markdown (MD) files Leveraging

 To better achieve your goal of leveraging your GitHub code and Markdown (MD) files for your WordPress blog or LinkedIn articles, consider t...