c
c-----------------------------------------------------------------------
c=======================================================================
c-----------------------------------------------------------------------
c
      subroutine bdnref(sz,vz,rz,chis,chil,chic)
c
c-----------------------------------------------------------------------
c
c PLEASE NOTE: the relative azimuth angle is reversed from its
c normal orientation in these tables, so that when the relative azimuth
c is used in the normal sense (i.e., 0 deg is forward reflection,
c 180 deg is backward reflection), it must first be subtracted from
c 180 deg before calling this subroutine, i.e.
c
c     CALL BDNREF(SZ,VZ,180.-RZ,CHIS,CHIL,CHIC)
c
c-----------------------------------
c
c VIS anisotropic factors for ocean, land, and cloud from solar zenith,
c viewing zenith, and azimuth angles by doing a three way interpolation
c
c     sz = solar zenith angle in degrees  0.=<sz<=90.
c     vz = viewing zenith angle in degrees 0.=<vz<=90.
c     rz = relative azimuth angle in degrees 0.=<rz<=180.
c   chis = anisotropic factor or chi for clear ocean
c   chil = anisotropic factor or chi for clear land
c   chic = anisotropic factor or chi for clouds
c
c-----------------------------------
c
c      s = solar zenith angle bin centers in cosine of angle (ocean,
c          cloud)
c      v = viewing angle bin centers in degrees, also used for solar
c          zenith angle bin centers (land)
c      r = relative azimuth angle bin centers in degrees
c  d1,d2 = interpolated viewing angle factors
c    a() = interpolated azimuth angle factors
c    c() = interpolated solar angle factors
c
c-----------------------------------------------------------------------
c
c Allocate memory
c
c   Arrays
c
      integer iaz(2),ivz(2),icsz(2),isz(2)
      real chi(3),c(2),a(2)
      real v(10),r(13),s(11)
      real surfc(21,10,13),cloud(11,10,13)
c
c   Scalars
c
      logical first
      integer gu
c
c   Save
c
      save first,r,v,s,surfc,cloud
c
c   Data
c
      data first /.true./
      data r /0.,15.,30.,45.,60.,75.,90.,105.,120.,135.,150.,165.,180./
      data v /0.,10.,20.,30.,40.,50.,60.,70.,80.,90./
      data s /1.00,0.95,0.85,0.75,0.65,0.55,0.45,0.35,0.25,0.15,0.05/
c
      gu=1
c
c-----------------------------------------------------------------------
c
c Read GOES bi-directional model on first call only
c
      if (first) then
         open(unit=gu,
     &        file='bdnref.ascii',
     &        form='formatted',
     &        status='old')
         read(gu,100) surfc,cloud
         close(gu)
         first=.false.
      endif
  100 format(10(x,f7.4))
c
c-----------------------------------------------------------------------
c
c Set up angular bins
c
      iaz(1)=int(1.+rz/15.)
      iaz(2)=min0(1+iaz(1),13)
      ivz(1)=int(1.+vz/10.)
      ivz(2)=min0(1+ivz(1),9)
      icsz(1)=int(11.5-10.*cosd(sz))
      icsz(2)=min0(1+icsz(1),11)
      isz(1)=int(1.+sz/10.)
      isz(2)=min0(1+isz(1),10)
c
c-----------------------------------------------------------------------
c
c Interpolate
c
      do kscene=1,3
         do ksz=1,2
            c(ksz)=0.
            do kaz=1,2
               if (kscene.eq.1) then
                  d1=surfc(icsz(ksz),ivz(1),iaz(kaz))
                  d2=surfc(icsz(ksz),ivz(2),iaz(kaz))
               else
                  if (kscene.eq.2) then
                     d1=surfc(11+isz(ksz),ivz(1),iaz(kaz))
                     d2=surfc(11+isz(ksz),ivz(2),iaz(kaz))
                  else
                     d1=cloud(icsz(ksz),ivz(1),iaz(kaz))
                     d2=cloud(icsz(ksz),ivz(2),iaz(kaz))
                  endif
               endif
               if (ivz(1).ne.ivz(2)) then
                  call lineq(v(ivz(1)),v(ivz(2)),d1,d2,slop,cept)
                  a(kaz)=slop*vz+cept
               else
                  a(kaz)=d1
               endif
            enddo
            if (iaz(1).ne.iaz(2)) then
               call lineq(r(iaz(1)),r(iaz(2)),a(1),a(2),slop,cept)
               c(ksz)=slop*rz+cept
            else
               c(ksz)=a(1)
            endif
         enddo
         if (kscene.ne.2) then
            if (icsz(1).ne.icsz(2)) then
               call lineq(s(icsz(1)),s(icsz(2)),c(1),c(2),slop,cept)
               chi(kscene)=slop*cosd(sz)+cept
            else
               chi(kscene)=c(1)
            endif
         else
            if (isz(1).ne.isz(2)) then
               call lineq(cosd(v(isz(1))),cosd(v(isz(2))),c(1),c(2),
     &                                                        slop,cept)
               chi(kscene)=slop*cosd(sz)+cept
            else
               chi(kscene)=c(1)
            endif
         endif
         if (chi(kscene).lt.0)
     &                write(6,*) 'bdnref<0 ',chi(kscene),kscene,sz,vz,rz
      enddo
c
c-----------------------------------------------------------------------
c
c Transfer values and return
c
      chis=chi(1)
      chil=chi(2)
      chic=chi(3)
      return
c
c-----------------------------------------------------------------------
c
c Finem respice
c
      end
c
c-----------------------------------------------------------------------
c=======================================================================
c-----------------------------------------------------------------------
c
      subroutine lineq(x1,x2,y1,y2,a,b)
      real x1,x2,y1,y2,a,b
      if ((x2-x1).eq.0.) stop
      a=(y2-y1)/(x2-x1)
      b=y1-a*x1
      return
      end