SUBROUTINE interp_s(ik,levels,sza,o3column,s,jo2, 1,1
     $                    sdat,o2jdat,sza_tab,o3_tab)

      USE SC_data, ONLY : nlam,nsza,numo3

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: ik,levels
      REAL, INTENT(IN) :: sza,o3column(levels)  
      REAL (KIND=4), INTENT(IN) :: sdat(nsza,numo3,levels,nlam)
      REAL (KIND=4), INTENT(IN) :: o2jdat(nsza,numo3,levels)
      REAL (KIND=4), INTENT(IN)  :: sza_tab(nsza)
      REAL (KIND=4), INTENT(IN)  :: o3_tab(numo3,levels)
      REAL, INTENT(OUT) :: s(nlam), jo2

      INTEGER :: ijj,ikk,ikkm,il,is
      REAL :: omt,omu,t,u
C+
C NAME:
C   interp_s
C PURPOSE:
C   Interpolate s values for each wavelength in table to specified O3
c   col and zenith angles
C CATEGORY:
C CALLING SEQUENCE:
C   Call interp_s(nlam,sza,o3column,s)
C INPUTS:
C   nlam -- number of wavelength intervals used
C   sza --  input array (lat) zenith angle values for which s is desired
C   o3column --input array (lat,levels) overhead o3 column value 
C OPTIONAL INPUT PARAMETERS:
C OUTPUTS:
C   s -- array of s values (nlam,levels,lat) for each wavelength 
C             at model p-level interpolated to o3column and sza values
C   jo2 -- array (lat,levels) of J(O2) values interped as above (REAL)
C INTERNAL VARIABLES
C   sza_tab -- values of sza corresponding to sdat table values
C   o3_tab -- array of overhead O3 values at each p-level (numo3s,np_ctm)
C                used to index sdat
C   sdat -- input array of values of radiative source function 
C           (nzens,numo3,np_ctm,nlam) gridded to ctm p layers
C   o2jdat -- array (nzens,numo3s,levels) of J(O2) values from srcfun output
C COMMON BLOCKS:
C SIDE EFFECTS:
C PROCEDURE:
C   bi-linear interpolation, for sza>94 s=0, for o3 out of range use min/max
C RESTRICTIONS:
C REQUIRED ROUTINES:
C MODIFICATION HISTORY: 
C   Created 930825 - SR Kawa
C   Modified 960710 for 28 levels and to handle J(O2) separately
C
       
c... For each input solar zenith angle, find the first element of
c... tabled sza_tab values that is greater than it.  Use this
c... table element and previous table element to determined
c... interpolated value.
      do is=1,nsza
         ijj = is 
         if(sza_tab(is).gt.sza) goto 333 
      enddo 
 333  continue
      
c... point is dark, set s/jo2=0        
      if(sza.gt.sza_tab(nsza)) then
         do il=1,nlam
            s(il) = 0. 
         enddo 
         jo2 = 0.
      else  
         t = (sza-sza_tab(ijj-1))/(sza_tab(ijj)-sza_tab(ijj-1))
         omt = 1.-t
         
c...  For each input overhead o3 column find the first element
c...  of tabled o3_tab values that is > than it.  Use this
c...  table element and previous table element to determine
c...  interpolated value
         do is=1,numo3
            ikk = is 
            if (o3_tab(is,ik).gt.o3column(ik)) goto 334 
         enddo 
 334     continue
         ikkm = ikk-1 
         if((ikk.gt.1).and.(o3column(ik).le.o3_tab(numo3,ik))) then 
            u = (o3column(ik)-o3_tab(ikkm,ik))/
     &           (o3_tab(ikk,ik)-o3_tab(ikkm,ik))
            omu = 1.-u
c...  do bilinear interpolation at ik for each wavelength
c...  from numerical recipes, p.96
            do il=1,nlam       
               s(il) = omt*omu*sdat(ijj-1,ikkm,ik,il)
     &              +t*omu*sdat(ijj,ikkm,ik,il)
     &              +t*u*sdat(ijj,ikk,ik,il)
     &              +omt*u*sdat(ijj-1,ikk,ik,il)
            enddo
            jo2 = omt*omu*o2jdat(ijj-1,ikkm,ik)
     &           +t*omu*o2jdat(ijj,ikkm,ik)
     &           +t*u*o2jdat(ijj,ikk,ik)
     &           +omt*u*o2jdat(ijj-1,ikk,ik)       
c...  extrap before table
         else if (ikk.eq.1) then
            do il=1,nlam
               s(il) = omt*sdat(ijj-1,1,ik,il)+t*sdat(ijj,1,ik,il)
            enddo
            jo2 = omt*o2jdat(ijj-1,1,ik)+t*o2jdat(ijj,1,ik)
c...  extrap past table
         else 
            do il=1,nlam
               s(il) = omt*sdat(ijj-1,numo3,ik,il)+
     &              t*sdat(ijj,numo3,ik,il)
            enddo 
            jo2 = omt*o2jdat(ijj-1,numo3,ik)+t*o2jdat(ijj,numo3,ik)
         endif  
      endif
      
      RETURN
      END SUBROUTINE interp_s