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

      subroutine aerbwi(wl,extinc,wa,ga) 2,9
c 
c purpose
c     At a given wavelength find the extinction efficiency,
c     single scattering albedo and asymmetry factor of boundary
c     layer aerosols by interpolation
c
c input
c  wl         wavelength
c 
c output
c  extinc     extinction efficiency
c  wa         single scattering albedo
c  ga         asymmetry factor
c  pa         legendre moments of scattering phase function 

      
      parameter (nstrms=40)
      parameter (naerw=47)
      parameter (naerwm=naerw*nstrms)
      real*8 zip
      parameter (zip=-1.d0)

      character*80 errmes(10)

      real wlbaer(naerw),qbaer(naerw),wbaer(naerw),gbaer(naerw)
      real pmaer(*)
c
      common /extd/ awl(naerw),dummy(naerw,48),aerstr(naerw,3,4)
c
      common /aerblk/
     &    wlaer(naerw),aerext(naerw),aerabs(naerw),aerasm(naerw),beta
      
      data maerw/naerw/
      save maerw

c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

      !write(*,*) 'db_aerbwi.1'

      if(wlaer(1).eq.0.) then   ! no boundary layer aerosol
        extinc=0.
        wa=0.
        ga=0.
        return
      endif        

      call locate(wlaer,maerw,wl,l)
      if(wl.le.wlaer(1)) then 
        extinc=aerext(1)*(wlaer(1)/wl)**beta
        wa=1.-(aerabs(1)/aerext(1))
        ga=aerasm(1)
      elseif(wl.ge.wlaer(maerw)) then 
        extinc=aerext(maerw)*(wlaer(maerw)/wl)**beta
        wa=1.-(aerabs(maerw)/aerext(maerw))
        ga=aerasm(maerw)
      else
        wt=log(wl/wlaer(l))/log(wlaer(l+1)/wlaer(l))

        extinc=aerext(l)*(aerext(l+1)/aerext(l))**wt
        if(aerabs(l).gt.0. .and. aerabs(l+1).gt.0. ) then 
          absorp=aerabs(l)*(aerabs(l+1)/aerabs(l))**wt
        else
          absorp=aerabs(l)*(1.-wt)+aerabs(l+1)*wt
        endif
        if(extinc.gt.0.) wa=max(0.,min(1.-absorp/extinc,1.))
        ga=(1.-wt)*aerasm(l)+wt*aerasm(l+1)  
        ! linear interpolation because aerasm can be negative
      endif
      !write(*,*) 'wlaer: ',wlaer(1),wlaer(2)
      !write(*,*) 'extinc: ',extinc
      !write(*,*) 'aerext: ',aerext(1),aerext(2)
      !write(*,*) 'wl,l ',wl,l

      return

c.......................................................................

      entry aerswi(ja,wl,extinc,wa,ga)

c 
c purpose
c     compute extinction efficiency, single scattering albedo and
c     asymmetry factor of STRATOSPHERIC aerosols by interpolation
c
c input
c  ja         stratospheric aerosol type
c  wl         wavelength
c 
c output
c  extinc     extinction efficiency of stratospheric aerosol of type ja
c  wa         single scattering albedo
c  ga         asymmetry factor
c

      call locate(awl,naerw,wl,l)
      if(wl.le.awl(1)) then 
        extinc=aerstr(1,1,ja)*(awl(1)/wl)**beta
        wa=1.-(aerstr(1,2,ja)/aerstr(1,1,ja))
        ga=aerstr(l,3,ja)
      elseif(wl.ge.awl(maerw)) then 
        extinc=aerstr(naerw,1,ja)*(awl(1)/wl)**beta
        wa=1.-(aerstr(naerw,2,ja)/aerstr(naerw,1,ja))
        ga=aerstr(naerw,3,ja)
      else
        wt=log(wl/awl(l))/log(awl(l+1)/awl(l))
        extinc=aerstr(l,1,ja)*(aerstr(l+1,1,ja)/aerstr(l,1,ja))**wt
        absorp=aerstr(l,2,ja)*(aerstr(l+1,2,ja)/aerstr(l,2,ja))**wt
        if(extinc.gt.0.) wa=max(0.,min(1.-absorp/extinc,1.))
        ga=(1.-wt)*aerstr(l,3,ja)+wt*aerstr(l+1,3,ja)  
        ! linear interpolation because aerasm can be negative
      endif
      !write(*,*) 'awl: ',awl(1),awl(2)
      !write(*,*) 'extinc: ',extinc
      !write(*,*) 'aerext: ',aerext(1),aerext(2)
      !write(*,*) 'wl,l ',wl,l

      return

c.......................................................................

      entry usraer(wlbaer,qbaer,wbaer,gbaer,pmaer,abaer,q55,imoma)
c
c purpose 
c       Replaces data in wlaer,aerext,aerabs,aerasm with user specified
c       aerosol spectral data.  If pmaer is defined, it is saved for
c       later use in PHAERW
c
c       A similar function is performed by subroutine STDAER, which
c       fills the same variables with spectral properties of standard
c       aerosol models
c
c input:
c
c  wlaer    wavelengths at which aerosol model is specified
c  qbaer    aerosol extinction efficiency
c  wbaer    aerosol single scattering albedo
c  gbaer    aerosol asymmetry factor
c  abaer    power law slope used to extrapolate optical depth
c           outside of wlaer range  tauaer=(wlaer(j)/wl)**abaer
c           where j is 1 or n (n = total number of assigned values
c           of wlaer).  Power law extrapolation is not used when
c           ABAER.le.0
c  pmaer    user defined aerosol phase function (may be undefined)
c
c input/output:
c
c  q55      extinction efficiency at 0.55 um
c
c output:
c
c  imoma    0 => aerosol phase function not defined (use Henyey-Greenstein)
c           1 => user defined aersol phase function
c

c check input parameters

      beta=max(abaer,0.)

      nwlbaer=numset(zip,wlbaer,naerw)
      nqbaer=numset(zip,qbaer,naerw)
      nwbaer=numset(zip,wbaer,naerw)
      ngbaer=numset(zip,gbaer,naerw)
      npmaer=numset(zip,pmaer,naerwm)

      ne=0

      if(nwlbaer.eq.0) then         ! wavelength points not set

        qbaer(1)=1.
        nqbaer=1
        if(nwbaer.ne.1) then
          ne=ne+1
          errmes(ne)='specify one value of wbaer when wlbaer not set'
        endif

      elseif(nwlbaer.eq.1) then     ! only one wavelength point

        if(nqbaer.gt.1) then 
          ne=ne+1
          errmes(ne)='number of elements must match: wlbaer, qbaer'
        endif          

      else                          ! many wavelength point

        if(nwlbaer.ne.nqbaer) then
          ne=ne+1
          errmes(ne)='number of elements must match: wlbaer, qbaer'
        endif
        if(nwlbaer.ne.nwbaer) then
          ne=ne+1
          errmes(ne)='number of elements must match: wlbaer, wbaer'
        endif
        if(ngbaer.eq.0) then 
          if(npmaer.ge.1) then 
            if(mod(npmaer,nwlbaer).ne.0) then
              ne=ne+1
              errmes(ne)='incorrect number of phase function moments'
            endif
          else
            ne=ne+1
            errmes(ne)='must specify either gbaer or pmaer'
          endif
        else
          if(ngbaer.ne.nwbaer) then
            ne=ne+1
            errmes(ne)='number of elements must match: wlbaer, gbaer'
          endif
        endif
      
      endif

      if((ngbaer.eq.0).eqv.(npmaer.eq.0)) then 
        ne=ne+1
        errmes(ne)='must specify either gbaer or pmaer, not both'
      endif
       
      if(ne.gt.0) then 
        write(*,*) 'Error in user specified aerosols (iaer=5)'
        write(*,'(/,1x,5a8)')
     &          'nwlbaer','nqbaer','nwbaer','ngbaer','npmaer'
        write(*,'(5i8,/)')  nwlbaer,nqbaer,nwbaer,ngbaer,npmaer
        do i=1,ne
          write(*,'(2a)') 'Error in USRAER -- ',errmes(i)
        enddo
        stop
      endif
      
c duplicate input to two wavelength points if a 
c single wavelength specified in the input file

      !write(*,*) 'db usraer 1'
      if(nwbaer.eq.1) then  
        if(nwlbaer.eq.0) wlbaer(1)=.55 
        wlbaer(2)=2*wlbaer(1)
        nwlbaer=2
        if(abaer.ge.0.) then
          qbaer(2)=qbaer(1)*(wlbaer(1)/wlbaer(2))**beta
        else
          qbaer(2)=qbaer(1)
        endif
        wbaer(2)=wbaer(1)
        gbaer(2)=gbaer(1)

c duplicate phase function entries

        !write(*,*) 'db usraer 2'
        if(npmaer.gt.0) then
          do i=npmaer,1,-1
            pmaer(2*i)=pmaer(i)
            pmaer(2*i-1)=pmaer(i)
          enddo
          npmaer=2*npmaer
        endif
        nqbaer=2
      endif

      !write(*,'(a10,/,1p10e11.3)')  'wlbaer',(wlbaer(i),i=1,2)
      !write(*,'(a10,/,1p10e11.3)')  'qbaer',(qbaer(i),i=1,2)
        
      !write(*,*) 'db usraer 3'
      if(nwlbaer.gt.naerw) then
        write(*,*) 'Too many spectral points in usraer'
        stop
      endif
      do i=1,nwlbaer 
        wlaer(i)=wlbaer(i)
        aerext(i)=qbaer(i)
        aerabs(i)=(1.-wbaer(i))*aerext(i)
        aerasm(i)=gbaer(i)
      enddo

      !write(*,*) 'db usraer 4'
      !write(*,'(a10,/,1p10e11.3)')  'wlaer',(wlaer(i),i=1,nwlbaer)
      !write(*,'(a10,/,1p10e11.3)')  'aerext',(aerext(i),i=1,nwlbaer)
      !write(*,'(a10,/,1p10e11.3)')  'aerabs',(aerabs(i),i=1,nwlbaer)
      !write(*,'(a10,/,1p10e11.3)')  'aerasm',(aerasm(i),i=1,nwlbaer)
      !write(*,*) 'pmaer'
      !do j=1,npmaer/nwlbaer
      !  write(*,'(1p10e11.3)') (pmaer((j-1)*nwlbaer+i),i=1,nwlbaer)
      !enddo
        
      call locate(wlbaer,nwlbaer,.55,j)

      !write(*,*) 'db usraer 5'
      !write(*,*) 'wlbaer ',j,wlbaer(j),wlbaer(j+1)
      f=log(0.55/wlbaer(j))/log(wlbaer(j+1)/wlbaer(j))
      !write(*,*) 'f ',f
      if(0.55.lt.wlbaer(1)) then
        !write(*,*) 'db usraer 5.1',beta
        q55=qbaer(1)*(wlbaer(1)/0.55)**beta
      elseif(0.55.gt.wlbaer(nwlbaer)) then
        !write(*,*) 'db usraer 5.2',beta
        q55=qbaer(nwlbaer)*(wlbaer(nwlbaer)/0.55)**beta
      else
        !write(*,*) 'db usraer 5.3',qbaer(j),f
        q55=qbaer(j)*(qbaer(j+1)/qbaer(j))**f
      endif

      maerw=nwlbaer

c if defined, save aerosol phase function for use in phaerw

      !write(*,*) 'db usraer 6'
      imoma=0
      if(npmaer.gt.0) then
        call phaers(nqbaer,wlbaer,pmaer)
        imoma=1
      endif

      !write(*,*) 'db usraer 7'
      return

      end