c -------------------------------------------------------------------------
c  This file contains the subroutines for the emissivity and reflectivity.
c  These routines are called if an external file containing the spectral
c    emissivity or reflectivity is desired.  The data are then interpolated
c    to the necessary spectral resolution.
c
c  Use of external emissivity or reflectivity files is activated by
c    setting the first of the quadratic coefficients in TAPE5 to "-99"
c    (for emissivity and/or reflectivity)
c
c  This file contains the following subroutines:
c     EXTEMIS  = for computing emissivity from external file
c     EXTREFL  = for computing reflectivity from external file
c     RDERFILE = to read the external file
c     XINTARR  = to linearly interpolate input data to desired grid
c     LOCATE   = to locate the starting point for interpolation
c
c  Routines are called from the subroutines:
c     EMINIT, EMBND, FLINIT, RADNN
c
c  Requires the input files:
c     EMISIV.IN for emissivity
c     REFLEC.IN for reflectivity
c
c  Input file format:
c     line 1:  80 character header
c     line 2:  V1, V2, NPTS ( e13.6,e13.6,i4 )
c         V1 = starting wavenumber of data
c         V2 = ending wavenumber of data
c         NPTS = number of data points in file (max = MXPTS)
c     line 3 - n:  coefficients ( 6(e11.4,2x) )
c        coefficients are listed, 6 values per line, with equal
c         wavenumber spacing between the coefficients
c        HOWEVER, if v1 < 0 then both the coefficients and wavenumber
c         values are specified: v(1),c(1),v(2),c(2),v(3),c(3)
c
c  Call format:  VF, VL, INTVLS
c     VF = first wavenumber
c     VL = last wavenumber
c     INTVLS = number of spectral intervals = total points - 1
c       (used to match format of internal emis/refl calls)
c       maximum number given by the parameter MXINTV
c
c  Data are interpolated and put in the array EMIS or REFL which
c   are passed via the common block "EMRFVAL"
c
c  Written by:  Ned Snell, AER, Inc.   April 1995
c    Modified:  October 1995
c       added more double precision statements to fix roundoff problems with
c       the interpolation at high wavenumbers (tested 10000-13000 cm-1)
c
c -------------------------------------------------------------------------

c --- emissivity -----------------------------------------------------
      SUBROUTINE EXTEMIS(VF,VL,INTVLS)
      PARAMETER (MXPTS=2000, MXINTV=2400)
      DOUBLE PRECISION VF,VL,VCOEF,COEF
      COMMON /EMRFVAL/EMIS(MXINTV),REFL(MXINTV)
      COMMON /IFIL/ IRD,IPR,IPU,NOPR,NFHDRF,NPHDRF,NFHDRL,NPHDRL,
     &    NLNGTH,KFILE,KPANEL,LINFIL,NFILE,IAFIL,IEXFIL,
     &    NLTEFL,LNFIL4,LNGTH4 
      DIMENSION VCOEF(MXPTS),COEF(MXPTS)
      CHARACTER*80 HDR
      SAVE INITE,COEF,VCOEF,V1,V2,NPTS,DVCOE,IFRM
      DATA INITE /1/

c --- if inite=1 then the routine is in the first call
c         => open emissivity file and read data
      IF (INITE.EQ.1) THEN
          INITE=0
          OPEN(77,FILE='EMISIV.IN',TYPE='OLD',ERR=901)

c read txt header and write to TAPE6
          READ(77,1000) HDR
          WRITE(IPR,1001) HDR

          READ(77,1010) V1,V2,NPTS

c check form of input:
c   v1 < 0 => will specify a wavenumber for each coefficient
c      (coefficients in range abs(v1) - v2)
c   v1 > 0 => fixed spacing of coefficients between v1 and v2
          IF (V1.LE.0) THEN
              IFRM=2
              V1=ABS(V1)
          ELSE
              IFRM=1
c determine spectral spacing of coefficients
              DVCOE=(V2-V1)/FLOAT(NPTS-1)
c fill wavenumber array
              DO I=1,NPTS
                  VCOEF(I)=V1+FLOAT(I-1)*DVCOE
              ENDDO
          ENDIF

c read data from file
          CALL RDERFILE(77,IFRM,NPTS,VCOEF(1),COEF(1))

          CLOSE(77)

      ENDIF
c --- done with initialization

c check to be sure that stored data is in range of file
      IF ((VF.LT.V1).OR.(VL.GT.V2)) THEN
          WRITE(IPR,2000) V1,VF,V2,VL
          WRITE(*,2000) V1,VF,V2,VL
          STOP 'ERROR ON EMISSIVITY INPUT RANGE'
      ENDIF
      IF (INTVLS.GT.MXINTV) THEN
          WRITE(*,*) 'INTVLS, MAX = ',INTVLS,MXINTV
          STOP 'ERROR'
      ENDIF

c determine output dv
      IF (INTVLS.GT.0) THEN
          DVOUT=(VL-VF)/FLOAT(INTVLS)
      ELSE
          IF (INTVLS.LT.0) THEN
              WRITE(*,*) 'NLIM2 < NLIM1'
              STOP 'ERROR: INTVLS < 0'
          ENDIF
      ENDIF
ccccct
        write(ipr,*)'BEFORE XINTARR',vcoef(1),coef(1),npts,vf,  dvout, intvls
      
c interpolate stored coefficients to correct resolution
      CALL XINTARR(VCOEF(1),COEF(1),NPTS,VF,DVOUT,EMIS(1),INTVLS)
ccccct
        write(ipr,*)'AFTER XINTARR',vcoef(1),coef(1),npts,vf,  dvout, intvls

      RETURN

c error messages
  901 WRITE(IPR,1901)
      STOP 'ERROR OPENING EMISSIVITY FILE'

 1000 FORMAT(A80)
 1001 FORMAT(2X,'READING EMISSIVITY FROM A FILE:',/4X,A80)
 1010 FORMAT(E13.6,E13.6,I4)
 1901 FORMAT(//2X,'ERROR OPENING EMISSIVITY FILE:  EMISIV.IN')
 2000 FORMAT(//2X,'ERROR ON EMISSIVITY INPUT RANGE:',
     &    /2X,'V1(FILE): ',1PE13.6,2X,'V1(DESIRED): ',1PE13.6,
     &    /2X,'V2(FILE): ',1PE13.6,2X,'V2(DESIRED): ',1PE13.6)

      END

c --- reflectivity ---------------------------------------------------
      SUBROUTINE EXTREFL(VF,VL,INTVLS)
      PARAMETER (MXPTS=2000, MXINTV=2400)
      DOUBLE PRECISION VF,VL,VCOEF,COEF
      COMMON /EMRFVAL/EMIS(MXINTV),REFL(MXINTV)
      COMMON /IFIL/ IRD,IPR,IPU,NOPR,NFHDRF,NPHDRF,NFHDRL,NPHDRL,
     &    NLNGTH,KFILE,KPANEL,LINFIL,NFILE,IAFIL,IEXFIL,
     &    NLTEFL,LNFIL4,LNGTH4 
      DIMENSION VCOEF(MXPTS),COEF(MXPTS)
      CHARACTER*80 HDR
      SAVE INITR,COEF,VCOEF,V1,V2,NPTS,DVCOE,IFRM
      DATA INITR /1/

c --- if init=1 then the routine is in the first call
c         => open reflectivity file and read data
      IF (INITR.EQ.1) THEN
          INITR=0
          OPEN(77,FILE='REFLEC.IN',TYPE='OLD',ERR=901)

c read txt header and write to TAPE6
          READ(77,1000) HDR
          WRITE(IPR,1001) HDR

          READ(77,1010) V1,V2,NPTS

c check form of input:
c   v1 < 0 => will specify a wavenumber for each coefficient
c       (coefficients in the range abs(v1) - v2)
c   v1 > 0 => fixed spacing of coefficients between v1 and v2
          IF (V1.LE.0) THEN
              IFRM=2
              V1=ABS(V1)
          ELSE
              IFRM=1
c determine spectral spacing of coefficients
              DVCOE=(V2-V1)/FLOAT(NPTS-1)
c fill wavenumber array
              DO I=1,NPTS
                  VCOEF(I)=V1+FLOAT(I-1)*DVCOE
              ENDDO
          ENDIF

c read data from file
          CALL RDERFILE(77,IFRM,NPTS,VCOEF(1),COEF(1))

          CLOSE(77)

      ENDIF
c --- done with initialization

c check to be sure that stored data is in range of file
      IF ((VF.LT.V1).OR.(VL.GT.V2)) THEN
          WRITE(IPR,2000) V1,VF,V2,VL
          WRITE(*,2000) V1,VF,V2,VL
          STOP 'ERROR ON REFLECTIVITY INPUT RANGE'
      ENDIF
      IF (INTVLS.GT.MXINTV) THEN
          WRITE(*,*) 'INTVLS, MAX = ',INTVLS,MXINTV
          STOP 'ERROR'
      ENDIF

c determine output dv
      IF (INTVLS.GT.0) THEN
          DVOUT=(VL-VF)/FLOAT(INTVLS)
      ELSE
          IF (INTVLS.LT.0) THEN
              WRITE(*,*) 'NLIM2 < NLIM1'
              STOP 'ERROR: INTVLS < 0'
          ENDIF
      ENDIF

c interpolate stored coefficients to correct resolution
      CALL XINTARR(VCOEF(1),COEF(1),NPTS,VF,DVOUT,REFL(1),INTVLS)

      RETURN

c error messages
  901 WRITE(IPR,1901)
      STOP 'ERROR OPENING REFLECTIVITY FILE'

 1000 FORMAT(A80)
 1001 FORMAT(2X,'READING REFLECTIVITY FROM A FILE:',/4X,A80)
 1010 FORMAT(E13.6,E13.6,I4)
 1901 FORMAT(//2X,'ERROR OPENING REFLECTIVITY FILE:  EMISIV.IN')
 2000 FORMAT(//2X,'ERROR ON REFLECTIVITY INPUT RANGE:',
     &    /2X,'V1(FILE): ',1PE13.6,2X,'V1(DESIRED): ',1PE13.6,
     &    /2X,'V2(FILE): ',1PE13.6,2X,'V2(DESIRED): ',1PE13.6)

      END

c --- utility subroutines (common to emissivity and reflectivity) ---
c -------------------------------------------------------------------------
      SUBROUTINE RDERFILE(IUER,IFRM,NPTS,VCOEF,COEF)
c
c subroutine to read data from the emissivity or reflectivity file
c
c  written by:  Ned Snell, AER, Inc.   April 1995
c
      DOUBLE PRECISION VCOEF,COEF
      DIMENSION VCOEF(*),COEF(*)

c number of points per line dictated by form of input
c   => number of lines in file depends on form
c 6 values per line if no wavenumbers
c 3 values per line if wavenumber-coefficient pairs are given
      NRD=(FLOAT(NPTS*IFRM)/6.0)+0.9999

      J=1

      IF (IFRM.EQ.1) THEN
          DO I=1,NRD
              READ(IUER,1000,END=999) (COEF(K),K=J,J+5)
              J=J+6
          ENDDO
      ELSE
          DO I=1,NRD
              READ(IUER,1000,END=999) VCOEF(J),COEF(J),
     &            VCOEF(J+1),COEF(J+1),VCOEF(J+2),COEF(J+2)
              J=J+3
          ENDDO
      ENDIF

      RETURN

  999 STOP 'ERROR ON EMISSIVITY OR REFLECTIVITY INPUT FILE'

 1000 FORMAT(6(E11.4,2X))
      END
c -------------------------------------------------------------------------
      SUBROUTINE XINTARR(VC,COEF,NPTS,VF0,DVOUT,COUT,INTVLS)
c
c  subroutine to linearly interpolate between pairs of coefficients
c
c  written by:  Ned Snell, AER, Inc.   April 1995
c

      DOUBLE PRECISION VF0,VC,VF,COEF
      DIMENSION VC(*),COEF(*),COUT(*)

      VF=VF0
c start interpolation at vf - find pair of coefficients that 
c  bracket the starting point
      CALL LOCATE(VC(1),NPTS,VF,J)

c check for error (should have been stopped before this point,
c  but just in case....
      IF ((J.EQ.0).OR.(J.EQ.NPTS)) THEN
          PRINT *,VC(1),VC(NPTS),NPTS,VF,J
          STOP 'ERROR IN XINTARR - FILE OUT OF RANGE'
      ENDIF

c if intvls = 0, then want coefficient for a single point
      IF (INTVLS.EQ.0) THEN
          VC1=VC(J)
          VC2=VC(J+1)
          DCDV=(COEF(J+1)-COEF(J))/(VC2-VC1)
          COUT(1)=COEF(J)+DCDV*(VF-VC1)
          RETURN
      ENDIF

c starting interpolation range is between coefficient pair
c  at wavenumbers vc(j) and vc(j+1)
c initialize cout array counter
      ILOW=1

      INTVLS=INTVLS+1

  100 VC1=VC(J)
      VC2=VC(J+1)
      IPTS=((VC2-VF)/DVOUT)+1
c ipts is the number of points between vf and vc2 (where vc2 is the
c  last point before a new pair of coefficients is required)
c intvls is the number of points needed before returning to the main
c  part of the program
c icnt is set to either ipts or intvls depending on the location of vf
c  relative to vc2 (if icnt < intvls, then choose next set of coefficients
c  and finish loading the array)
      ICNT=MIN(IPTS,INTVLS)
      ICNT=ICNT-1

      DCDV=(COEF(J+1)-COEF(J))/(VC2-VC1)
      
c interpolation:
c    coef(i)=coef(j)+dcdv*deltv
c    deltv=v(i)-vc1=vf+(i-1)*dvout-vc1
c thus,
c    dcdv*deltv=dcdv*(v(i)-vc1)=dcdv*(vf-vc1)+dcdv*(dvout*(i-1))
c                               ---- a ------ ---- b ----
      A=DCDV*(VF-VC1)
      B=DCDV*DVOUT

c interpolation loop - fill cout array
      DO I=ILOW,ILOW+ICNT
          COUT(I)=COEF(J)+A+(B*FLOAT(I-ILOW))
      ENDDO

c see if more points are needed (move to next pair of coefficients)
      ILOW=ILOW+ICNT
      IF (ILOW.LT.INTVLS) THEN
          J=J+1
          VF=VF0+FLOAT(ILOW-1)*DVOUT
          INTVLS=INTVLS-ILOW-1
          ILOW=ILOW+1
c if intvls = 0, then want coefficient for a single point
          IF (INTVLS.EQ.0) THEN
              VC1=VC(J)
              VC2=VC(J+1)
              DCDV=(COEF(J+1)-COEF(J))/(VC2-VC1)
              COUT(ILOW)=COEF(J)+DCDV*(VF-VC1)
              RETURN
          ENDIF
          GOTO 100
      ENDIF

      RETURN
      END
c -------------------------------------------------------------------------
      SUBROUTINE LOCATE(XX,N,X,J)
c
c Given the array xx of length n, and given the value x, this routine
c  will return a value of j such that x is between xx(j) and xx(j+1)
c  if j=0 or j=n then the value x is out of range.
c The lengthy if-statement is used so that the array xx may be either
c  monotonically increasing or decreasing, though for fase
c  emissivity or reflectivity it must always be increasing.
c
c  written by:  Ned Snell, AER, Inc.   April 1995
c  based on the routine in "Numerical Recipes" by Press et.al.
c

      DOUBLE PRECISION XX,X
      DIMENSION XX(*)

c check for exact match
      DO I=1,N
          IF (X.EQ.XX(I)) THEN
              J=I
              RETURN
          ENDIF
      ENDDO

c initialize limits
      JL=0
      JU=N+1

c search
  10  IF ((JU-JL).GT.1) THEN
          JM=(JU+JL)/2

c    "if ((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm))) then"
c  .eqv. may not be standard fortran
c true if both logical operands are either both true or
c  both false => replace with the following
          IF ( ((XX(N).GT.XX(1)).AND.(X.GT.XX(JM))).OR.
     &         ((XX(N).LE.XX(1)).AND.(X.LE.XX(JM))) ) THEN
              JL=JM
          ELSE
              JU=JM
          ENDIF
          GOTO 10
      ENDIF

c set output
      J=JL

      RETURN
      END
c -------------------------------------------------------------------------