      subroutine albval2(alb,ib,tau_in,csz_in,wind_in,chl_in,
     1                   nstart,nend,nval)
c***********************************************************************
c J.Cole    ...April 17, 2018.  Limit the cosine of solar zenith angle to 
c                               max and min in table.
c M.Lazare  ...Jul 20, 2014. Modified to pass in start and stop indices.
c L.Solheim ...Apr 2,2004.   Prevous version albval.
c
c Given nval values for optical depth (tau_in), cosine of solar zenith
c angle (csz_in), wind speed (wind_in) and chlorophyll concentration
c (chl_in), return nval surface albedo values for band ib. These albedo
c values are determined using a lookup table produces by Z. Jin and
c Thomas Charlock at NASA.
c
c Quadra-linear interpolation is done over the nodal values for tau,
c csz, wind and chl. Extrapolation is done on values outside of the
c nodal range. This produces reasonable values as long as the point
c is "close" to an end point and is efficient because no checks need
c to be made on out of range values. However, linear extrapolation is
c generally not good for points not "close" to the end points. Therefore
c care should be taken so that this routine is not called with
c interpolation points that fall outside of the nodal point range.
c
c Input:
c  ib            ...frequency band for which albedos are to be evaluated
c  tau_in (nval) ...optical depths
c  csz_in (nval) ...cosines of solar zenith angle
c  wind_in(nval) ...wind speeds (m/s)
c  chl_in (nval) ...chlorophyll concentrations (mg/m**3)
c  nstart        ...starting index number.
c  nend          ... end     index number
c  nval          ...number of albedo values requested
c
c Output:
c  alb(nval) ...surface albedo values corresponding to input values of
c               optical depth, zenith angle, wind speed and chlorophyll
c               concentration for band ib
c***********************************************************************

      !--- input/output parameters
      integer       :: ib,nval,nstart,nend
      real,  dimension(nval) :: tau_in, csz_in, wind_in, chl_in
      intent(in)  :: ib, tau_in, csz_in, wind_in, chl_in
      real,  dimension(nval) :: alb
      intent(out) :: alb

      !--- local work space
      real   wt(2), ws(2), ww(2), wc(2), csz_loc(nval)

      !--- nodal values found in lookup table

      !--- table dimensions
      ! nb ...number of bands
      !       4 currently: 0.20-0.69, 0.69-1.19, 1.19-2.38, 2.38-4.00 um
      ! nt ...number of aerosol/cloud optical depths
      ! ns ...number of solar zenith angles
      ! nw ...number of wind speeds
      ! nc ...number of chlorophyll concentration values
      parameter (nb=4, nt=16, ns=15, nw=7, nc=5)

      !--- aerosol/cloud optical depth: tau<1.0 aerosol, tau>1.0 cloud
      real,   parameter, dimension(nt) :: tau=
     & (/ 0.00, 0.05, 0.10, 0.16, 0.24, 0.35,  0.5,  0.7,
     &    0.99, 1.30, 1.80, 2.50, 5.00, 9.00, 15.0, 25.0 /)

      !--- cosine of solar zenith angle
      real,   parameter, dimension(ns) :: csz=
     & (/ 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 /)

      !--- wind speeds in m/s
      real,   parameter, dimension(nw) :: wind=
     & (/ 0.0, 3.0, 6.0, 9.0, 12.0, 15.0, 18.0 /)

      !--- ocean chlorophyll concentrations in mg/m**3
      real,   parameter, dimension(nc) :: chl=
     & (/ 0.0, 0.1, 0.5, 2.0, 12.0 /)

      !--- common space for lookup table data
      common /salbtab/ atable(nb,nc,nw,ns,nt), init_table
c

c.....Abort if atable has not been initialized
      if (init_table.ne.1) then
        write(6,*)'albval: Lookup table has not been initialized'
        call xit('ALBVAL',-1)
      endif

c.....Ensure ib is valid
      if (ib.lt.1.or.ib.gt.nb) then
        write(6,*)'albval: ib is out of range. ib = ',ib
        call xit('ALBVAL',-2)
      endif

c....Bound input csz to be within min and max of lookup table
      do i=nstart,nend
         csz_loc(i) = MAX(MIN(csz_in(i),csz(ns)),csz(1))
      end do

      do i=nstart,nend
c.......Calculate weigths
        it=mvidx(tau,nt,tau_in(i))
        is=mvidx(csz,ns,csz_loc(i))
        iw=mvidx(wind,nw,wind_in(i))
        ic=mvidx(chl,nc,chl_in(i))
        wt(2) = (tau_in(i)-tau(it))/(tau(it+1)-tau(it))
        wt(1) = 1.0-wt(2)
        ws(2) = (csz_loc(i)-csz(is))/(csz(is+1)-csz(is))
        ws(1) = 1.0-ws(2)
        ww(2) = (wind_in(i)-wind(iw))/(wind(iw+1)-wind(iw))
        ww(1) = 1.0-ww(2)
        wc(2) = (chl_in(i)-chl(ic))/(chl(ic+1)-chl(ic))
        wc(1) = 1.0-wc(2)

c.......calculate alb(i) for band ib, using quadra-linear interpolation
        alb(i)=0.0
        do itt=it,it+1
          do iss=is,is+1
            do iww=iw,iw+1
              do icc=ic,ic+1
                wtt=wt(itt-it+1)*ws(iss-is+1)*ww(iww-iw+1)*wc(icc-ic+1)
                alb(i) = alb(i) + wtt*atable(ib,icc,iww,iss,itt)
              enddo
            enddo
          enddo
        enddo
      enddo
      return
      end
