        PROGRAM  AEFILTER
C
C       This program reads an AERI radiance file (SSEC format)
C       and applies a filter function to match a given radiometer channel.
C       Created  17 April 1994  ROK
C       Modified 19 April 1994 ROK
C       Modified 19 March 1996 ROK
C       Modified 20 March 1996 ROK to include proper Eq B.T. for filter
C       Modified 22 March 1996 ROK to include decimal day of month
C       Modified 24 March 1996 ROK to fix bug in FNORM.
C
        PARAMETER (MAXFLT=100, MAXINT=5000, MAXDAT=5000,MAXNBT=2000)
        REAL    FLTWNS(MAXFLT),FLTPTS(MAXFLT),FLTINT(MAXINT)
        REAL    BTABLE(MAXNBT),HDR(100), RAD(MAXDAT)
        INTEGER NUMFLT,NUMINT
        INTEGER*2 DosSleep, ISTAT
        INTEGER*4 ARG
        CHARACTER*40 FLTFIL, FLTLAB, FLTOUT, INFILE, OUTFIL
        CHARACTER*78 COMMENT
        LOGICAL FILE_EXIST
        DATA  INP /1/, IOUT /2/, IFLT /3/
        DATA INIARG/1024/, MAXARG/16384/
        WRITE(*,*)
        WRITE(*,*) '   PROGRAM  AEFILTER'
        WRITE(*,*)
        WRITE(*,*) '   Version 1.4 24 March 1996'
        WRITE(*,*)
C ...........................................................
  2     WRITE(*,*) ' ENTER filename containing Filter:'
        READ(*,'(A)') FLTFIL
        OPEN(IFLT,FILE=FLTFIL,STATUS='UNKNOWN',ERR=2)
        WRITE(*,*) ' ... successfully opened file ',FLTFIL
        READ(IFLT,'(A)') FLTLAB
        WRITE(*,*) ' Filter Label: ', FLTLAB
        READ(IFLT,*) NUMFLT
        WRITE(*,*) ' Number of Filter Points: ',NUMFLT
        WRITE(*,*) ' <<< Press Enter >>> '
        READ(*,*)
        DO I = 1, NUMFLT
          READ(IFLT,*) FLTWNS(I), FLTPTS(I)
        ENDDO
        CLOSE(IFLT)
        WRITE(*,*) ' Number of points =',NUMFLT
        WRITE(*,*) ' (Wavenumber, Filter Value, DWN, PSUM)'
        PSUM = 0
        DWN = 0.5*(FLTWNS(2)-FLTWNS(1))
        PSUM = PSUM + FLTPTS(1)*DWN
        WRITE(*,90) FLTWNS(1),FLTPTS(1), DWN, PSUM
        DO I = 2, NUMFLT-1
          DWN = (FLTWNS(I) - FLTWNS(I-1))/2 +
     *          (FLTWNS(I+1) - FLTWNS(I))/2
          PSUM = PSUM + FLTPTS(I)*DWN
          WRITE(*,90) FLTWNS(I),FLTPTS(I), DWN, PSUM
        ENDDO
        DWN = 0.5*(FLTWNS(NUMFLT)-FLTWNS(NUMFLT-1))
        PSUM = PSUM + FLTPTS(NUMFLT)*DWN
        WRITE(*,90) FLTWNS(NUMFLT),FLTPTS(NUMFLT), DWN, PSUM
  90    FORMAT(1X,4F12.5)
        WRITE(*,*) ' Approx. Filter Normalization =',PSUM
        WRITE(*,*)
C ...........................................................
        WRITE(*,*) ' ENTER File Name Containing Radiance Data:'
        READ(*,'(A)') INFILE
C ... Initially open file to read 100 word header.
c ...    Check that file exists before opening first time
10       inquire(file=infile,exist=file_exist,err=15)
         if( .not.file_exist )goto 15
         ARG = INIARG
         GO TO 19
15       CONTINUE
         ARG=ARG*2
         IF(ARG.GT.MAXARG)ARG=MAXARG
         WRITE(*,*)'Waiting for creation of file ', infile
         ISTAT=DosSleep(ARG)
         GO TO 10
19       CONTINUE
C ...   OPEN input file
        OPEN(INP,FILE=INFILE,ACCESS='DIRECT',RECL=400)
        WRITE(*,*) ' ... succesfully opened file ',INFILE
C ...       READ first record to get starting Wavenumber and DV
        IREC = 1
        READ(INP,REC=IREC,ERR=999) (HDR(I),I=1,100)
        CLOSE(INP)
        DV   = HDR(33)
        WNA  = HDR(34)
        WNB  = HDR(35)
        NWPR = HDR(36)
        NHDRW= HDR(37)
        NPTS = NWPR - NHDRW
C ...........................................................
        WRITE(*,*) ' Set up to interpolate.'
C ...       Set up interpolated filter function
        V0 = WNA
        V1 = WNB
        CALL INTFLT(NUMFLT,FLTWNS,FLTPTS,V0,V1,DV,NUMINT,FLTINT)
        PSUM = 0.0
        DO I = 1, NPTS
          PSUM = PSUM + FLTINT(I)*DV
        ENDDO
        FNORM = PSUM
C ...       Write interpolated filter to a data file for later plotting.
        FLTOUT = 'FLTINT.DAT'
        OPEN(IFLT,FILE=FLTOUT,STATUS='UNKNOWN')
        WRITE(*,*) ' ... writing interpolated filter to ',FLTOUT
        WRITE(IFLT,*) 'INTERPOLATED FILTER FROM ',FLTLAB
        DO I = 1, NPTS
          V = V0 + DV*FLOAT(I-1)
          WRITE(IFLT,*) V, FLTINT(I)
        ENDDO
        CLOSE(IFLT)
C ...........................................................
C ...   Compute Table of filtered radiance versus Eq. B.T.
        WRITE(*,*) '... Compute table of Equiv. B.T.'
        NBTPTS =  200
        TLO    = 125.00
        DT     =   1.0
C ...               Loop over temperatures from LO to HI (Kelvin)
C ...               IBTPTR is the temperature index pointer
        DO IBTPTR = 1, NBTPTS
          T = TLO + DT*FLOAT(IBTPTR-1)
          RSUM = 0.0
          DO I = 1, NPTS
            V = V0 + DV*FLOAT(I-1)
            IF (FLTINT(I) .NE. 0.0)
     *      RSUM = RSUM + HPLANW(V,T)*FLTINT(I)*DV
          ENDDO
          PRADAVG  = RSUM/FNORM
          BTABLE(IBTPTR) = PRADAVG
        ENDDO
C ...........................................................
        WRITE(*,*) ' ENTER File Name for OUTPUT File:'
        READ(*,'(A)') OUTFIL
C ...   OPEN output file
        NBPR = 4*NWPR
        OPEN(INP,FILE=INFILE,ACCESS='DIRECT',RECL=NBPR)
        OPEN(IOUT,FILE=OUTFIL,STATUS='UNKNOWN')
C ...........................................................
        WRITE(*,*) 'Please ENTER a comment line:'
        READ(*,'(A)') COMMENT
        WRITE(IOUT,93) COMMENT
        WRITE(IOUT,93) 
     * '  DATE HRMISE  MoDyFract  UTC HOUR',
     * '     WNAVG  RADIANCE   BT(K)  BT(C)'
 93       FORMAT(80A)
        WRITE(*,*) ' ... applying filter to input radiance'
        WRITE(*,*) '  REC   TABLE(C)   PLANCK(C)   TABLE-PLANCK(C) '
        IREC = 1
        ARG = INIARG
 100    READ(INP,REC=IREC,ERR=999) 
     *           (HDR(I),I=1,100),(RAD(I),I=1,NPTS)
          DATE = HDR(22)
          TIME = HDR(23)
          IHR   = TIME/10000
          IMI   = (TIME-IHR*10000)/100
          ISE   = TIME - IHR*10000 - IMI*100
          ITIME = TIME
          IDATE = DATE
          IYR   = DATE/10000
          IMO   = (DATE-IYR*10000)/100
          IDY   = DATE - IYR*10000 - IMO*100
          DECHR = IHR + IMI/60. + ISE/3600.
          DECDAT= FLOAT(IMO*100+IDY)+DECHR/24.
C ...........................................................
C  RSUM = Partial sum of weighted radiance average
C  VSUM = Partial sum of weighted central wavenumber
C ...........................................................
C          WRITE(*,*) ' ... applying FILTER to Record #',IREC
          RSUM = 0.0
          VSUM = 0.0
          DO I = 1, NPTS
            V = V0 + DV*FLOAT(I-1)
            VSUM = VSUM + V*FLTINT(I)*DV
            RSUM = RSUM + RAD(I)*FLTINT(I)*DV
          ENDDO
          WNAVG   = VSUM/FNORM
          RADAVG  = RSUM/FNORM
C ...               Use look up table to get B.T. in Kelvin
C          WRITE(*,*) ' ... look up BT: RAD =',RADAVG
          BTAVGK  = GETBT(BTABLE,TLO,NBTPTS,DT,RADAVG)
          BTAVGC  = BTAVGK - 273.15
C ...               Use Inverse Planck function at mean wavenumber
          HBTAVK = HBRITW(WNAVG,RADAVG)
          HBTAVC = HBTAVK - 273.15
          WRITE(*,'(I5,3F11.4)') IREC, BTAVGC, HBTAVC, BTAVGC-HBTAVC

          WRITE(IOUT,95) IDATE,ITIME,DECDAT,
     *                   DECHR,WNAVG,RADAVG,BTAVGK,BTAVGC
 95       FORMAT(I6,1X,I6.6,F11.5,F10.5,F10.3,F10.4,F8.2,F7.2)
          IREC = IREC + 1
          ARG = INIARG
        GOTO 100
C ...      SLEEP
 999    CLOSE(IOUT)
        WRITE(*,*) '  ###  normal termination ###'
C ...........................................................
        END

        SUBROUTINE INTFLT(NUMFLT,FLTWNS,FLTPTS,V0,V1,DV,NUMINT,FLTINT)
C
C       Linearly Interpolate Points in FLTPLTS into array FLTINT
C
        INTEGER NUMFLT, NUMINT
        REAL    V0,V1,DV
        REAL    FLTWNS(*),FLTPTS(*),FLTINT(*)

        NUMINT = INT((V1-V0)/DV) + 1
        DO J = 1, NUMINT
          FLTINT(J) = 0.0
        ENDDO
        DO J = 1, NUMINT
          VINT = V0 + DV*FLOAT(J-1)
          DO I =  1, NUMFLT-1
            IF (VINT .GE. FLTWNS(I) .AND. VINT .LE. FLTWNS(I+1)) THEN
              FLTINT(J) = ((FLTPTS(I+1)-FLTPTS(I))/
     *                     (FLTWNS(I+1)-FLTWNS(I)))*(VINT-FLTWNS(I))
     *                    + FLTPTS(I)
            ENDIF
          ENDDO
        ENDDO
        RETURN
        END
C ... $TITLE:'HBRITW.FOR, 13 JAN 94, HBH'
C ... $DEBUG
      FUNCTION HBRITW(WN,RAD)                                                 
C *** HAL WOOLF'S SOFTWARE. . .
C ... Modified 3/24/96 to avoid divide by zero. ROK
C $ 'HIS' ... RADIANCE TO BRIGHTNESS TEMPERATURE                              
      DATA H/6.6237E-27/,C0/2.99791E+10/,BK/1.38024E-16/                      
      DATA INIT/1/                                                            
      TNB(X,Y,Z)=Y/ALOG(X/Z+1.)                                               
      IF(INIT.EQ.0) GO TO 100                                                 
      C1=2.*H*C0**2                                                           
      C2=H*C0/BK                                                              
      INIT=0                                                                  
  100 F1=C1*WN**3                                                             
      F2=C2*WN                                                                
      IF(F1 .GT. 1.E-5 .AND. RAD .GT. 1.E-5) THEN
        TBB=TNB(F1,F2,RAD)
      ELSE
        TBB=0.
      ENDIF
      HBRITW=TBB                                                              
      RETURN                                                                  
      END                                                                     
C 
c  $TITLE:'HPLANW.FOR,3-8-88'
      FUNCTION HPLANW(WN,TEM)
C *** HAL WOOLF'S SOFTWARE. . .
C $ 'HIS' ... TEMPERATURE TO PLANCK RADIANCE
      DATA H/6.6237E-27/,C0/2.99791E+10/,BK/1.38024E-16/
      DATA INIT/1/
      BNT(X,Y,Z)=X/(EXP(Y/Z)-1.)
      IF(INIT.EQ.0) GO TO 100
      C1=2.*H*C0**2
      C2=H*C0/BK
      INIT=0
  100 F1=C1*WN**3
      F2=C2*WN
      RAD=BNT(F1,F2,TEM)
      HPLANW=RAD
      RETURN
      END
C
      FUNCTION  GETBT(BTABLE,TLO,NBTPTS,DT,RAD)
C
C     This function interpolates the table value for the
C     input radiance and returns the corresponding temperature.
C
      PARAMETER (SMALL=2.E-5)
      REAL*4 BTABLE(*)
      
      IF (RAD .LT. BTABLE(1))      RESULT=BTABLE(1)
      IF (RAD .GT. BTABLE(NBTPTS)) RESULT= BTABLE(NBTPTS)
      DO IBTPTR = 1, NBTPTS-1
       IF (((RAD+SMALL) .GE. BTABLE(IBTPTR)) .AND.
     *     ((RAD+SMALL) .LT. BTABLE(IBTPTR+1))) THEN
         ISAVE=IBTPTR
         GOTO 10
       ENDIF
      ENDDO
10    CONTINUE
      T0 = TLO + DT*FLOAT(ISAVE-1)
      T1 = TLO + DT*FLOAT(ISAVE-1+1)
      RESULT=((T1-T0)/(BTABLE(ISAVE+1)-BTABLE(ISAVE)))*
     *       (RAD-BTABLE(ISAVE)) + T0
      GETBT = RESULT
      RETURN
      END
 
