* for disk operation
*********1*********2*********3*********4*********5*********6*********7**
*      PROGRAM DRIP
* MODIFICATION OF PROGRAM BY S. CHAKRABARTI, G. GIBSON,
* T. SASSEEN, D. COTTON AND R. GLADSTONE
* WRITTEN AUGUST 1989 BY GREG FRUTH
* READS ECOM-721 DATA FROM DISK FILE AND CRUNCHES IT FOR PLOTTING
* BY PROGRAM ZAP
* modification by Dan Cotton to make time and zenith angles variable
*  8/6/94
*********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT	INTEGER*2 (I-N)
      PARAMETER	(MAXLNS = 11, MAXBNS = 340,  maxizbns=36)
      REAL	WLINE(MAXLNS)
      INTEGER*2	WL1(MAXLNS),WH1(MAXLNS),
     &	WL(MAXLNS), WH(MAXLNS),
     &	WL2(MAXLNS),WH2(MAXLNS)
      REAL	SB1,SB,SB2,DN1,DN2
      REAL	DTA(124),SCB(300),SCE(300)
      INTEGER*2	ISPEC(1308)
      REAL	PTTIME(MAXBNS),PTGLON(MAXBNS,2),PTGLAT(MAXBNS,2),
     &	PTBLAT(MAXBNS,2),PTSZEN(MAXBNS,2),
     &	S(128),WAV(128),PCNT(MAXBNS,maxizbns,2),
     &	PBK1(MAXBNS,maxizbns,2),PLIN(MAXBNS,maxizbns,2),
     &  PBK2(MAXBNS,maxizbns,2),ptzen(MAXBNS,maxizbns,2),
     &  ptazim(maxbns,maxizbns,2),ptsaz(maxbns,maxizbns,2)
      REAL	sqglon(MAXBNS,2),sqglat(MAXBNS,2),sqblat(MAXBNS,2),
     &	sqszen(MAXBNS,2),sqazim(maxbns,maxizbns,2),
     &  sqsaz(maxbns,maxizbns,2),sqzen(maxbns,maxizbns,2)
      REAL	sdglon(MAXBNS,2),sdglat(MAXBNS,2),sdblat(MAXBNS,2),
     &	sdszen(MAXBNS,2),sdazim(maxbns,maxizbns,2),
     &  sdsaz(maxbns,maxizbns,2),sdzen(maxbns,maxizbns,2)
      real      fcnt,ozen
      LOGICAL	L1,L2,L3,OK
      INTEGER*4	FLG,COUNT,GOOD,TGOOD,BAD,TBAD,INDEX
      CHARACTER	ATAP*8,RV*4,DY*4,YR*4,GRATFILE*50,JUNK*90,YN*1
      COMMON	/CORFAC/ FACTOR(128)
      COMMON	/DEVS/ TAPEFILE,LOG
      DATA	NSEC /60/
      DATA	NUMLNS,NUMBNS /0,0/
      DATA	IBL1,IBL2 /-1,-1/
      DATA	IRO,IGO /-1,-1/
      DATA	GOOD,TGOOD,BAD,TBAD,INDEX /0,0,0,0,1/
* info on how to read disk file
      character tname*18,base*10
      integer*4 nhed,nscan,nevt,nem
* LOGICAL UNIT ASSIGNMENTS:
      INTEGER	TAPEFILE,LOG
      DATA	TAPEFILE,LOG /9,10/
*********1*********2*********3*********4*********5*********6*********7**
* EXPLANATION OF SOME OF THE ARRAYS AND VARIABLES
*
* PCNT - WILL BE NUMBER OF SPECTRA IN EACH BIN
* PLIN - WILL CONTAIN COUNTS IN A LINE
* PBK1, PBK2 - CONTAIN COUNTS IN BACKGROUND RANGES, 1 AND 2
*
* ALL OF THESE ARRAYS CAN CONTAIN UP TO MAXBNS TIME BINS FOR UP TO
* MAXLNS LINES. THE COUNTS ARE SEPARATED INTO maxizbns zenith ANGLE BINS,
* THE FIRST zenith BIN GOING FROM 0 DEGREES UP TO BUT NOT INCLUDING
* zbnw DEGREES. THE REST ARE SIMILAR, EXCEPT THAT BIN maxizbns INCLUDES
* zenith ANGLE 180 DEGREES; inc is a flag which indicates whether the
* particular bin is during an increasing (1) or decreasing (2)
* zenith angle swing or.
* THE BIN WIDTH NSEC IS 60 SECONDS by default.
*
* FILES REQUIRED:
*	tapename - DISK FILE CONTAINING DATA ORIGNALLY STORED ON
*		   1/2 INCH TAPE (E.G. hula0375)
*	grat.n - FILE THAT CONTAINS THE WAVELENGTHS OF INTEREST
*		 THAT CAN BE SEEN WITH GRATING STEP N AND THE
*		 BIN NUMBERS FOR THE LINE WIDTH AND BACKGROUND
*		 REGIONS (E.G. grat.33)
*	caldta.dta - FILE THAT CONTAINS SENSITIVITIES VS.
*	             WAVELENGTHS FOR VARIOUS GRATING STEPS
* FILES PRODUCED:
*	tapename.two - PROCESSED LINE AND SPECTRUM DATA
*		       (E.G. hula0375.two)
*	tapename.log - PRINTOUT OF TAPE HEADER FILES, STATUS, ETC.
*		       (E.G. hula0375.log)
* NOTE: THE LOG FILE IS OPENED IN SUBROUTINE PRHEAD
*	caldta.dta is accessed only in subroutine convrt
*********1*********2*********3*********4*********5*********6*********7**
* ZERO ARRAYS
*********1*********2*********3*********4*********5*********6*********7**
      do inc=1,2
        DO ITBN = 1,MAXBNS
	  PTTIME(ITBN) = 0.0
	  PTGLON(ITBN,inc) = 0.0
          PTGLAT(ITBN,inc) = 0.0
	  PTBLAT(ITBN,inc) = 0.0
	  PTSZEN(ITBN,inc) = 0.0
	  sqglon(ITBN,inc) = 0.0
          sqglat(ITBN,inc) = 0.0
	  sqblat(ITBN,inc) = 0.0
	  sqszen(ITBN,inc) = 0.0
	  DO izbn = 1,maxizbns    
	    PCNT(ITBN,izbn,inc) = 0.0
            PBK1(ITBN,izbn,inc) = 0.0
            PLIN(ITBN,izbn,inc) = 0.0
            PBK2(ITBN,izbn,inc) = 0.0
	    ptzen(itbn,izbn,inc) = 0.0
	    sqzen(itbn,izbn,inc) = 0.0
	    ptazim(itbn,izbn,inc) = 0.0
	    sqazim(itbn,izbn,inc) = 0.0
	    ptsaz(itbn,izbn,inc) = 0.0
	    sqsaz(itbn,izbn,inc) = 0.0
          END DO
        END DO
      end do
      CALL GETCOR
*********1*********2*********3*********4*********5*********6*********7**
*
* READ HEADER RECORDS AND DATA
*
*********1*********2*********3*********4*********5*********6*********7**
      write(*,*) 'Enter the time bin in seconds (0 = 60 - default)'
      read(*,*) NSEC
      if(NSEC .le. 0) NSEC = 60
      print *, 'nsec = ',nsec
      write(*,*) 'Enter the zenith angle bin width (min width = 5)'
      read(*,*) izbnw
      nzbns = int(180./izbnw)
      if(nzbns .gt. maxizbns) then
        write(*,*) '***** Not enough zenith angle bins ****'
	write (*,'(''EDIT PROGRAM TO SET MAXBNS TO AT LEAST '',
     &                 I4,'' AND RECOMPILE'')') nzbns
	write(*,*) 'ABORTED.'
	stop
      endif
      write(*,*) izbnw
      write(*,*) 'Enter the name of the tape file to be processed'
      read(*,*) tname
      write(*,*) tname
c   10 OPEN (UNIT=TAPEFILE,FILE='/veebs2/scratch/bbush/ECOM/'//TNAME,
   10 OPEN (UNIT=TAPEFILE,FILE=TNAME,
     &      STATUS='OLD',FORM='UNFORMATTED')
      read(tapefile) nhed,nscan,nevt,nem
*********1*********2*********3*********4*********5*********6*********7**
* PRHEAD READS THE HEADER RECORD OF THE TAPE AND RETURNS THE TAPE NAME,
* THE REV #, THE START AND END TIMES, THE DAY OF THE YEAR AND THE YEAR.
*********1*********2*********3*********4*********5*********6*********7**
   20 CALL PRHEAD(nhed,ATAP,RV,TSTART,TEND,DY,YR)
      WRITE(*,*) 'THIS IS TAPE ',ATAP
      WRITE(*,*) 'REVOLUTION ',RV,'  DAY ',DY,'  YEAR ',YR
      WRITE(*,*) 'TAPE STARTS AT: ',TSTART
      WRITE(*,*) '       ENDS AT: ',TEND
      WRITE(*,*)
*********1*********2*********3*********4*********5*********6*********7**
* GET INPUT PARAMETERS
*********1*********2*********3*********4*********5*********6*********7**
      WRITE(*,*) 'ENTER START TIME FOR ANALYSIS',
     &			' (ENTER -1 FOR BEGINNING OF TAPE)'
      READ(*,*) TML
      WRITE(*,*) 'ENTER END TIME FOR ANALYSIS',
     &			' (ENTER -1 FOR END OF TAPE)'
      READ(*,*) TMH
      IF (TML .EQ. -1.0) TML = TSTART
      IF (TMH .EQ. -1.0) TMH = TEND
      WRITE(LOG,*) '*** ANALYSIS STARTS AT: ',TML,' SEC ***'
      WRITE(LOG,*) '***            ENDS AT: ',TMH,' SEC ***'
      WRITE(LOG,*)
      WRITE(*,*) '*** ANALYSIS STARTS AT: ',TML,' SEC ***'
      WRITE(*,*) '***            ENDS AT: ',TMH,' SEC ***'
      WRITE(*,*)
      NUMBNS = ( INT(TMH - TML) / NSEC) + 1
      IF (NUMBNS .GT. MAXBNS) THEN
	WRITE (*,*) '***** NOT ENOUGH TIME BINS! *****'
	WRITE (*,'(''EDIT PROGRAM TO SET MAXBNS TO AT LEAST '',
     &                 I4,'' AND RECOMPILE'')') NUMBNS
	WRITE(*,*) 'ABORTED.'
	STOP
      END IF
      WRITE(*,*) 'ENTER THE GRATING STEP (ENTER -1 FOR VALUE FROM TAPE)'
      READ(*,*) IG
      WRITE(*,*)
*********1*********2*********3*********4*********5*********6*********7**
* PRSCAN READS THE SCAN RECORD OF THE TAPE AND RETURNS ARRAYS SCB AND
* AND SCE, WHICH CONTAIN TIMES DURING WHICH THE DATA IS BAD
*********1*********2*********3*********4*********5*********6*********7**
   30 CALL PRSCAN(nscan,SCB,SCE)
*********1*********2*********3*********4*********5*********6*********7**
* PREVNT READS THE EVENT RECORD OF THE TAPE
*********1*********2*********3*********4*********5*********6*********7**
   40 CALL PREVNT(nevt)
*********1*********2*********3*********4*********5*********6*********7**
* SKIP EM DATA
*********1*********2*********3*********4*********5*********6*********7**
   50 do i=1,nem
	READ (TAPEFILE) JUNK
      end do
*********1*********2*********3*********4*********5*********6*********7**
* START READING DATA
*********1*********2*********3*********4*********5*********6*********7**
      ITBN = 0
      write(*,*) 'SEARCHING FOR TIME RANGE...'
   60 CONTINUE
      OK = .FALSE.
*********1*********2*********3*********4*********5*********6*********7**
* READDTA READS THE ASPECT PART OF THE DATA INTO DTA.
* THE FLAG THAT IS RETURNED IS:
*   0 - END OF FILE
* 540 - GOOD READ
*********1*********2*********3*********4*********5*********6*********7**
      CALL READDTA(DTA,FLG)
      IF (FLG .EQ. 0) GO TO 140
      IF (FLG .ne. 540) GO TO 130
      DO I = 1,124
        CALL IBMTOIEEE( DTA(I) )
      END DO
*********1*********2*********3*********4*********5*********6*********7**
* READISPEC READS THE SPECTRA PART OF THE DATA INTO ISPEC.
* THE FLAG THAT IS RETURNED IS:
*    0 - END OF FILE
* 2640 - GOOD READ
*********1*********2*********3*********4*********5*********6*********7**
   80 CALL READISPEC(ISPEC,FLG)
      IF (FLG .EQ. 0) GO TO 140
      IF (FLG .ne. 2640) GO TO 130
      CALL CRSPC(ISPEC)
*********1*********2*********3*********4*********5*********6*********7**
* SKIP DATA THAT IS NOT IN THE DESIRED TIME RANGE
*********1*********2*********3*********4*********5*********6*********7**
      IF (DTA(1) .GT. TMH) GO TO 140
      IF (DTA(1) .LT. TML) then
       zen=90.0-dta(124)
       GO TO 60
      end if
*********1*********2*********3*********4*********5*********6*********7**
* DTA CONTAINS ASPECT INFORMATION, ISPEC CONTAINS COUNT VALUES
* FOR TEN SPECTRA, FOR WHICH THERE ARE 128 CHANNELS EACH.
* SCB AND SCE CONTAIN TIME RANGES DURING WHICH DATA IS KNOWN TO BE BAD.
* THIS WILL SKIP BAD DATA.                                 
*********1*********2*********3*********4*********5*********6*********7**
   90 DO K= 1,300
        IF (DTA(1) .GT. SCB(K) .AND. DTA(1) .LT. SCE(K)) GO TO 130
        IF (DTA(1) .LT. SCE(K)) GO TO 100
      END DO
  100 CONTINUE
*********1*********2*********3*********4*********5*********6*********7**
* CONVERT WAVELENGTH AND FIND SENSITIVITY. NOTE THAT THE RSHIFT HAS
* BEEN ADDED. LSHIFT WAS THE ONLY FUNCTION AVAILABLE ON THE PERKEN ELMER
* AND COULD ACCEPT NEGATIVE SHIFTS. RSHIFT AND LSHIFT CAN'T TAKE
* NEGATIVE ARGUMENTS; THAT IS WHY THIS LOOKS REALLY STUPID.
* SEE ECOM-721 ENGINEERING REPORT FOR MEANING OF DSS FLAGS.
*********1*********2*********3*********4*********5*********6*********7**
      DSS1 = AND( (RSHIFT( ISPEC(1305),3 )),1 )
      DSS2 = AND( (RSHIFT( ISPEC(1306),6 )),1 )
      IF (DSS1 .EQ. DSS2) THEN
	WRITE(*,*) '** ECOM DSS ERROR! **'
	WRITE(*,*) '** ERROR IGNORED. **'
	WRITE(LOG,*) '>** ECOM DSS ERROR! **'
	WRITE(LOG,*) '>** ERROR IGNORED. **'
	GO TO 130
      END IF
*********1*********2*********3*********4*********5*********6*********7**
* IG2 CONTAINS GRATING STEP READ FROM TAPE. IG2 SOMETIMES GIVES 32,
* WHICH IS ASSUMED TO MEAN 33. SEE DOCUMENTATION.
* IR CONTAINS RESOLUTION MODE READ FROM TAPE.
* RESOLUTION MODE		RESOLUTION	DETECTOR
*	1			LOW		BOTH
*	2			HIGH		FILTERED
*	3			HIGH		UNFILTERED
*********1*********2*********3*********4*********5*********6*********7**
      IG2 = INT( 2.0 * FLOAT( ISPEC(1301) ) / 14.57143 + 1.07 )
      IF (IG2 .EQ. 32) IG2 = 33
      IF (IG .EQ. -1) IG = IG2
      IR = AND( RSHIFT( ISPEC(1306),3 ),3 )
      IF (IR .LT. 2) IR = 1
      L1 = (IG .EQ. IGO) .AND. (IR .EQ. IRO)
      L2 = (IBL1 .EQ. ISPEC(1305)) .AND. (IBL2 .EQ. ISPEC(1306))
      L3 = (IGO .EQ. -1) .AND. (IRO .EQ. -1)
      IF (L1 .AND. L2) GO TO 110
      IF (L3) THEN
	IF (IG .NE. IG2) THEN
	  WRITE(*,*) '** INPUT GRATING STEP AND INITAL GRATING STEP',
     &		' FROM TAPE DO NOT MATCH! **'
	  WRITE(*,*) '   INITIAL STEP FROM TAPE = ',IG2
	  WRITE(*,*) '   USE TAPE VALUE? (Y/N)'
	  READ(*,*) YN
	  IF ((YN .EQ. 'Y') .OR. (YN .EQ. 'y')) IG = IG2
	END IF
	WRITE(LOG,*) '>THE GRATING-RESOLUTION SETTING IS:'
	WRITE(LOG,'('' '',5X,I2,2X,I2)') IG,IR
	WRITE(*,*) 'THE GRATING-RESOLUTION SETTING IS:'
	WRITE(*,'('' '',5X,I2,2X,I2)') IG,IR
      ELSE
	IF (.NOT. L1) THEN
          WRITE(*,'(''** GRATING-RESOLUTION CHANGED TO '',I2,
     &		1X,I2,'' !?! **'')') IG2,IR
          WRITE(LOG,'(''>** GRATING-RESOLUTION CHANGED TO '',I2,
     &		1X,I2,'' !?! **'')') IG2,IR
	END IF
	IF (.NOT. L2) THEN
	  WRITE(*,*) '** ECOM DSS FLAGS CHANGED !?! **'
	  WRITE(LOG,*) '>** ECOM DSS FLAGS CHANGED !?! **'
        END IF
	WRITE(*,*) '** CHANGE INGORED. **'
	WRITE(LOG,*) '>** CHANGE IGNORED. **'
      END IF
      IGO = IG2
      IRO = IR
      IBL1 = ISPEC(1305)
      IBL2 = ISPEC(1306)
      IF (.NOT. L3) GO TO 110
*********1*********2*********3*********4*********5*********6*********7**
* GET WAVELENGTH ARRAY WAV, AND SENSITIVITY ARRAY S, CORRESPONDING
* TO THE PROPER GRATING STEP AND RESOLUTION MODE. THESE ARE STORED IN
* FILE caldta.dta.
*********1*********2*********3*********4*********5*********6*********7**
      CALL CONVRT( IG,IR,S,WAV )
*********1*********2*********3*********4*********5*********6*********7**
* OPEN THE grat.n FILE CONTAINING THE BIN NUMBERS FOR THE LINES AND
* BACKGROUND AREAS WE WANT.
* WLINE   = LINE WAVELENGTH
* WL1,WH1 = LOW AND HIGH CHANNEL LIMITS FOR LEFT (LOW WAVLNGTH) BACKGND
* WL,WH   = LOW AND HIGH CHANNEL LIMITS FOR LINE
* WL2,WH2 = LOW AND HIGH CHANNEL LIMITS FOR RIGHT (HI WAVLNGTH) BACKGND
*********1*********2*********3*********4*********5*********6*********7**
      CALL FOPEN(2,GRATFILE(IG))
      NUMLNS = 0
      DO NW = 1,MAXLNS
  105	READ (2,*,END = 107) WLINE(NW),WL1(NW),WH1(NW),
     &                        WL(NW),WH(NW),WL2(NW),WH2(NW)
        IF (WLINE(NW) .NE. 0.0) THEN
	  NUMLNS = NUMLNS + 1
	ELSE
	  GO TO 105
	END IF
      END DO
  107 WRITE(LOG,'(''>'',I2,'' LINES AVAILABLE WITH THIS GRATING STEP''
     &		  )') NUMLNS
      CLOSE (UNIT = 2)
      do nw=1,numlns
       write(*,*) nw,wline(nw)
      end do
      write(*,*)
      write(*,*) 'Enter the number of the wavelength to be processed'
      read(*,*) nline
*********1*********2*********3*********4*********5*********6*********7**
* CHANNEL LIMITS ARE STORED IN WL,WH, ETC. NOW FIND THE AVERAGE
* SENSITIVITY IN EACH RANGE, SB1, SB, SB2.
* DN1, AND DN2 ARE NORMALIZATION FACTORS WHICH TAKE INTO ACCOUNT THAT
* THE LINE AND BACKGROUND REGION MAY DIFFER IN SIZE.
*********1*********2*********3*********4*********5*********6*********7**
	WRITE(LOG,*) '>',WLINE(nline),WL1(nline),WH1(nline),WL(nline),
     &                           WH(nline),WL2(nline),WH2(nline)
        SB1 = 0.0
	SB = 0.0
        SB2 = 0.0
	DO I = WL1(nline),WH1(nline)
	  SB1 = SB1 + S(I) 
	END DO
	DO I = WL(nline),WH(nline)
	  SB = SB + S(I)
	END DO
	DO I = WL2(nline),WH2(nline)
	  SB2 = SB2 + S(I)
	END DO
        SB1 = SB1 / FLOAT( WH1(nline) - WL1(nline) + 1 )
        SB = SB / FLOAT( WH(nline) - WL(nline) + 1 )
	SB2 = SB2 / FLOAT( WH2(nline) - WL2(nline) + 1 )
	IF (SB1 .NE. 0.0)
     &    DN1 = FLOAT( WH(nline) - WL(nline) + 1 ) /
     &    (FLOAT( WH1(nline) - WL1(nline) + 1 ) * 2 * SB1 )
        IF (SB2 .NE. 0.0)
     &    DN2 = FLOAT( WH(nline) - WL(nline) + 1 ) /
     &    (FLOAT( WH2(nline) - WL2(nline) + 1 ) * 2 * SB2 )
      WRITE(*,*)
      WRITE(*,*) 'TIME RANGE FOUND. NOW READING DATA...'
      WRITE(LOG,*) 
*********1*********2*********3*********4*********5*********6*********7**
* SUM COUNTS AND BACKGROUND
*********1*********2*********3*********4*********5*********6*********7**
  110 CONTINUE
*********1*********2*********3*********4*********5*********6*********7**
* FOR EACH SPECTRUM, FIND SATELLITE POSITION, LOOK DIRECTION (ZENITH
* ANGLE), AND TIME OF OBSERVATION. IREP KEEPS TRACK OF WHICH SPECTRUM
* (OUT OF THE 10 IN THIS SET) IS THE CURRENT ONE.
*********1*********2*********3*********4*********5*********6*********7**
      DO IREP = 1,10
        ozen=zen
	ZEN =  90.0 - DTA(IREP * 10 + 24)
        AZI =  DTA(IREP * 10 + 23)
	TIME = DTA(IREP * 10 + 15)
        GLON = DTA(12)
        GLAT = DTA(11)
        BLAT = ASIND(DTA(10) / DTA(4) )
        SZEN = 90.0 - DTA(18)
        SAZI = DTA(17)
CCC        if(irep .eq. 1) write(*,*) time,szen,glon,glat
c        amag = dta(20)
c        zmag = 90.0 - dta(21)
c        zenm = acosd(cosd(zmag) * cosd(zen) + sind(zmag) *
c     &         sind(zen) *  cosd(azi - amag))
        PHASE = COSD(ZEN) * COSD(SZEN) + SIND(ZEN) *
     &          SIND(SZEN) * COSD(AZI - SAZI)
c        szen = zmag
*********1*********2*********3*********4*********5*********6*********7**
* NOW BLAT, GLAT, GLON, SZEN, TIME HAVE VALUES OF MAG LAT,
* GEO LAT, GEO LONG, SOLAR ZENITH ANGLE AND CORRESPONDING TIMES.
* ZEN AND AZI ARE THE LINE OF SIGHT ZENITH AND AZIMUTH ANGLES.
*********1*********2*********3*********4*********5*********6*********7**
	ITBN = ( INT(TIME - TML) / NSEC ) + 1
	IF (ITBN .GT. MAXBNS) ITBN = MAXBNS
	IF (ITBN .LT. 1) ITBN = 1
        if (zen .ge. ozen) then
         inc=1
        else
         inc=2
        end if
	izbn = ( INT(zen) / izbnw ) + 1
c	izbn = ( INT(zenm) / izbnw) + 1
	IF (izbn .GT. maxizbns) izbn = maxizbns
	IF (izbn .LT. 1) izbn = 1
        if (IZBN .le. 18) then
           f1200 = 0.083
        else
           f1200 = 0.083 + 0.002 * float(IZBN - 18)
        end if
*********1*********2*********3*********4*********5*********6*********7**
* NOW ADD SUMMED TOTALS TO THE PROPER ARRAY
* NW = LINE NUMBER, ITBN = TIME BIN, izbn = zenith ANGLE BIN
*********1*********2*********3*********4*********5*********6*********7**
	FNT1 = SUMF( WL1(nline),WH1(nline),IREP,ISPEC )
	FINT = SUMF( WL(nline),WH(nline),IREP,ISPEC )
	FNT2 = SUMF( WL2(nline),WH2(nline),IREP,ISPEC )
        if((wline(nline) - 1200.0) .eq. 0.0) then
          fnt1 = fnt1 * f1200
          fnt2 = fnt2 * f1200
        end if
	PCNT(ITBN,izbn,inc) = PCNT(ITBN,izbn,inc) + 1
        PBK1(ITBN,izbn,inc) = PBK1(ITBN,izbn,inc) + FNT1
        PLIN(ITBN,izbn,inc) = PLIN(ITBN,izbn,inc) + FINT
	PBK2(ITBN,izbn,inc) = PBK2(ITBN,izbn,inc) + FNT2
	ptzen(itbn,izbn,inc)=ptzen(itbn,izbn,inc)+zen
	sqzen(itbn,izbn,inc)=sqzen(itbn,izbn,inc)+zen*zen
        if ((azi .lt. -90.0) .and. (ptazim(itbn,izbn,inc)
     &   .gt. 0.0)) azi=azi+360.0
        if ((azi .gt. 90.0) .and. (ptazim(itbn,izbn,inc)
     &   .lt. 0.0)) azi=azi-360.0
	ptazim(itbn,izbn,inc)=ptazim(itbn,izbn,inc)+azi
	sqazim(itbn,izbn,inc)=sqazim(itbn,izbn,inc)+azi*azi
        if ((sazi .lt. -90.0) .and. (ptsaz(itbn,izbn,inc)
     &   .gt. 0.0)) sazi=sazi+360.0
        if ((sazi .gt. 90.0) .and. (ptsaz(itbn,izbn,inc)
     &   .lt. 0.0)) sazi=sazi-360.0
	ptsaz(itbn,izbn,inc)=ptsaz(itbn,izbn,inc)+sazi
	sqsaz(itbn,izbn,inc)=sqsaz(itbn,izbn,inc)+sazi*sazi
*********1*********2*********3*********4*********5*********6*********7**
* PTBLAT, PTGLAT, PTGLON, PTSZEN, WILL SUM UP BLAT, GLAT, ETC.
* FOR AVERAGING OVER EACH TIME BIN. PTTIME WILL BE SET TO THE
* MIDDLE OF THE TIME BIN LATER ON. (THIS ENSURES THAT PTTIME
* IS ALWAYS FILLED IN.)
*********1*********2*********3*********4*********5*********6*********7**
  120	continue
        if ((glon .lt. -90.0) .and. (ptglon(itbn,inc) .gt. 0.0)) then
         glon=glon+360.0
        end if
        if ((glon .gt. 90.0) .and. (ptglon(itbn,inc) .lt. 0.0)) then
         glon=glon-360.0
        end if
  	PTGLON(ITBN,inc) = PTGLON(ITBN,inc) + GLON
	PTGLAT(ITBN,inc) = PTGLAT(ITBN,inc) + GLAT
	PTBLAT(ITBN,inc) = PTBLAT(ITBN,inc) + BLAT
	PTSZEN(ITBN,inc) = PTSZEN(ITBN,inc) + SZEN
	sqglon(ITBN,inc) = sqglon(ITBN,inc) + GLON*glon
	sqglat(ITBN,inc) = sqglat(ITBN,inc) + GLAT*glat
	sqblat(ITBN,inc) = sqblat(ITBN,inc) + BLAT*blat
	sqszen(ITBN,inc) = sqszen(ITBN,inc) + SZEN*szen
      END DO
      OK = .TRUE.
*********1*********2*********3*********4*********5*********6*********7**
* GO BACK FOR ANOTHER DATA SET, REPORTING GOOD/BAD SET AND ROUGHLY
* HOW MUCH DATA HAS BEEN READ. INDEX COUNTS THE DATA IN 5% CHUNKS.
* ONE DATA SET = ONE ASPECT+SPECTUM SET = 10 SPECTRA =~ 1.92 SEC
*********1*********2*********3*********4*********5*********6*********7**
  130 IF (.NOT. OK) THEN
	IF (GOOD .NE. 0) THEN
	  WRITE(LOG,'(''READ '',I6,'' GOOD DATA SETS.'')') GOOD
	  GOOD = 0
        END IF
        BAD = BAD + 1
	TBAD = TBAD + 1
      END IF
      IF (OK) THEN
	IF (BAD .NE. 0) THEN
     	  WRITE(LOG,'(''SKIPPED '',I6,'' BAD DATA SETS.'')') BAD
          BAD = 0
	END IF
	GOOD = GOOD + 1
	TGOOD = TGOOD + 1
      END IF
      IF (ITBN .GE. NUMBNS * INDEX / 20) THEN
	WRITE(*,'(I3,''% OF DATA HAS BEEN PROCESSED'')')
     &		index * 5
	INDEX = INDEX + 1
      END IF
      GO TO 60
*********1*********2*********3*********4*********5*********6*********7**
* FINISHED READING DATA
* NOW AVERAGE PTGLON, PTGLAT, ETC. OVER EACH TIME BIN.
* PCNT(LINE#,TIMEBIN,ZENBIN) = NUMBER OF SPECTRA IN THAT BIN
* SINCE PCNT IS INDEPENDENT OF LINE #, WE USE LINE 1
* ALSO SET PTTIME TO THE MIDDLE OF EACH TIME BIN.
*********1*********2*********3*********4*********5*********6*********7**
  140 IF (GOOD .NE. 0) THEN
	WRITE(LOG,'(''READ '',I6,'' GOOD DATA SETS.'')') GOOD
      END IF
      IF (BAD .NE. 0) THEN
	WRITE(LOG,'(''SKIPPED '',I6,'' BAD DATA SETS.'')') BAD
      END IF
      DO ITBN = 1,NUMBNS
	PTTIME(ITBN) = TML + (ITBN - 1) * NSEC + (NSEC / 2.0)
        do inc=1,2
	  COUNT = 0
	  DO izbn = 1,nzbns
	    COUNT = COUNT + PCNT(ITBN,izbn,inc)
	  END DO
	  IF (COUNT .NE. 0) THEN
	    IF (COUNT .NE. 1) THEN
             fcnt=float(count)
             sdglon(itbn,inc)=sqrt(abs((fcnt*sqglon(itbn,inc)-
     &        ptglon(itbn,inc)**2)/fcnt/(fcnt-1.0)))
             sdglat(itbn,inc)=sqrt(abs((fcnt*sqglat(itbn,inc)-
     &        ptglat(itbn,inc)**2)/fcnt/(fcnt-1.0)))
             sdblat(itbn,inc)=sqrt(abs((fcnt*sqblat(itbn,inc)-
     &        ptblat(itbn,inc)**2)/fcnt/(fcnt-1.0)))
             sdszen(itbn,inc)=sqrt(abs((fcnt*sqszen(itbn,inc)-
     &        ptszen(itbn,inc)**2)/fcnt/(fcnt-1.0)))
            else
             sdglon(itbn,inc)=0.0
             sdglat(itbn,inc)=0.0
             sdblat(itbn,inc)=0.0
             sdszen(itbn,inc)=0.0
	    END IF
	    PTGLON(ITBN,inc) = PTGLON(ITBN,inc) / FLOAT(COUNT)
	    PTGLAT(ITBN,inc) = PTGLAT(ITBN,inc) / FLOAT(COUNT)
	    PTBLAT(ITBN,inc) = PTBLAT(ITBN,inc) / FLOAT(COUNT)
     	    PTSZEN(ITBN,inc) = PTSZEN(ITBN,inc) / FLOAT(COUNT)
	  END IF
        end do
      END DO
*********1*********2*********3*********4*********5*********6*********7**
* WRITE DATA TO FILE
*********1*********2*********3*********4*********5*********6*********7**
      OPEN (UNIT = 7, FILE =  ATAP // '.two',STATUS = 'UNKNOWN')
      OPEN (UNIT = 8, FILE =  ATAP // '.asp',STATUS = 'UNKNOWN')
C      OPEN (UNIT = 7, FILE =  ATAP // '.two',
C     & ACCESS = 'SEQUENTIAL', FORM = 'UNFORMATTED')
      WRITE(7,*) ATAP
      WRITE(7,*) RV
      WRITE(7,*) DY
      WRITE(7,*) YR
      WRITE(7,*) TML,TMH,NUMBNS,nzbns,NUMLNS,ig,ir,
     &  izbnw,NSEC
      WRITE(8,*) ATAP,RV,DY,YR
	WRITE(7,*) WLINE(nline),WL1(nline),WH1(nline),
     &		 WL(nline),WH(nline),WL2(nline),WH2(nline),
     &		 SB1,SB,SB2,DN1,DN2
      do inc=1,2
        DO ITBN = 1,NUMBNS
	  WRITE(7,*) PTTIME(ITBN),PTGLON(ITBN,inc),PTGLAT(ITBN,inc),
     &		   PTBLAT(ITBN,inc),PTSZEN(ITBN,inc)
	  WRITE(8,*) PTTIME(ITBN),PTGLON(ITBN,inc),PTGLAT(ITBN,inc),
     &		   PTBLAT(ITBN,inc),PTSZEN(ITBN,inc)
 	  WRITE(7,*) sdglon(ITBN,inc),sdglat(ITBN,inc),
     &		   sdblat(ITBN,inc),sdszen(ITBN,inc)
	  DO izbn = 1,nzbns
           fcnt=pcnt(itbn,izbn,inc)
           if (fcnt .ge. 2) then
             sdzen(itbn,izbn,inc)=
     &        sqrt(abs((fcnt*sqzen(itbn,izbn,inc)
     &        -ptzen(itbn,izbn,inc)**2)/fcnt/(fcnt-1.0)))
             sdazim(itbn,izbn,inc)=
     &        sqrt(abs((fcnt*sqazim(itbn,izbn,inc)
     &        -ptazim(itbn,izbn,inc)**2)/fcnt/(fcnt-1.0)))
             sdsaz(itbn,izbn,inc)=
     &        sqrt(abs((fcnt*sqsaz(itbn,izbn,inc)
     &        -ptsaz(itbn,izbn,inc)**2)/fcnt/(fcnt-1.0)))
           else
             sdazim(itbn,izbn,inc)=0.0
             sdsaz(itbn,izbn,inc)=0.0
           end if
           if (fcnt .ne. 0) then
             ptzen(itbn,izbn,inc)=ptzen(itbn,izbn,inc)/fcnt
             ptazim(itbn,izbn,inc)=ptazim(itbn,izbn,inc)/fcnt
             ptsaz(itbn,izbn,inc)=ptsaz(itbn,izbn,inc)/fcnt
           end if
	    WRITE(7,*) PCNT(ITBN,izbn,inc),PBK1(ITBN,izbn,inc),
     &	             PLIN(ITBN,izbn,inc),PBK2(ITBN,izbn,inc)
     	    WRITE(7,*) ptzen(itbn,izbn,inc),sdzen(itbn,izbn,inc),
     &		     ptazim(itbn,izbn,inc),sdazim(itbn,izbn,inc),
     &		     ptsaz(itbn,izbn,inc),sdsaz(itbn,izbn,inc)
	  END DO
        END DO
      end do

      DO ITBN = 1, NUMBNS
         DO IZBN =1, NZBNS
     	    WRITE(8,*) pttime(itbn), ptazim(itbn,izbn,1), ptazim(itbn,izbn,2)
         END DO
      END DO


      WRITE(7,*) S,WAV
      WRITE(*,*)
      WRITE(*,*) 'FINISHED WITH THIS TAPE.'
      WRITE(*,'(''DATA STORED IN FILE '',A20)')
     &   ATAP // '.two'
      WRITE(LOG,*)
      WRITE(LOG,*) 'FINISHED WITH THIS TAPE.'
      WRITE(LOG,'(''TOTAL NUMBER OF BAD DATA SETS: '',I5)') TBAD
      WRITE(LOG,'(''TOTAL NUMBER OF GOOD DATA SETS: '',I5)') TGOOD
      WRITE(LOG,'(''DATA STORED IN FILE '',A20)')
     &  '../DATA/' // ATAP // '.two'
      STOP
      END

*********1*********2*********3*********4*********5*********6*********7**
      FUNCTION SUMF( N1,N2,IREP,ISPEC )
      IMPLICIT INTEGER*2 (I-N)
      DIMENSION ISPEC(1308)
      COMMON /CORFAC/ FACTOR(128)
      SUM = 0.0
      DO I = N1,N2
	IND = (IREP - 1) * 128 + I + 4
        SUM = SUM + FLOAT( ISPEC(IND) ) * FACTOR(I)
      END DO
      SUMF = SUM
      RETURN
      END

*********1*********2*********3*********4*********5*********6*********7**
      SUBROUTINE FOPEN(UNITN,FNAME)
* SUBROUTINE TO OPEN EXISTING FILES
*********1*********2*********3*********4*********5*********6*********7**
      INTEGER UNITN,IERR
      CHARACTER*50 FNAME
      OPEN (UNIT = UNITN, FILE = FNAME,
     &      STATUS = 'OLD', IOSTAT = IERR)
      IF (IERR .NE. 0) THEN
        WRITE(*,*) 'IOSTAT ERROR WITH FILE ', FNAME,
     &			', IOSTAT = ', IERR
        STOP
      END IF
      RETURN
      END

*********1*********2*********3*********4*********5*********6*********7**
      CHARACTER*50 FUNCTION GRATFILE(IGRAT)
* FIND THE FILENAME CORRESPONDING TO THE GRATING STEP. YES, IT'S DUMB,
* BUT IT WORKS AND IT'S EASY TO UNDERSTAND.
*********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT INTEGER*2 (I-N)
      IF (IGRAT .EQ. 1) THEN
	GRATFILE = 'grat.1'
      ELSE IF (IGRAT .EQ. 2) THEN
	GRATFILE = 'grat.2'
      ELSE IF (IGRAT .EQ. 10) THEN
	GRATFILE = 'grat.10'
      ELSE IF (IGRAT .EQ. 13) THEN
	GRATFILE = 'grat.13'
      ELSE IF (IGRAT .EQ. 16) THEN
	GRATFILE = 'grat.16'
      ELSE IF (IGRAT .EQ. 19) THEN
	GRATFILE = 'grat.19'
      ELSE IF (IGRAT .EQ. 23) THEN
	GRATFILE = 'grat.23'
      ELSE IF (IGRAT .EQ. 33) THEN
	GRATFILE = 'grat.33'
      END IF
      RETURN
      END

*********1*********2*********3*********4*********5*********6*********7**
      SUBROUTINE ACCUM(irgn,SC,NSPEC,ACDATA)                                 
*                                                                       
*     BY SC ON AUG 22,1978. LAST CHANGE 15 AUG 1989 BY G. FRUTH
*                                                                       
*     THIS SUBROUTINE WILL ACCUMULATE COUNTS FROM SOME GRATING          
*     POSITION AND RESOLUTION MODE.                                     
*                                                                       
*     PARAMETERS :                                                      
*                                                                       
*     SC : THE PHOTON COUNT ARRAY OF SIZE 128 TO BE ACCUMULATED.               
*     IRGN : THE REGION WE ARE IN
*                                                                       
*     VARIABLES :                                                       
*        ACDATA  = THE ACCUMULATED DATA ARRAY                                
*	 NSPEC   = THE NUMBER OF SPECTRA ACCUMULATED
*********1*********2*********3*********4*********5*********6*********7** 
      IMPLICIT INTEGER*2 (I-N)                                          
      PARAMETER (MAXRGNS = 4)
      INTEGER*4 NSPEC(MAXRGNS)
      REAL SC(128),ACDATA(MAXRGNS,128)                                     
      COMMON /CORFAC/ FACTOR(128)
      NSPEC(IRGN) = NSPEC(IRGN) + 1
      DO ICHAN = 1,128                                                  
	ACDATA(IRGN,ICHAN) = ACDATA(IRGN,ICHAN) +
     &		SC(ICHAN) * FACTOR(ICHAN)
      END DO
      RETURN                                                            
      END                                                               

*********1*********2*********3*********4*********5*********6*********7**
      SUBROUTINE GETCOR
* SET UP THE CORRECTION FACTOR TABLE. THIS ROUTINE SETS UP ARRAY
* ELEMENTS 1-9 AND 30-128. ELEMENTS 10-29 ARE SET UP BY BLOCK DATA.
*********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT INTEGER*2 (I-N)
      COMMON /CORFAC/ FACTOR(128)
      DO I = 1,9
        FACTOR(I) = 0.0
      END DO
      DO I = 30,128
        FACTOR(I) = 1.0
      END DO
      RETURN
      END

*********1*********2*********3*********4*********5*********6*********7**
      BLOCK DATA
* INSERT VALUES INTO ELEMENTS 10-29 OF FACTOR
*********1*********2*********3*********4*********5*********6*********7**
      DIMENSION FACTR(20)
      COMMON /CORFAC/ FACTOR(128)
      EQUIVALENCE (FACTR(1),FACTOR(10))
      DATA FACTR /0.2532,0.3817,0.4444,0.5405,0.5882,0.6452,0.6897,
     &            0.7143,0.7576,0.7692,0.8000,0.8130,0.8333,0.8696,
     &            0.8696,0.8929,0.9091,0.9259,0.9259,0.9434/
      END                                                               

*********1*********2*********3*********4*********5*********6*********7** 

