SUBROUTINE RFLUXINTERP

C   Received from Randy Kawa 1/24/95 - adapted from photcodes7.f 

#include "comphot.h"

      	real tempuse(nlat,nz,nx,nlam),tuse(nlat)
        integer locsza(nlat)
 
c       interpolate radiative flux function values to model conds. 

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 ij=1,nlat
           do is=1,nsza
              if (sza_tab(is).gt.sza2d(ij)) go to 333 
           end do 
 333       locsza(ij) = is 
        end do
            
        do ij=1,nlat
           if (locsza(ij).eq.(nsza+1)) then 
c             Point is in darkness.  Set s's to zero. 
              do ik=1,nz 
                 o2jint(ij,ik) = 0.
              do il=1,nlam
                 sflux(ij,ik,il) = 0. 
              end do 
              end do 
           else  
              ijj = locsza(ij)                       
              tuse(ij)=(sza2d(ij)-sza_tab(ijj-1))/
     1         (sza_tab(ijj)-sza_tab(ijj-1)) 
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 ik=1,nz 
                 do is=1,no3 
                    if (o3_tab(is,ik).gt.ovho3(ij,ik)) go to 334 
                 end do 
 334             ikk = is 
                 ikkm = ikk-1 
                 if ((ikk.gt.1).and.(ikk.le.no3)) then 
                    u=(ovho3(ij,ik)-o3_tab(ikkm,ik))/
     1               (o3_tab(ikk,ik)-o3_tab(ikkm,ik))          

c                   do bilinear interpolation at ik for each wavelength
c                   from numerical recipes, p.96

	            do il=1,nlam       
	               sflux(ij,ik,il)=(1.-tuse(ij))*(1.-u)*
     1                  stab(ijj-1,ikkm,ik,il)+tuse(ij)*(1.-u)*
     1		        stab(ijj,ikkm,ik,il)+tuse(ij)*u*
     1                  stab(ijj,ikk,ik,il)+
     1	         	(1.-tuse(ij))*u*stab(ijj-1,ikk,ik,il)
                    end do

	    o2jint(ij,ik)=(1.-tuse(ij))*(1.-u)*o2jdat(ijj-1,ikkm,ik)
     1			+tuse(ij)*(1.-u)*o2jdat(ijj,ikkm,ik)
     2			+tuse(ij)*u*o2jdat(ijj,ikk,ik)
     3	         	+(1.-tuse(ij))*u*o2jdat(ijj-1,ikk,ik)
                                
                 else if (ikk.eq.1) then 
c                   write (6,33) ij,ik,ovho3(ij,ik),o3_tab(1,ik)
 33                 format(' Ovhd o3 col(',i3,i3,') of ',e10.3,
     1               ' too thin! Lowest tabled val=',e10.3)  
      
                    do il=1,nlam
	             sflux(ij,ik,il)=(1.-tuse(ij))*stab(ijj-1,1,ik,il)+
     1                  tuse(ij)*stab(ijj,1,ik,il)
	            end do        
	        o2jint(ij,ik)=(1.-tuse(ij))*o2jdat(ijj-1,1,ik)+
     1                  tuse(ij)*o2jdat(ijj,1,ik)

                 else 
c                    write (6,34) ij,ik,ovho3(ij,ik),o3_tab(no3,ik)
  34                 format(' Ovhd o3 col(',i3,i3,') of ',e10.3,
     1               ' too thick! Highest tabled val=',e10.3) 
                    do il=1,nlam
	          sflux(ij,ik,il)=(1.-tuse(ij))*stab(ijj-1,no3,ik,il)+
     1                  tuse(ij)*stab(ijj,no3,ik,il)
	            end do 
	          o2jint(ij,ik)=(1.-tuse(ij))*o2jdat(ijj-1,no3,ik)+
     1                  tuse(ij)*o2jdat(ijj,no3,ik)

                 end if  
              end do 
	    end if
           end do 

        return
        end