      SUBROUTINE gfluxinterp(oz, sz, pr, rs, lat, ilam0, ilam1, 
     &                      flux, success)
c------------------------------------------------------------------------
c Calling parameters
c
c  Input:
      REAL*8 oz,      ! ozone (DU)
     &       sz,      ! solar zenith angle (deg)
     &       pr,      ! terrain pressure (atm)
     &       rs,      ! surface reflectivity (0 .le. rs .le. 1)
     &       lat      ! latitude (deg) [not used in current implement'n]
      INTEGER*4 ilam0,  ! first wavelength index to use in RT tables
     &          ilam1   ! last wavelength index to use.

c  Output:
      REAL*8 flux(0:ilam1-ilam0)  ! interpolated instantaneous spectral flux
      LOGICAL success             ! problem flag. .TRUE. if everything ok.

c========================================================================
c GFLUXINTERP 
c
c  Performs interpolation on continuous scan tables with respect to
c  oz (ozone), sz (solar zenith angle, degrees), pr (terrain pressure),
c  and rs (surface reflectivity) over the range of wavelength indices 
c  (ilam0..ilam1).
c
c  Returns the global spectral flux (flux).  This array
c  must be dimensioned long enough for results at
c  all the wavelengths in the range (ilam0..ilam1).  Also returns a logical
c  variable, success, to indicate successful completion of the calculation.
c
c  Interpolation is performed using Lagrange interpolation.  Since several
c  quantities are being interpolated with the same coefficients, these
c  coefficients are pre-calculated (by routine int_range).  Interpolation
c  of the fluxes with respect to solar zenith angle is
c  actually done by interpolating ln(flux) with respect to 
c  ln(cos(sza)).
c
c  Current implementation uses only mid-latitude ozone/temperature profile
c  fluxes.
c  
c------------------------------------------------------------------------
c Common blocks
c
      INCLUDE 'params.inc'
      INCLUDE 'rttbls_cs.cmn'

c------------------------------------------------------------------------
c Local variables
c
      REAL*8 radeg      /57.295779513082323/    !deg-rad conversion
      INTEGER*4 inrange, sz0, szdeg, oz0, ozdeg
      REAL*8 poz, psz, ppr, prod
      REAL*8 lcs, ozn
      REAL*8 szden(0:3), ozden(0:3), prden(0:1)
      REAL*8 flux_p

      INTEGER*2 printcall /0/

c-----------------------------------------------------------------------
c Print calling parameters.
c
      IF (printcall .EQ. 1) THEN
        WRITE(6,94) oz, sz, pr, rs, lat, ilam0, ilam1
      END IF

c-----------------------------------------------------------------------
c Idiot-check ilam0, ilam1
c
      IF(ilam0 .LT. 0 .OR. ilam0 .GT. n_waveln-1) GOTO 901
      IF(ilam1 .LT. 0 .OR. ilam1 .GT. n_waveln-1) GOTO 902
      IF(ilam0 .GT. ilam1) GOTO 903

c-----------------------------------------------------------------------
c Set up for SZA interpolation
c
      lcs= log(cos(sz/radeg))
d     WRITE(6,*)'sz:',sz,'  lcs:',lcs
      CALL int_rang(lcs, lncossza, 10, sz0, szdeg, 
     &                                            szden, inrange)

d     WRITE(6,990)lcs, n_sza, (lncossza(i),i=0,n_sza-1)
d     WRITE(6,991) sz0, szdeg, inrange, (szden(i),i=0,szdeg)

c-----------------------------------------------------------------------
c Set up for OZONE interpolation
c
      CALL int_rang(oz, ozone(6), 10, oz0, ozdeg, 
     &                                            ozden, inrange)

d     WRITE(6,992) oz, (ozone(j),j=6,15)
d     WRITE(6,993)oz0, ozdeg, inrange, (ozden(j), j=0,ozdeg)

c-----------------------------------------------------------------------
c Set up for PRESSURE interpolation
c  (comment out one of the following two...)
c .. for interpol'n wrt linear pressure

c     prden(0)= press(0) - press(1)
c     prden(1)= -prden(0)

c .. for interpol'n wrt logarithmic pressure

      prden(0)= log(press(0)) - log(press(1))
      prden(1)= -prden(0)

c-----------------------------------------------------------------------
c Loop through the desired wavelengths, do the interpolations in the 
c local variables, *_x, and stuff the results in the output arrays.
c
      idest=0
      DO ilam= ilam0, ilam1
        flux_x=  0.d0

        DO ipr= 0, 1
          flux_p= 0.d0

c .. use the following line for interpol'n wrt linear pressure

c         ppr= (pr - press(1-ipr))/prden(ipr)

c .. or for logarithmic pressure, use

          ppr= (log(pr) - log(press(1-ipr)))/prden(ipr)



          DO ioz= 0, ozdeg
            poz= 1./ozden(ioz)
            DO koz= 0, ozdeg
              IF (koz .NE. ioz) poz=poz*(oz-ozone(6+oz0+koz))
            END DO

            DO isz= 0, szdeg
              psz= 1./szden(isz)
              DO ksz= 0, szdeg
                IF (ksz .NE. isz) psz=psz*(lcs-lncossza(sz0+ksz))
              END DO

              prod= poz*psz

c     WRITE(6,97)sz0, isz, oz0, ioz, ilam, ipr, rs, poz, psz, ppr,
c    &       fs(sz0+isz,6+oz0+ioz,ilam,ipr),
c    &       gg(sz0+isz,6+oz0+ioz,ilam,ipr),
c    &       sb(6+oz0+ioz, ilam, ipr),
c    &       (fs(sz0+isz,6+oz0+ioz,ilam,ipr) + 
c    &        gg(sz0+isz,6+oz0+ioz,ilam,ipr)  )/
c    &         (1.0 - rs*sb(6+oz0+ioz, ilam, ipr)),
c    &           log (( fs(sz0+isz,6+oz0+ioz,ilam,ipr) +
c    &                  gg(sz0+isz,6+oz0+ioz,ilam,ipr)  )/
c    &                     (1.0 - rs*sb(6+oz0+ioz, ilam, ipr)))


              flux_p= flux_p + prod* 
     &           log (( fs(sz0+isz,6+oz0+ioz,ilam,ipr) + 
     &                  gg(sz0+isz,6+oz0+ioz,ilam,ipr)  )/
     &                     (1.0 - rs*sb(6+oz0+ioz, ilam, ipr))) 

            END DO  !isz
          END DO    !ioz
          flux_x= flux_x + ppr*exp(flux_p)
        END DO      !ipr

        flux(idest)=flux_x 
        idest= idest+1

      END DO        !ilam

      success=.TRUE.
      RETURN

c-----------------------------------------------------------------------
c Handle error conditions
c
  901 WRITE(6,91) ilam0, n_waveln
      GOTO 999
  902 WRITE(6,92) ilam1, n_waveln
      GOTO 999
  903 WRITE(6,93) ilam0, ilam1
  999 CONTINUE
      success=.FALSE.
      RETURN

   91 FORMAT(' !!! fluxinterp: called with ilam0=',i6,'. ilam0 must be ',
     &  'in the range [0, ',i4,'].')
   92 FORMAT(' !!! fluxinterp: called with ilam1=',i6,'. ilam1 must be ',
     &  'in the range [0, ',i4,'].')
   93 FORMAT(' !!! fluxinterp: called with ilam0 > ilam1.'/
     &   '   ilam0=',i6,'  ilam1=',i6)
   94 FORMAT('>gfluxinterp: ',5f10.2,2i5)
c  97 FORMAT(6i5,4(1pe10.2)/5(1pe10.2))
d 990 FORMAT('SZA interpolation precalcs:'/
d    & '  lcs:', 1pe10.2,'  n_sza:',i4,5(/5(1pe10.2)))
d 991 FORMAT('  sz0:', i3,'  szdeg:',i3,'  inrange:', i3,
d    & 5(1pe10.2))
d 992 FORMAT('OZ: ',F7.2,'  TBL: ',10F7.2)
d 993 FORMAT('OZN interpolation precalcs: ',3i3, 4f10.2)

      END
