C  This is an example to use the ocean albedo look-up-table (ocnalbtab.bin) 
C   to obtain the albedo at your specified spectral band, optical depth,
C   cosine of solar zenith, wind and chlorophyll concentration.
C  
	parameter (nb=24, nt=16, ns=15, nw=7, nc=5)
        real*4 alb(24),rflx(24)
	real   taunode(16), szanode(15), windnode(7), chlnode(5),
     $         wlnode(25), albb(24),wtflx(24),
     $         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/
        data wlnode /0.25, 0.30, 0.33, 0.36, 0.40, 0.44, 0.48, 0.52,
     &               0.57, 0.64, 0.69, 0.752, 0.780, 0.867, 1.0, 1.096,
     &               1.19, 1.276, 1.534, 1.645, 2.128, 2.381, 2.907, 
     &               3.425, 4.0/

C     --------------------------- INPUT ---------------------------------------
C   
c    specify the parameters for albedo here:

        tau = 0.40              !aerosol/cloud optical depth
        sza = 0.50             !cosine of solar zenith angle
        wind = 10.00             !wind speed in m/s
        chl = 0.10              !chlorophyll concentration in mg/m3
        wls = 0.69              !start wavelength (um) of your band
        wle = 1.19              !end wavelength (um) of your band
c
C   Output albedo for these input is 0.068. If not, your system may use
C   different record unit and you need change the record length.
C     -------------------------------------------------------------------------

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 24*4 to 24 (recl=24).	

        open(1,file='ocnalbtab24bnd.bin'
     &  ,form='Unformatted',access='Direct',recl=24*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)

	if(wls .gt. wle)stop 'Err: Start wavelength should be smaller.'
	if(wls .lt. wlnode(1))wls=wlnode(1)
	if(wle .gt. wlnode(25))wle=wlnode(25)
	call locate(wlnode,25,wls,ib1)
	call locate(wlnode,25,wle,ib2)

C             ** get alb(ib1:ib2) by 4 dimensional linear interpolaton **
	do ib=ib1,ib2
	  albb(ib) = 0.0
	  wtflx(ib) = 0.0
	enddo
	do 100 itt=it,it+1
	do 100 iss=is,is+1
	 do 200 iww=iw,iw+1
	 do 200 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)
	  do ib=ib1,ib2
	    albb(ib) = albb(ib) + wtt*alb(ib)
	  enddo
200	 continue
c	          *** get 24 band down flux weights ***
	 irr = 8400 + (itt-1)*15 + iss
         read(1,rec=irr) rflx
         wtf = rt(itt-it+1)*rs(iss-is+1)
c	 do ib=ib1,ib2
	 do ib=1,24
	   wtflx(ib) = wtflx(ib) + wtf*rflx(ib)
	 enddo
100	continue

C             ** get albedo in the specified band by weighted sum **
	twtb = 0.0
	albedo = 0.0
	do ib=ib1,ib2
	  if(ib .eq. ib1)then
	    wtb = (wlnode(ib1+1)-wls)/(wlnode(ib1+1)-wlnode(ib1))
	  else if (ib .eq. ib2)then
	    wtb = (wle-wlnode(ib2))/(wlnode(ib2+1)-wlnode(ib2))
	  else
	    wtb = 1.0
	  endif
     	  albedo = albedo + wtb*wtflx(ib)*albb(ib)
	  twtb = twtb + wtb*wtflx(ib)
	enddo
	if (twtb.eq.0. .and. ib1.eq.ib2)then
	   albedo = albb(ib1)
	else
	   albedo = albedo/twtb
	endif

	write(*,'(6a6,/,6f6.2,/,a11,f6.3,/)')'Tau','COSUN','Wind',
     &        'Chl','WL1','WL2',tau,sza,wind,chl,wls,wle,
     &        'Albedo =',albedo
	
        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=======================================================================