CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TABLE_PM_EDB_AVERAGE( OUTPUT, UNIT, WHO_CALLED, AVERAGE_EPOCH )
C+
C
C FUNCTIONAL DESCRIPTION:
C
C     This subroutine computes the average of the counts that were accumulated
C     by the routine TABLE_PM_EDB and conditionally writes them to an output
C     file.  It is intended to be called at the end of an averaging interval. 
C     It is also called at the end of the program to generate a partial
C     average, if needed.  The accumulated counts are cleared before return.
C
C FORMAL ARGUMENTS:
C
C     OUTPUT:
C	L*4	  I	If true the accumulated values will be averaged and,
C			depending on the value of PRINT, written to the
C			output file.  If false no averaging or writing will
C			be done; the variables used for accumulation will be
C			cleared, however.
C
C     UNIT:
C	I*4	  I	The Fortran unit number to which the averaged values
C			will be written.
C
C     WHO_CALLED:
C	I*4	  I	An integer informing this routine which application
C			is calling it.  Allowable values are 1 (PMRATES) and
C			2 (PMSW).
C
C     AVERAGE_EPOCH:
C	R*8	  I	The Epoch time of the start of the averaging interval.
C
C INTERFACE:
!
!	INTERFACE TABLE_PM_EDB_AVERAGE
!	  SUBROUTINE TABLE_PM_EDB_AVERAGE( OUTPUT, UNIT, WHO_CALLED, AVERAGE_EPOCH )
!	    LOGICAL, INTENT(IN) :: OUTPUT
!	    INTEGER, INTENT(IN) :: UNIT, WHO_CALLED
!	    REAL(KIND=8), INTENT(IN) :: AVERAGE_EPOCH
!	  END SUBROUTINE TABLE_PM_EDB_AVERAGE
!	END INTERFACE TABLE_PM_EDB_AVERAGE
C
C IMPLICIT INPUTS:
C
C     Contents of common block SUM, which contains the sums and the time of
C     the first summed item.
C
C     The array W is assumed to contain zeroes upon entry.
C
C IMPLICIT OUTPUTS:
C
C     None
C
C SIDE EFFECTS:
C
C     A line is written to the output file.
C
C     The arrays PM_EQ, PM_RATE, UNC_PM_RATE, PM_COUNTER, PM_T1, and TCOUNT,
C     as well as several other variables in common block SUM and the local
C     array W, are cleared to zero before returning.
C
C COMMON BLOCKS:
C
C     SUM	Contains the variables used to hold accumulated sums.
C     FORM	Contains the format statement and control variables to be used
C		to write to the output file.
C
C INCLUDE FILES:
C
C     FORM.INC	defines the variables in the common block FORM
C
C DESIGN:
C
C     All arguments were originally local variables in the routine TABLE_PM_EDB
C     except OUTPUT.  This is set true (the normal case) if the accumulated
C     rates are to be averaged and printed.  This was added to allow
C     TABLE_PM_EDB_AVERAGE to also clear the arrays in which sums were kept. 
C     Since this operation is done unconditionally, leaving its code in
C     TABLE_PM_EDB would have required another set of arguments or another
C     COMMON block to make those arrays, which are local to
C     TABLE_PM_EDB_AVERAGE, accessible to TABLE_PM_EDB. 
C
C     Under some circumstances averages will be computed but not written.  This
C     is controlled by the argument OUTPUT, which is set in the calling module,
C     and the variable PRINT, which is part of the common block SUM and is set
C     both in the calling module and in this module.
C
C This version:   1.5    7-MAR-2000
C
C CREATED:
C
C     13-SEP-1999    Larry Bleau
C
C DEVELOPMENT HISTORY:
C
C	 Date     | Name  | Description
C ----------------+-------+-----------------------------------------------------
C  14-SEP-1999 1.0   LRB    split off from TABLE_PM_EDB, have separate
C                           routine to just average and print
C  15-SEP-1999       LRB    move reset of PM_QUALITY to TABLE_PM_EDB_AVERAGE
C  15-SEP-1999       LRB    separate TABLE_PM_EDB_AVERAGE from TABLE_PM_EDB,
C			    put into different file; add comments
C  12-NOV-1999 1.3   LRB    add comments
C  31-Dec-1999       SEL    changed fmt to I2.2 for year
C  12-JAN-2000 1.4   LRB    changed fmt 100 to I2.2 for year
C   7-MAR-2000       JAP    Temporary fix for high speeds
C-
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CDEC$ IDENT '1.5'
C
      IMPLICIT NONE
      INTEGER UNIT, WHO_CALLED
      LOGICAL OUTPUT
      REAL*8 AVERAGE_EPOCH
C
      INTEGER I, J, FLAG, PM_QUALITY, AVGCNTR
      REAL V_RATIO, FIT_OVER_MOM
      REAL DENSITY_COM, VSW_COM, VTH_COM, ANGLE
      REAL VSW, VTH, MSW, ALPHA_TO_PROTON, DENSITY, CHI2, VSW_CHI_MIN, VTH_CHI_MIN
      REAL SUMW, SUMWX, SUMWX2, THETA, SIGMA_THETA, W(0:19)
      REAL VSW_UNC, VTH_UNC, T96, TVAC, CF
      REAL*8 VSW_NOTUSED, VTH_NOTUSED, DENSITY_NOTUSED, DENSITY_MOM
      REAL*8 VSW_MOM, VTH_MOM
      REAL*8 MOMENTS(2)
      REAL*8 TIME_6MAR/1204837800.0/, TIME_10JUN/1213124076.0/, TIME_18APR98/1271548800.0/
      STATIC PM_QUALITY
C
C The following variables and the COMMON block SUM hold the accumulated sums
C
      INTEGER*4 PM_RATE(6), YEAR, DOY, HR, MIN, SEC, MINRATE, MAXRATE
      INTEGER*4 MMCOUNTER, UNC_PM_RATE(6)
      REAL*4 PM_EQ(6), PM_COUNTER(6), PM_T1(120), TCOUNT(120)
      LOGICAL PRINT
      COMMON /SUM/ PM_EQ, PM_RATE, YEAR, DOY, HR, MIN, SEC, PM_COUNTER, PRINT,
     &             MINRATE, MAXRATE, MMCOUNTER, PM_T1, TCOUNT, UNC_PM_RATE
C
C Include the definition of the format statement to use
C
      INCLUDE 'INC_DIR:FORM.INC'
      DATA FMT(1:39)  /'(I2.2,1X,I3,1H:,I2.2,1H:,I2.2,1H:,I2.2,'/
      DATA FMT(40:93) /'F7.1,F7.2,F7.2,F7.2,2F6.2,2F7.1,2F7.2,2F7.2,F5.1,F6.2,'/
      DATA FMT(94:)   /'F7.2,I6,I8,2F6.1,2F7.2,I6)'/
      DATA NPARAMS/22/
      DATA NUM_REAL /15, 4, 8*0/
      DATA NUM_INT  /2, 1, 8*0/
      DATA ANGLE_INDEX  /4/
C
C Start code, test if we should compute and print average or if we should just
C reset the counters
C
      IF (OUTPUT) THEN
C
C Divide the sums by the number of readouts accumulated over, giving averages,
C PROVIDED that we got at least 1 value!
C
        DO 80  I = 1, 6
C
C If we got no values for this PM_EQ, set it to a fill value. Otherwise...
C
          IF(PM_COUNTER(I) .EQ. 0)THEN
            PM_EQ(I) = -1.0
C
C ...divide the sums by the number of readouts accumulated over, giving
C averages.
C
          ELSE
            PM_EQ(I) = PM_EQ(I)/PM_COUNTER(I)
          ENDIF
  80    CONTINUE
C
        DO 90  I=1, 6
C
C If we got no values for this RATE, set it to a fill value. Otherwise...
C
          IF(PM_COUNTER(I) .EQ. 0)THEN
            UNC_PM_RATE(I) = -1
            PM_RATE(I) = -1
            PM_QUALITY = -1
C
C ...divide the sums by the number of readouts accumulated over, giving
C averages.  
C
          ELSE
            PM_RATE(I) = PM_RATE(I)/PM_COUNTER(I)
            UNC_PM_RATE(I) = UNC_PM_RATE(I)/PM_COUNTER(I)
          ENDIF
  90    CONTINUE
C
C If we got no values for the maxima or minima set them to a fill value.
C Otherwise... 
C
        IF(MMCOUNTER .EQ. 0)THEN
          MAXRATE = -1
          MINRATE = -1
C
C ...divide the sums by the number of readouts accumulated over, giving
C averages.  
C
        ELSE
          MAXRATE = MAXRATE/MMCOUNTER
          MINRATE = MINRATE/MMCOUNTER
        ENDIF    ! MMCOUNTER .EQ. 0
C
C If WHO_CALLED = 1, then this routine was ultimately called by PMRATES.
C Otherwise it was called by  PMSW
C
        IF(WHO_CALLED .EQ. 1)THEN
C
C Write out the data we have just finished averaging, together with the times
C from the first readout in the averaging interval.  Don't write out any lines
C that have only zeros or -1's (missing), though.
C
          PRINT=.FALSE.
          DO 95 I = 1, 6
            IF(PM_EQ(I) .GT. 0) PRINT=.TRUE.
  95      CONTINUE
          IF(PRINT)WRITE(UNIT,100) YEAR, DOY, HR, MIN, SEC, (UNC_PM_RATE(I),I=1,6), (PM_EQ(I), I=1,6)
C
 100      FORMAT(I2.2,1X,I3,':',I2.2,':',I2.2,':',I2.2,3X,6(1X,I7),3X,6(1X,F5.2))
        ELSE
C
C Since we plan to use the theta distribution as weighting factors for the
C average 
C
          DO 120  J = 1, 19
            DO 110 I = 1, 6
              IF(TCOUNT(I) .NE. 0)THEN
                W(J) = W(J) + PM_T1(J+(I-1)*20 + 1)
                SUMW = SUMW + PM_T1(J+(I-1)*20 + 1)
              ELSE
                W(I) = -1
C
C We don't have a complete set; eject from the loop
C
                GOTO 140
              ENDIF
 110        CONTINUE
 120      CONTINUE
C
C Now have 20 #'s, summed over 6 step values. Find weighted average and sigma
C
          DO 130 I = 1, 19
            SUMWX = SUMWX + W(I)*I
            SUMWX2 = SUMWX2 + W(I)*I*I
 130      CONTINUE
C
          IF(SUMW.NE.0)THEN 
            THETA = SUMWX/SUMW
            SIGMA_THETA = SQRT((SUMWX2/SUMW - THETA**2))
          ELSE
            THETA = -1.0
            SIGMA_THETA = -1.0
          ENDIF
C
C Don't attempt to fit if any of the rates are fill data
C
 140      IF(PM_QUALITY .EQ. 0)THEN
C
C Correct rates for background
C
            CALL BACKGROUND_C(PM_RATE,FLAG)
C
C Find SW parameters
C
            CALL SWPARAMETERS(PM_RATE,VSW,MSW,ALPHA_TO_PROTON,DENSITY,VTH,CHI2,AVGCNTR,VSW_CHI_MIN,VTH_CHI_MIN)
C
C Save uncorrected (but Chi sq averaged) Vsw
C
            VSW_UNC = VSW
            VTH_UNC = VTH
C
C Call modified Hefti routine(some improvements, but unweighted rates).  This
C is used for the solar wind speed from a moment calculation (VSW_MOM) only.
C
	    CALL HEFFALUMP_H(PM_RATE, AVERAGE_EPOCH, VSW_MOM, VTH_NOTUSED, DENSITY_NOTUSED)
C
C Line segment speed correction
C
            CALL CORRECT_VSW_H(VSW_MOM)
C
C Call Fred Ipavich/Simon Hefti routine (rates weighted with V**-2, other
C changes) to get the thermal speed and the density from moment calculations
C (VTH_MOM, DENSITY_MOM). 
C
	    CALL FISH(PM_RATE, AVERAGE_EPOCH, VSW_NOTUSED, VTH_MOM, DENSITY_MOM, MOMENTS)
C
C Correction to DENSITY_MOM based on velocity from HEFFALUMP_H
C
            DENSITY_MOM = DENSITY_MOM / ( 27.68 * VSW_MOM**(-0.55339) )
C
C Renormalize density if time is before Mar 6 1996 at 2110, when bias changed
C (PM MCP bias=level 150 before 6 MAR 96, 155 thereafter). Since the bias
C increased, and we are normalizing the early data to the data after Mar 6, the
C density from before Mar 6 must be increased 
C
            IF(AVERAGE_EPOCH .LT. TIME_6MAR)DENSITY = DENSITY * 1.28897
C
C On Jun 10, the threshold was changed. Since we are normalizing the data after
C Jun 10 to the data in the interval from Mar 6 to Jun 10, after the threshold
C change, density after Jun 10 must be increased.
C
            IF(AVERAGE_EPOCH .GE. TIME_10JUN)DENSITY = DENSITY * 0.620065
C
C On Apr 16, 1998 MTOF was turned off due to loss of Vf referenced voltages and
C a high main current.  It remained off for all of the 17th and most of the
C 18th.  When the PM was turned back on, the MCP bias was increased to 160.  As
C a result, desities after that date must be decreased. 
C
            IF(AVERAGE_EPOCH .GE. TIME_18APR98)DENSITY = DENSITY / 1.9
C
C By subtracting the epoch time of Jan 1 00:00:00 from the current epoch time,
C get seconds since the start of 1996
C
            T96 = AVERAGE_EPOCH - 1199145600.0
C
C Before the "SOHO vacation" we have a time decay represented by the double
C exponetial below.  After the vacation there is a recovery, and we need a
C different equation to describe the time decay.  MTOF was turned on for the
C first time after the "vacation" on 1998 DOY 301, so DOY 302 is used as the
C demarcation point.  
C
            IF(AVERAGE_EPOCH  .LT. 1288310400.0)THEN  ! 1288310400.0 = epoch time of 1998 302 0000
C
C Remove decay in density caused by loss of MCP gain.  At first, the decay was
C exponential, with a relatively short time constant.  Then the rate of decay
C slowed.  As a result the function chosen to fit the decay was of the form
C A*EXP(-T/B) + C*EXP(-T/D) + E, where the additive contant is inserted to
C guarantee that the MCP gain does not reach zero. 
C
              DENSITY = DENSITY / (1.6331*EXP(-T96/5.1858E6) + 0.54036*EXP(-T96/8.4945E7) + 0.05)
            ELSE
C
C Compute the time in seconds since the end of the "vacation"; i.e. since 1998
C 302 0000.  The density correction is an exponential falloff beginning at that
C time, plus a constant to prevent gain from ever reaching zero.
C
              TVAC = AVERAGE_EPOCH - 1288310400.0
              DENSITY = DENSITY / (0.40954*EXP(-TVAC/5.9519E7) + 0.03)
            ENDIF
C
C Correct Vth, based on Vsw
C
            CALL CORRECT_VTH(VTH,VSW)
C
C Correct Vsw (this is no longer a rate correction; its based on speed itself)
C
            CALL CORRECT_VSW(VSW)
C
C Correct to remove residual density dependence on Max PM rate
C
            IF((MAXRATE .GE. 100).AND.(MAXRATE .LE. 750000))THEN
              DENSITY = DENSITY / (1.0127*EXP(-7.874E-7*FLOAT(MAXRATE)))
            ELSE
              PRINT = .FALSE.
            ENDIF
C
C More density corrections
C
            CF=0.7727+2.9154*(VTH/VSW)
            IF(CF .GT. 1.2)CF=1.2
            DENSITY = DENSITY /CF
C
C Remove residual dependances; do this last to avoid messing up earlier
C corrections 
C
            VTH = VTH/(0.88147+3.7559E-3*VTH)
            VTH = VTH/(1.0788-1.8905E-4*VSW)
C
C Correction for BOTH thermal velocities!
C
            IF(CHI2 .LT. 13.)THEN
              CF = 1.0
            ELSEIF(CHI2 .LT. 80.)THEN
              CF = 0.62 * (CHI2)**0.185
            ELSE
              CF = 1.4
            ENDIF
            VTH = VTH / CF
            VTH_MOM = VTH_MOM / CF
C
C Second line segment corrections.  While both corrections could, in principle,
C be handled by one routine, the creation of a second routine is easier than
C returning to uncorrected data and deriving a new version of CORRECT_VSW or
C CORRECT_VSW_H. 
C
            CALL CORRECT_VSW_JUN97(VSW)
            CALL CORRECT_VSW_H_JUN97(VSW_MOM)
C
C Correct density to reflect new speed corrections.  Limit the correction to
C the range of data used to derive it.  This correction replaces two previous
C ones. 
C
            V_RATIO = VSW_UNC/VSW
            IF(V_RATIO .LT. 0.85) V_RATIO = 0.85
            IF(V_RATIO .GT. 1.015) V_RATIO = 1.015
            CF = 27.296 - 54.411*V_RATIO + 28.104*V_RATIO**2
            DENSITY = DENSITY / CF
C
C Combine results of 2 methods.  Restrict correction factor to a reasonable
C value by restricting the range of the ratio of fit to mom
C
            VSW_COM = -1.0
            FIT_OVER_MOM = VSW/VSW_MOM
            IF(FIT_OVER_MOM .LT. 0.85)FIT_OVER_MOM = 0.85
            IF(FIT_OVER_MOM .GT. 1.1)FIT_OVER_MOM = 1.1 
            IF(FIT_OVER_MOM .LE. 1.0)THEN
              IF(VSW_MOM .NE. 0.0)VSW_COM = VSW_MOM/(0.73285 + 0.27034 * (FIT_OVER_MOM))
            ELSE
              IF(VSW_MOM .NE. 0.0)VSW_COM = VSW_MOM/(1.0114 - 0.0081712 * (FIT_OVER_MOM))
            ENDIF 
C
C Temporary fix, because our combined speed can't go high enough when the speed as high (based on 
C comparison to SWE)
C
            IF((VSW_COM .GT. 600.0).AND.(VSW_COM .LE. 750.))THEN
              VSW_COM = VSW_COM + (VSW_COM - 600.)/3.
            ELSEIF(VSW_COM .GT. 750.)THEN
              VSW_COM = VSW_COM + 50.
            ENDIF
C
            DENSITY_COM = -1.0
            FIT_OVER_MOM = DENSITY/DENSITY_MOM
            IF(FIT_OVER_MOM .LT. 0.8)FIT_OVER_MOM = 0.8 
            IF(FIT_OVER_MOM .GT. 1.6)FIT_OVER_MOM = 1.6
            IF(DENSITY_MOM .NE. 0.0)DENSITY_COM = DENSITY/(0.63284 + 0.35416 * (FIT_OVER_MOM))
C
            VTH_COM = -1.0
            FIT_OVER_MOM = VTH/VTH_MOM
            IF(FIT_OVER_MOM .LT. 0.8)FIT_OVER_MOM = 0.8 
            IF(FIT_OVER_MOM .GT. 1.2)FIT_OVER_MOM = 1.2 
            IF(VTH_MOM .NE. 0.0)VTH_COM = VTH_MOM/(1.2251 - 0.2303 * (FIT_OVER_MOM)) 
C
C Derive angle from theta bin
C
            CALL CORRECT_ANGLE(THETA, ANGLE)
C
            IF(PRINT)THEN
              WRITE(UNIT,FMT(1:FMTLEN)) YEAR, DOY, HR, MIN, SEC, VSW_COM,
     &             DENSITY_COM, VTH_COM, ANGLE, THETA, SIGMA_THETA, VSW,
     &             VSW_MOM, DENSITY, DENSITY_MOM, VTH, VTH_MOM, MSW,
     &             ALPHA_TO_PROTON, CHI2, MINRATE, MAXRATE, VSW_CHI_MIN,
     &             VSW_UNC, VTH_CHI_MIN, VTH_UNC, AVGCNTR 
            ELSE
              CALL WRITE_FLAG_RECORD(UNIT, YEAR, DOY, HR, MIN, SEC, -1.)
            ENDIF
          ENDIF    ! PM_QUALITY .EQ. 0
        ENDIF    ! WHO_CALLED .EQ. 1
      ENDIF    ! OUTPUT
C
C Clear PM_RATE arrays and PM_COUNTER
C
      DO 210 I = 1, 6
        PM_EQ(I) = 0.0
        PM_RATE(I) = 0
        UNC_PM_RATE(I) = 0
        PM_COUNTER(I) = 0.0
        DO 200 J=1,20
          PM_T1(J+20*(I-1)) = 0.0
          TCOUNT(J+20*(I-1)) = 0.0
 200    CONTINUE
 210  CONTINUE
C
C Reset other variables
C
      MAXRATE = 0
      MINRATE = 0
      MMCOUNTER = 0
C
      DO 215 J=1,19
        W(J) = 0.0
 215  CONTINUE
      SUMW = 0.0
      SUMWX = 0.0
      SUMWX2 = 0.0
C
      PM_QUALITY = 0
C
      RETURN
      END
