function cldtau(sza,az,sat,pt,rs,rad)
c
c author: colin seftor (from ed celarier's taucloud)
c
c date:   11 january 2002
c
c purpose:
c   **** use linear interpolation to ****
c   calculate the putative cloud optical thickness for a particular
c   viewing geometry, and for a terrain of a certain pressure and
c   reflectivity, given the toms-measured radiance in one of the 
c   reflectivity channels.
c
c calling parameters:
c   name	type	i/o	units	description
c   --------	----	---	-----	--------------
c   cldtau 	r8	o		calculated putative cloud optical
c					thickness
c   sza		r8	i	deg	solar zenith angle
c   az		r8	i	deg	solar-fov-satellite azimuth
c   sat		r8	i	deg	satellite zenith angle
c   pt		r8	i	atm	terrain surface pressure
c   rs		r8	i	frac	surface reflectivity
c   rad		r8	i	i/f	toms-measured radiance
c
c external routines required:
c   radlow      calculates pairs (tcl, radiance) for low tcl.
c   radhigh     calculates pairs (tcl, radiance) for high tcl.
c
c Added some lines to take into account the case there 
c rad falls in between the radhigh and radlow tables
c Aapo Tanskanen Feb 4, 2003
c 
c cldtau should not be interpolated in case rad is higher 
c than the upper boundary of radhigh, ie. cldtau=100
c Both tables need to be retrieved
c Aapo Tanskanen Feb 5, 2003
c
c============================================================================
      implicit none
      integer*4 ntcl_l
      parameter(ntcl_l=5)
      integer*4  ntcl_h
      parameter(ntcl_h=6)
c     -- calling parameters
      real*8 cldtau,sza,az,sat,pt,rs,rad
c     -- local variables
      real*8 outrad(0:ntcl_l+ntcl_h-1)
      real*8 outrad_l(0:ntcl_l-1), outrad_h(0:ntcl_h-1)
      real*8 outtcl_l(0:ntcl_l-1), outtcl_h(0:ntcl_h-1)
      integer*4 i0, nout_l, nout_h, i
      real*8 frac
      logical success
c
c     -- get radiances for given solar, azimuth, and satellite angles
      call radlow(rs,sza,az,sat,pt,outrad_l,outtcl_l,nout_l,success)
      if(.not.success) return
      call radhigh(sza,az,sat,rs,outrad_h,outtcl_h,nout_h,success)
      if (.not.success) return
      if (rad .le. outrad_l(nout_l-1)) then
        if (rad .le. outrad_l(0)) then
          cldtau= 0.d0
          return
        else
          do i = 0,nout_l-1
            if ((rad.gt.outrad_l(i) .and. rad.le.outrad_l(i+1)) .or.
     1          (rad.lt.outrad_l(i) .and. rad.ge.outrad_l(i+1))) then
              frac = (rad - outrad_l(i)) / (outrad_l(i+1) - outrad_l(i))
              cldtau = (1.-frac)*outtcl_l(i) + frac*outtcl_l(i+1)
              goto 100
            endif
          enddo
100       continue
        endif
      else
        if (rad .gt. outrad_h(nout_h-1)) then
          cldtau= 100.d0
          return
        else 
          do i = 0,nout_h-1
            if ((rad.gt.outrad_h(i) .and. rad.le.outrad_h(i+1)) .or.
     1          (rad.lt.outrad_h(i) .and. rad.ge.outrad_h(i+1))) then
              frac = (rad - outrad_h(i)) / (outrad_h(i+1) - outrad_h(i))
              cldtau = (1.-frac)*outtcl_h(i) + frac*outtcl_h(i+1)
              goto 200
            endif
          enddo
200       continue
        endif
      endif
      if ((rad.gt.outrad_l(nout_l-1) .and. rad.lt.outrad_h(0)) .or.
     1     (rad.lt.outrad_l(nout_l-1) .and. rad.gt.outrad_h(0))) then
         frac = (rad - outrad_l(nout_l-1)) / 
     1        (outrad_h(0) - outrad_l(nout_l-1))
         cldtau = frac*outtcl_h(0) + (1.-frac)*outtcl_l(nout_l-1)
      endif
      end