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