c=======================================================================


      subroutine phaers(nwa,wbaer,pmaer) 1,1
c
c purpose
c     save legendre moments of scattering phase function of 
c     boundary layer aerosols.  A single phase function is used
c     for all boundary layer aerosols
c
c input
c     nstr      number of radiation streams
c     nwa       number of wavelength bins
c     wbaer     wavelength at which phase information is available
c     pmaer     legendre moments of phase function
c
      parameter (naerw=47)
      parameter (nstrms=40)
      real wbaer(*),pmaer(*),wb(naerw),pm(naerw*nstrms),p(0:nstrms)

      save nw,wb,pm
      
      nw=nwa

      do i=1,nw
        wb(i)=wbaer(i)
      enddo

      nn=nstrms*nw
      do i=1,nn
        pm(i)=pmaer(i)
      enddo
c      write(*,'(a,/,1p(8e11.3))') 'wb',(wb(i),i=1,8)
c      write(*,'(a,/,1p(8e11.3))') 'pm',(pm(i),i=1,8*nstrms)

      return

c-----------------------------------------------------------------------
      entry phaerw(w,ta,tb,p)
c
c purpose:
c   Interpolate phase function moments to wavelength w.
c   Must use linear interpolation because moments can be pos or neg.
c   if w>wmax then p(w)=ptable(wmax)
c   if w<wmin then p(w)=ptable(wmin)
c
c input
c   w      wavelength in microns
c   ta     total aerosol optical depth increment in layer
c   tb     boundary layer aerosol optical depth increment in layer
c   p      legendre moments of scat phase function due to all 
c          aerosols except boundary layer aerosols
c output
c   p      legendre moments of scattering phase function (0:nstrms)
c          due to all aersols
c
      call locate(wb,nw,w,l)
      wt=(w-wb(l))/(wb(l+1)-wb(l))
      wt=max(0.,min(wt,1.))
      if(ta.gt.0.) then
        wa=tb/ta
      else
        wa=1.
      endif

c      write(*,'(8a11)') 'l','w','wt','wa','wb(l)','wb(l+1)','ta','tb'
c      write(*,'(i11,1p8e11.3)') l,w,wt,wa,wb(l),wb(l+1),ta,tb
c      write(*,'(a,/,1p(10e11.3))') 'pmaer1',
c     &     (pm(l+i*nw),i=1,nstrms)
c      write(*,'(a,/,1p(10e11.3))') 'pmaer2',
c     &     (pm(l+1+i*nw),i=1,nstrms)

      do i=1,nstrms
        pin=p(i)
        j=l+(i-1)*nw
        p(i)=pm(j)*(1.-wt)+pm(j+1)*wt
        p(i)=pin*(1.-wa)+p(i)*wa
      enddo
      p(0)=1.

c      write(*,'(a,/,1p(10e11.3))') 'pmaer3',
c     &     (p(i),0=1,nstrms)

      return
      end