C  This is an example to read the ocean albedo look-up-table (ocnalbtab1bnd.bin).
C  
	parameter (nb=4, nt=16, ns=15, nw=7, nc=5)
        real*4 alb
	real   taunode(16), szanode(15), windnode(7), chlnode(5),
     $         albb,
     $         rt(2), rs(2), rw(2), rc(2)

	data taunode /0.0, 0.05, 0.1, 0.16, 0.24, 0.35, 0.5, 0.7, 0.99,
     &                1.3, 1.8, 2.5, 5.0, 9.0, 15.0, 25.0 /
        data szanode /0.05, 0.09, 0.15, 0.21, 0.27, 0.33, 0.39, 0.45,
     &                0.52, 0.60, 0.68, 0.76, 0.84, 0.92, 1.0 /
        data windnode /0., 3., 6., 9., 12., 15., 18. /
        data chlnode /0.0, 0.1, 0.5, 2.0, 12.0/

c    specify the parameters for albedo here:

	tau = 10.00		!aerosol/cloud optical depth
	sza = 0.100		!cosine of solar zenith angle
	wind = 5.00		!wind speed in m/s
	chl = 0.10		!chlorophyll concentration in mg/m3
c
C   Output albedo for these input is 0.054. If not, your system may use
C   different record unit and you need change the record length.

C    now find the albedo corresponding to the 4 parameters above:

C    !!!Note: if the record unit in your sestem is WORDS instead of BYTE,
C    !!!      you MUST change the "recl" (record length) in the following
C    !!!      OPEN statement from 4 to 1 (recl=1).

        open(1,file='ocnalbtab1bnd.bin'
     &  ,form='Unformatted',access='Direct',recl=4)

	if(tau.lt.0. .or. (sza.lt.0. .or. sza.gt.1.) .or. wind.lt.0.
     &     .or. chl.lt.0.)stop 'Err: input parameters wrong!'
	if(tau .gt. taunode(nt))tau=taunode(nt)
	if(wind .gt. windnode(nw))wind=windnode(nw)
	if(chl .gt. 45.)chl=45.
	call locate(taunode,nt,tau,it)
	call locate(szanode,ns,sza,is)
	call locate(windnode,nw,wind,iw)
	call locate(chlnode,nc,chl,ic)
	rt(2) = (tau-taunode(it))/(taunode(it+1)-taunode(it))
        rt(1) = 1.0-rt(2)
	rs(2) = (sza-szanode(is))/(szanode(is+1)-szanode(is))
        rs(1) = 1.0-rs(2)
	rw(2) = (wind-windnode(iw))/(windnode(iw+1)-windnode(iw))
        rw(1) = 1.0-rw(2)
	rc(2) = (chl-chlnode(ic))/(chlnode(ic+1)-chlnode(ic))
        rc(1) = 1.0-rc(2)

C             ** get albedo by 4 dimensional linear interpolaton **
	albb = 0.0
	do 100 itt=it,it+1
	do 100 iss=is,is+1
	do 100 iww=iw,iw+1
	do 100 icc=ic,ic+1
	  irec = (itt-1)*15*7*5 + (iss-1)*7*5 + (iww-1)*5 + icc
	  read(1,rec=irec) alb
	  wtt = rt(itt-it+1)*rs(iss-is+1)*rw(iww-iw+1)*rc(icc-ic+1)
	  albb = albb + wtt*alb
100	continue
	print*,"rec=",irec,alb

	write(*,'(4a6,/,4f6.2,/,a,f6.3,/)')'Tau','COSUN','Wind',
     &        'Chl',tau,sza,wind,chl,
     &        'SW Albedo(0.25-4.0um) :',albb
	
        close(1)
	stop
	end

c=======================================================================
      subroutine locate(xx,n,x,j)
c
c purpose:  given an array xx of length n, and given a value X, returns
c           a value J such that X is between xx(j) and xx(j+1). xx must
c           be monotonic, either increasing of decreasing. this function
c           returns j=1 or j=n-1 if x is out of range.
c
c input:
c   xx      monitonic table
c   n       size of xx
c   x       single floating point value perhaps within the range of xx
c
      integer j,n
      real x,xx(n)
      integer jl,jm,ju

      if(x.eq.xx(1)) then
        j=1
        return
      endif
      if(x.eq.xx(n)) then
        j=n-1
        return
      endif
      jl=1
      ju=n
10    if(ju-jl.gt.1) then
        jm=(ju+jl)/2
        if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      goto 10
      endif
      j=jl
      return
      end
c=======================================================================