CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        PROGRAM PM_SW
C+
C 
C ABSTRACT: 
C 
C     This program reads one or more CELIAS science data files, finds solar
C     wind parameters from the PM rate and E/Q step data from each science
C     record, and writes these parameters out in a column format suitable  
C     for use by Kaleidagraph.
C 
C USAGE:
C
C     RUN PMSW
C
C INPUTS:
C
C     Start time          (user supplied)
C     Stop time           (user supplied)
C     Averaging interval  (user supplied)
C     CELIAS science data file(s)
C
C OUTPUTS:
C
C     Disk file containing solar wind parameters; default filename is PMSW.DAT.
C
C COMMON BLOCKS:
C  
C     AVERAGE		contains averaging interval AVINT
C     PACKET_CONTROL	contains NPACKETS, PACKET_COUNTER
C     FORM		contains format control information
C
C INCLUDE FILES:
C
C     EDBDEF.INC
C     INPUTDEF.INC
C     FORM.INC		defines FORM common block
C
C DESIGN:
C  
C     {tbs}
C  
C This version:   4.3   29-MAR-2000
C 
C CREATED: 
C 
C     17-MAY-1996    John Paquette
C 
C DEVELOPMENT HISTORY:
C 
C  Author       Date		Description of Change 
C  --------	-----------	--------------------- 
C  J. PAQUETTE  17-MAY-1996     Created (from PM_RATES)
C  J. PAQUETTE  29-MAY-1996     Add capability to span files
C  L. Bleau	 5-JUN-1996	Call READ_ORBIT_INFO once per day (LZ file)
C  J. PAQUETTE  10-JUN-1996     Change header
C  J. PAQUETTE  13-JUN-1996     Change header
C  J. PAQUETTE  19-JUN-1996     Add code to span from LZ to REL files
C  J. PAQUETTE  17-JUL-1996     Undo Larry's format change of 12 July
C  J. PAQUETTE  22-JUL-1996     Deal with failure to open input file
C  J. PAQUETTE  24-JUL-1996     In the event that the SR start time is 0, use
C				the EDB(0) time
C  J. PAQUETTE  15-AUG-1996     When an LZ file is missing, use REL/QKL, but
C                               switch back to LZ's as soon as possible
C  J. PAQUETTE  15-AUG-1996     Place repetitive section of code into a new
C                               subroutine, CONSTRUCT_TOMORROWS_FILENAME
C  J. PAQUETTE  16-AUG-1996     Place the code that allows processing to go
C                               back to LZ's in a subroutine CHECK_MIDNIGHT
C  J. PAQUETTE  19-AUG-1996     Correct logic identifying REL/QKL files; move
C                               assignment to JUST_OPENED to a less misleading
C                               place
C  J. PAQUETTE  21-AUG-1996     Add STATUS to call to CONSTRUCT_FILENAME; STOP
C                               if its bad 
C  J. PAQUETTE  29-AUG-1996     Check for start time based on EDB time
C  J. PAQUETTE   5-SEP-1996     Save last good time; use that (less 1/2 hr) as
C                               argument for CONSTRUCT_TOMORROWS_FILENAME
C  J. PAQUETTE   2-OCT-1996     Add code to section that checks if midnight was
C                               spanned to deal with the case when file break
C                               is within 1 SR of midnight 
C  J. PAQUETTE   7-JAN-1997     Change header
C  J. PAQUETTE  15-JAN-1997     Change header again
C  J. PAQUETTE  17-JAN-1997     Change header again
C  J. PAQUETTE  22-JAN-1997     Change header again
C  L. BLEAU     22-JAN-1997     add COMMON block FORM containing runtime format;
C                               add code to write it out in header
C  J. PAQUETTE  29-JAN-1997     Changeover from PMSW to PMSW_A
C  J. PAQUETTE   3-FEB-1997     Major header change to accomadate 2 new columns
C  L. BLEAU      3-FEB-1997     add call to MOD_FORMAT to modify runtime format
C				to use Fw.0 in place of Iw format specifiers
C  J. PAQUETTE   5-FEB-1997     Major header change to accomadate 2 new columns
C  L. BLEAU      5-FEB-1997 2.4 remove common FORM declaraction, add INCLUDE
C                               FORM.INC 
C  J. PAQUETTE   6-FEB-1997     Add definition of N/S to output filer header
C  J. PAQUETTE   6-FEB-1997     Minor header change
C  L. BLEAU      6-FEB-1997 2.7 add call to SIGNON
C  J. PAQUETTE  13-FEB-1997     Header change, new reccomended print command
C                               for output file 
C  J. PAQUETTE  13-FEB-1997     Header change
C  J. PAQUETTE  12-MAR-1997     Header change
C  J. PAQUETTE  15-APR-1997     Change header
C  J. PAQUETTE  15-MAY-1997     Deal with case when first file cannot be opened
C                               (e.g., if the first file is being transferred
C                               when the program is running) 
C  J. PAQUETTE  15-MAY-1997     Change header
C  J. PAQUETTE  24-JUL-1997     Change back PMSW
C  J. PAQUETTE  20-NOV-1997     PMSWC
C  J. PAQUETTE   4-DEC-1997     Renamed from PMSWC to PMSW
C  J. PAQUETTE   6-MAY-1998     Change header. Again.
C  J. PAQUETTE   8-JAN-1999     Replace call to CONSTRUCT_TOMORROWS_FILENAME
C                               with a call to CONSTRUCT_NEXT_AVAIL_FILENAME
C  J. PAQUETTE   2-FEB-1999     Print out last partial average even when 
C                               CONSTRUCT_NEXT_AVAIL_FILENAME fails
C  J. PAQUETTE  19-FEB-1999     Call CONSTRUCT_NEXT_AVAILABLE_FILENAME even if
C                               program was not processing LZ's
C  L. Bleau      5-MAR-1999 3.7 add IF-THEN-ELSE to code after end of DO 100
C				loop to check if we should span files or not;
C				if we should, check for next file, which is
C				what the code now does; if not, add call to
C				TABLE_PM_SR to calculate final partial average
C  J. PAQUETTE  10-MAR-1999 3.8 If we call CONTRUCT_NEXT_AVAIL_FILENAME and we
C  L. BLEAU                     were processing LZ files, subtract 1800. sec
C                               from LAST_GOOD_SR_EPOCH as before; if we were
C                               processing REL/QKL files don't subtract
C                               anything!  This change avoids an infinite loop
C                               that sometimes occurs when only the first REL
C                               file of a new day has come in. 
C  L. Bleau     11-MAR-1999     update, reformat comments
C  L. Bleau     12-MAR-1999 3.9	added line to reset FIRST under proper conditions
C  J. PAQUETTE  13-APR-1999     Changed WRITE(7,10) to TYPE 10 in error message
C  L. Bleau	20-APR-1999	added IMPLICIT NONE, declarations
C  J. PAQUETTE   3-MAY-1999     Add more description of PM Max and Min to header
C  L. Bleau	14-JUN-1999	add code to check for backward time step, don't
C				overwrite LAST_GOOD_SR_EPOCH
C  J. PAQUETTE  17-JUN-1999     Use EDB.QUALITY = -2 for EDB's before the start time
C  L. Bleau	21-JUN-1999 4.0	remove file selection code (end of loop 50) and
C				replace it with a call to SELECT_FILE
C  L. Bleau	28-JUN-1999 4.1	add argument to SELECT_FILE call
C  L. Bleau	25-AUG-1999	add diagnostic WRITEs as D lines
C  L. Bleau	 8-NOV-1999	add comments
C  J. PAQUETTE   7-MAR-2000     Temporary header change, describing Vsw kluge
C  J. PAQUETTE  29-MAR-2000     Header change
C-

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CDEC$	IDENT	'4.3'
        INCLUDE 'INC_DIR:EDBDEF.INC'
        INCLUDE 'INC_DIR:INPUTDEF.INC'
        BYTE APID
        INTEGER*4 STATUS, HDR_LEN, I, J, AVINT, YEAR, DOY, HOUR, MINUTE, SECOND
        INTEGER*4 MSEC, MICRO, MONTH, DAY, KK
        INTEGER*4 NPACKETS, PACKET_COUNTER, QAC_OFFSET, CONTROL, NUM_SR
        REAL*8 NEWFILE_EPOCH, OLD_SR_EPOCH, LAST_GOOD_SR_EPOCH
        LOGICAL DONE/.FALSE./, FOUND_START/.FALSE./, FIRST/.TRUE./, OPENFILE
        LOGICAL JUST_OPENED, SPAN_MIDNIGHT/.FALSE./, NEXT_ONE/.FALSE./
	LOGICAL QUIT, PART_AVG
        CHARACTER*80 NEW_NAME
        CHARACTER*512 REL_FILE_HEADER
        RECORD /FILE_HDR/ FILE_HEADER
        RECORD /SCI_REC/ SR, EMPTY
        COMMON /AVERAGE/ AVINT
        COMMON /PACKET_CONTROL/ NPACKETS, PACKET_COUNTER
	RECORD /INPUTS/ USER
	INCLUDE 'INC_DIR:FORM.INC'
	INTEGER*4 FMTLEN2
	CHARACTER*400 FMT2
	DATA FMTLEN/0/
C
C Prompt user for all information needed to run this program
C
	CALL SIGNON
	CONTROL = X'110'
	USER.READOUTS_PER_MIN = 2.
	CALL SOHO_GET_INPUTS(CONTROL, 29, 'PMSW.DAT', 'PMSW', USER)
	AVINT = USER.NREADOUTS
C
C Check to make sure that input time is after January 20   
C
        IF(USER.START_EPOCH .LT. 1200850200.)THEN
          TYPE 10
  10      FORMAT(' Start time cannot be before Jan 20 17:30.')
          STOP
        ENDIF
C
C To deal with the rather arbitrary record boundaries (e.g. on CD the records
C are 32767 bytes long), open using a C routine instead of FORTRAN OPEN
C
  50    CONTINUE
D	WRITE(88,1003) USER.INPUT_NAME(1:USER.I_NAME_LEN)
D1003	FORMAT(' About to open file ',A)
	CALL OPEN_STREAM_IN (1, USER.INPUT_NAME(1:USER.I_NAME_LEN), STATUS)
        JUST_OPENED = .TRUE.
	IF (STATUS .NE. 1) THEN
	  TYPE *,'ERROR OPENING FILE ',USER.INPUT_NAME(1:USER.I_NAME_LEN)
	  TYPE *,'Attempting to proceed to next file.'
          OPENFILE = .FALSE.
D	WRITE(88,*) 'File open failed'
        ELSE
          OPENFILE = .TRUE.
	ENDIF
C
        IF(OPENFILE)THEN
C
C Call routine to access orbit information and read it into memory
C
c	CALL READ_ORBIT_INFO (USER.INPUT_NAME(1:USER.I_NAME_LEN))
C
C Open Output file.
C
          OPEN(UNIT=15, NAME=USER.OUTPUT_NAME, FORM='FORMATTED', CARRIAGECONTROL='LIST', TYPE='NEW', RECL=187)
C
C Write program name, run date, and input file name to file.
C
          IF(FIRST)CALL PRINT_HEADER(15, USER)
C
C To find the length of the header (which is a constant 48 bytes in CD files,
C but is variable in REL files) look for the start of the first packet.  A
C packet begins with 88 (APID), where APID is the APplication ID for the
C relevant type of packet.  For CELIAS science this is A9(hex) = 169(decimal).
C
          APID = 'A9'X
          CALL FIND_HEADER_LENGTH(USER.INPUT_NAME(1:USER.I_NAME_LEN), APID, HDR_LEN)
C
C The header length tells us what kind of file that we have.  If its 48 bytes,
C we have data file from the CD.  Get its (binary) file header, output it to
C file, and use the number of packets obtained from the header to calculate the
C number of SR in the file. 
C
          IF(HDR_LEN .EQ. 48)THEN
            CALL GET_FILE_HEADER(1, FILE_HEADER)
            IF(FIRST)CALL DUMP_FILE_HEADER(15, FILE_HEADER)
            NPACKETS = FILE_HEADER.NUM_PACKETS
            QAC_OFFSET = FILE_HEADER.QAC_OFFSET
C
C Otherwise, we have an REL file. Dump its (ASCII) header, extract the number
C of packets, and calculate the number of SR in the file.
C
          ELSE
            IF(FIRST)THEN
              CALL DUMP_REL_FILE_HEADER(1, 15, HDR_LEN, NPACKETS, QAC_OFFSET)
            ELSE
              CALL GET_REL_FILE_HEADER(1, HDR_LEN, NPACKETS, QAC_OFFSET)
            ENDIF
          ENDIF
C
C Either way, the number of science records is equal to the following
C
          NUM_SR = NPACKETS/(8*20)+1
C
C Open the input file again for access to the QAC information; skip through it
C until we get the first QAC list element.
C
          CALL FIND_QAC(QAC_OFFSET, 394, HDR_LEN, USER.INPUT_NAME(1:USER.I_NAME_LEN))
C
C Print header (first time only)
C Also write out number of parameters (columns of values) and the format used
C to write them.  Both of these come from COMMON block FORM, which is
C intialized at runtime to the proper values.
C
          IF(FIRST)THEN
	    FMTLEN = INDEX(FMT,')   ')
	    IF(FMTLEN .GT. 0) THEN
	      CALL MOD_FORMAT(FMT(1:FMTLEN),FMT2,FMTLEN2)
	      WRITE(15,89) NPARAMS,FMT2(1:FMTLEN2)
	    ENDIF
  89	    FORMAT(' NPARAMS=',I2.2/
     &	           ' FORMAT=',A)
	    WRITE(15,90)
  90        FORMAT('Vsw2 = solar wind speed corresponding to minimum Chi sq'/,
     &             'Vsw3 = Chi sq averaged Vsw, but without speed correction'/,
     &             'Vth2 = thermal speed corresponding to minimum Chi sq'/,
     &             'Vth3 = Chi sq averaged Vth, but without speed correction or power law Vth correction'/,
     &             'ANGLE = arrival direction in degrees (generally N/S, with + meaning FROM the south)'/,
     &             'N/S = average theta bin after excluding bin 0'/,
     &             'PM_Min and PM_Max are corrected for quantization error but not for deadtime'/,
     &             'PM_Max is the average of maxima, not the maximum of averaged rates. PM_Min is the average of minima.'/,
     &             'Uses 96% background correction (really); then 80%, 50%, ...'/
     &             'As of 8 Mar 2000: temporary kluge to increase speeds above 600 km/s'//,
     &             'YR DOY:HR:MN:S    Vsw   Dens    Vth   ANGLE   N/S   Sig  Vsw_f  Vsw_m  Dens_f Dens_m Vth_f  Vth_m'
     &             '  Msw Alph%  CHI2  PM_Min PM_Max Vsw2  Vsw3    Vth2  Vth3  #Poss')
	  ENDIF
C
C Fill and then dump each science record. The SR structure is zeroed after each
C dump to avoid problems when there's missing data.
C
          PACKET_COUNTER = 0
          I = 0
          DO 100 WHILE((PACKET_COUNTER .LE. NPACKETS) .AND. (.NOT. DONE) .AND. (.NOT. SPAN_MIDNIGHT))
            CALL GET_SR(1, SR)
C
C Check to see if the start time of the science record that we just got is
C within 300 seconds of the start time.
C
            IF(SR.EPOCH .GE. (USER.START_EPOCH-300.D0))THEN
C
C If it is, then step through the EDB's, comparing the time of each to the
C start time. 
C
              J = 0
              DO 95 WHILE((J .LT. 20).AND.(.NOT.FOUND_START))
C
C If the EDB time is greater than or equal to the start time, set the
C FOUND_START flag. If not, mark the quality flag for the EDB as bad, so that
C it will not contribute to the average. 
C
                IF(SR.EDB(J).EPOCH .GE. USER.START_EPOCH)THEN
                  FOUND_START = .TRUE.
                ELSE
                  SR.EDB(J).QUALITY = -2
                ENDIF
                J = J + 1
  95          CONTINUE
            ENDIF
C
C If we found a time past the start time, then assume that all subsequent times
C are past the start time also. 
C
            IF(FOUND_START)THEN
D	WRITE(88,1001) SR.EPOCH
D1001	FORMAT(' SR.EPOCH = ',T25,F15.0)
              CALL TABLE_PM_SR(15, SR, DONE, USER.STOP_EPOCH, 2)
C
C If we have an REL or QKL file, check if we just spanned midnight
C
              IF(((USER.INPUT_TYPE .EQ. 2) .OR. (USER.INPUT_TYPE .EQ. 3)) .AND. (.NOT.JUST_OPENED))THEN 
C
C When the start of a new file is within 1 SR of midnight, use the time of the
C first good SR in place of the SR start time. 
C
                IF(NEXT_ONE)THEN
                  NEXT_ONE =.FALSE.
                  DO 97 KK = 0, 19
                    IF(SR.EDB(KK).QUALITY .EQ. 0)OLD_SR_EPOCH = SR.EDB(KK).EPOCH
   97             CONTINUE
                ENDIF
C
C Check if this science record is in a different day (crossed midnight) and we
C should switch to using the LZ file for that day.
C
                CALL CHECK_MIDNIGHT(OLD_SR_EPOCH, SR.EPOCH, SR.EDB(0).EPOCH,
     &                              USER.INPUT_NAME(1:USER.I_NAME_LEN), 
     &                              USER.I_NAME_LEN, USER.INPUT_TYPE, 29,
     &                              SPAN_MIDNIGHT) 
              ENDIF
            ENDIF
C
C Save the Epoch time of this science record in OLD_SR_EPOCH unless said Epoch
C time is zero.  In that instance, access the Epoch time of EDB 0 and save it
C in OLD_SR_EPOCH instead.
C Also save the Epoch time of this science record in LAST_GOOD_SR_EPOCH unless
C the time stepped backwards by more than a day.  If that happens, leave 
C LAST_GOOD_SR_EPOCH alone.
C
            IF(SR.EPOCH .EQ. 0)THEN
              OLD_SR_EPOCH = SR.EDB(0).EPOCH
            ELSE
	      IF (OLD_SR_EPOCH - SR.EPOCH .LT. 86400.) THEN
                LAST_GOOD_SR_EPOCH = SR.EPOCH
	      ENDIF
              OLD_SR_EPOCH = SR.EPOCH
            ENDIF
C
C Clear SR structure, and increment counter
C
            SR = EMPTY          
            I = I + 1
            IF(JUST_OPENED)NEXT_ONE = .TRUE.
            JUST_OPENED = .FALSE.
 100      CONTINUE
C
        ENDIF   ! IF(OPENFILE)
D	WRITE(88,1002) LAST_GOOD_SR_EPOCH
D1002	FORMAT(' LAST_GOOD_SR_EPOCH = ',T25,F15.0)
C
C If we got out of the DO 100 loop, one of the following conditions is true:
C
C    Condition			Meaning
C 1) DONE is .TRUE.		we have detected a time that is beyond the user
C				specified stop time (USER.STOP_EPOCH), and
C				should wind down the program
C 2) PACKET_COUNTER > NPACKETS	we have simply run out of data in the current
C				file, and should obtain the name of the next
C				data file and process it
C 3) SPAN_MIDNIGHT is .TRUE.	the REL/QKL file we are reading contains
C				midnight, an LZ file exists on the next day,
C				and we should switch to it and stop processing
C				the current REL/QKL file
C
C Call routine SELECT_FILE to determine the next data file that should be
C processed, if any.  If QUIT is .TRUE. after SELECT_FILE returns, then stop
C processing.  If PART_AVG is .TRUE., call partial average code before exiting.
C SELECT_FILE sets both of these flags .FALSE. when it is called, so we needn't
C set them here.
C
	CALL SELECT_FILE(29, DONE, SPAN_MIDNIGHT, QUIT, PART_AVG, FIRST,
     &			OPENFILE, USER, LAST_GOOD_SR_EPOCH, NEWFILE_EPOCH)
C
C If QUIT is .FALSE., loop back to 50 to process the next file.  Its name has
C already been stored into components of the record USER.
C
	IF (.NOT. QUIT) GOTO 50
C
C We come here only if QUIT is .TRUE., which means we've either reached the
C stop time, reached the end of file in a one-file-only mode,  or there are no
C more data files that can be read. 
C Under these conditions we don't want to proceed to process the next data
C file, if it indeed exists.  Instead, calculate the partial average if the
C flag PART_AVG is set.
C
	IF (PART_AVG) THEN
	  CALL TABLE_PM_SR(15, SR, DONE, USER.STOP_EPOCH, 2)
	ENDIF
C
C All done. Send the user a message telling him how to print the output file
C
        CLOSE(15)
        TYPE 1000, USER.OUTPUT_NAME(1:USER.O_NAME_LEN), USER.OUTPUT_NAME(1:USER.O_NAME_LEN)
1000    FORMAT(' Output file is ',A/' Recomended print command is MP4 ',A)
        END
