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


      subroutine tauaero(wl,nz,z,v55,taer55,naerz,taerst,laer,jaer, 1,5
     &                   waer,gaer,taua,tauab,imoma)
c
c purpose:  
c     compute the optical depth and scattering parameters at all
c     atmospheric levels of both boundary layer and stratospheric
c     aerosols
c
c input:
c   wl      wavelength
c   nz      number of atmospheric levels
c   z       altitude (km)
c   wh      water vapor density (g/m3)
c   v55     visibility at 550nm 
c   taer55  optical depth at 550nm  
c
c           (v55 and taer55 are alternate ways of setting the aerosol 
c           optical depth. Note that v55 cant be used to set aerosol
c           optical depth when the aerosol density in the lowest layer
c           is set to zero through zbaer and dbaer)
c
c naer number of stratospheric aerosol scattering layers
c laer index array of stratospheric aerosols (see subroutine zlayer)
c taerst optical depth array of stratospheric aerosol scattering layers
c jaer array of stratospheric aerosol types
c
c             0 = no stratospheric aerosols
c             1 = background        (between 10 and 30 km altitude)
c             2 = aged volcanic     (between 10 and 30 km)
c             3 = fresh volcanic    (between 10 and 30 km)
c             4 = meteor dust       (above 30 km)
c
c output:
c   waer    single scatter albedo due to aerosols at each atmospheric layer
c   gaer    depends on the value of imoma
c           imoma=0 asymmetry factor due to aerosols at each atmospheric layer
c           imoma=1 optical depth of boundary layer aerosol. this output
c                   is used to give the proper weight to the boundary
c                   layer phase function in routine depthscl

c   taua   aerosol optical depth increment at each atmospheric layer
c   tauab  boundary layer aerosol optical depth increment
c
c NOTE:     visibility is defined as 3.912/extinction_coefficient(km-1)
c           for a horizontal path 
c           (exp(-3.912) = 0.02 = visible contrast threshold)
c
c exp(-3.912)=.02 
c
c revisions:
c  971028: corrected expression for net asymmetry factor including both 
c          stratospheric and boundary layer aerosols (suggested by Xiang Li)
c
      parameter (visfac=3.912)
      parameter (ndbug=1)
      dimension z(*),taerst(*),laer(*),jaer(*),waer(*),gaer(*)
      dimension taua(*),tauab(*)

      data init,sigma/1, 0./
      save init, sigma 
      ! sigma is a wavelength independent cross-section for
      ! scattering by boundary layer aerosols (bla)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

c integrate vertical profile of boundary layer aerosol density

      call aervint(v55,nz,z,taua) ! taua = n(i) dz 

c initialize value of sigma on first call of tauaero
c the value of sigma (which independent of wavelenth)
c is saved for later calls to tauaero
      
      if(init.eq.1) then
        init=0
        sigma=0.
        call aerbwi(0.55,ext55,w55,g55)
        if(ext55.gt.0.) then 
          if(taer55.ge.0) then
            vint=total(taua,1,nz)
            if(vint.ne.0) sigma=taer55/(ext55*vint)
          else
            call aeroden(            sigma=visfac/(ext55*v55*aden0) 
          endif
        endif
      endif
        
c compute scattering parameters (extinction efficiency , single 
c scattering albedo and asymmetry factor) for boundary layer 
c aerosols at a given wavelength. Note that extinction efficiency
c is normalized to 1 at 550nm.

      call aerbwi(wl,extinc,wa,ga)

      ! write(*,'(a5,9a11)')
      !&     'init','wl','sigma','extinc','taer55','ext55','vint'
      ! write(*,'(i5,1p9e11.3)')init,wl,sigma,extinc,taer55,ext55,vint

      do 30 i=1,nz
        taua(i)=extinc*sigma*taua(i)
        tauab(i)=taua(i)
        waer(i)=wa
        gaer(i)=ga
        if(imoma.eq.1) gaer(i)=0.
 30   continue

c
c add in stratospheric aerosols at specified scattering layers
c

      do 40 i=1,naerz
        if(jaer(i).ne.0 .and. taerst(i).gt.0.) then
          nl=laer(i)
          ja=jaer(i)

c find scattering parameters for stratospheric
c aerosols at a given wavelength

          call aerswi(ja,wl,extinc,wa,ga)
          dt=taerst(i)*extinc
          if(sigma.ne.0) then
            ga=(dt*ga*wa+taua(nl)*gaer(nl)*waer(nl))/
     &           (dt*wa+taua(nl)*waer(nl))
            wa=(dt*wa+taua(nl)*waer(nl))/(dt+taua(nl))
            dt=dt+taua(nl)
          endif
          waer(nl)=wa
          gaer(nl)=ga
          taua(nl)=dt
        endif
 40   continue

c      write(*,'(/,1p,a,/(10e11.3))') 'waer',(waer(j),j=1,33)
c      write(*,'(/,1p,a,/(10e11.3))') 'gaer',(gaer(j),j=1,33)
c      write(*,'(/,1p,a,/(10e11.3))') 'taua',(taua(j),j=1,33)

      return

      end