      SUBROUTINE CalculateMode (dataset, numrow, mode)
      INCLUDE 'STDHDR.FOR'
      REAL dataset(0:maxc)
      INTEGER numrow
      REAL mode

      INTEGER Num
      LOGICAL even

      even = .TRUE.
      IF (MOD(numrow,2) .NE. 0) even = .FALSE.
      CALL SortDataSet(dataset, numrow)
      IF (even) THEN
         Num = numrow / 2 - 1
         mode = (dataset(Num) + dataset(Num + 1)) / 2.0
      ELSE
          Num = numrow / 2
          mode = dataset(Num)
      END IF

      END  !Subroutine  CalculateMode



      SUBROUTINE StatsMinMax (dataset, numdat, minval, maxval, range)
      INCLUDE 'STDHDR.FOR'
      REAL dataset(0:maxc )
      INTEGER numdat,i
      REAL minval,maxval,range

      minval = dataset(0)
      maxval = dataset(0)
      DO i = 1 , numdat - 1
         IF (dataset(i) .LT. minval)  minval = dataset(i)
         IF (dataset(i) .GT. maxval)  maxval = dataset(i)
      END DO
      range = maxval - minval

      END  !Subroutine  StatsMinMax



      SUBROUTINE MatMean (datary, numrow, sumx, mean)
      INCLUDE 'STDHDR.FOR'
      REAL datary(0:maxc)
      INTEGER numrow,i
      REAL sumx, mean

      sumx = 0
      DO i = 0 , numrow - 1
         sumx = sumx + datary(i)
      END DO
      mean = sumx / numrow

      END  !Subroutine  matmean



      SUBROUTINE MatStd (datary, mean, numrow, variance, stddev)
      INCLUDE 'STDHDR.FOR'
      REAL datary(0:maxc)
      INTEGER numrow
      REAL variance, stddev,mean

      REAL xsqr
      INTEGER i

      xsqr = 0
      DO i = 0 , numrow - 1
         xsqr = xsqr + datary(i) ** 2
      END DO
      variance = (xsqr - numrow * mean ** 2) / (numrow - 1)
      stddev = SQRT(ABS(variance))

      END  !Subroutine  MatStd



      SUBROUTINE SortDataSet (dataset, numrow)
      INCLUDE 'STDHDR.FOR'
      REAL dataset(0:maxc)
      INTEGER numrow

      INTEGER k, j
      REAL TempX
      LOGICAL abort

      IF (numrow .GT. 1) THEN
         DO j = 0 , numrow - 1
            abort = .FALSE.
            TempX = dataset(j)
            k = j - 1
            DO WHILE ( ( abort .EQV. .FALSE.) .AND. (k .GE. 0))
               IF (TempX .LT. dataset(k)) THEN
                  dataset(k + 1) = dataset(k)
                  k = k - 1
               ELSE
                  abort = .TRUE.
               END IF
               dataset(k + 1) = TempX
            END DO

         END DO
      END IF

      END  !Subroutine  sortdataset



      SUBROUTINE SummaryStats (dataset, numobs, numcol, minima,
     + maxima, range, sumxx, mean, variance, stddev, semean, mode)
      INCLUDE 'STDHDR.FOR'
      REAL dataset(0:maxr,0:maxc)
      INTEGER numobs, numcol, i, j
      REAL minima(0:maxc), maxima(0:maxc), range(0:maxc)
      REAL sumxx(0:maxc), mean(0:maxc), variance(0:maxc)
      REAL stddev(0:maxc), semean(0:maxc), mode(0:maxc)
      REAL tempv(0:maxc)

      DO j = 0 , numcol - 1
         DO i = 0 , numobs - 1
            tempv(i) = dataset(i, j)
         END DO
         CALL CalculateMode(tempv, numobs, mode(j))
         CALL StatsMinMax(tempv, numobs, minima(j),maxima(j),range(j))
         CALL MatMean(tempv, numobs, sumxx(j), mean(j))
         CALL MatStd(tempv, mean(j), numobs, variance(j), stddev(j))
         semean(j) = stddev(j) / SQRT(numobs)
      END DO

      END  !Subroutine  summarystats

