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

      subroutine beamsct(nphi,phi,numu,umu,phi0,amu0,nz,dtaug,wreal, 1,1
     &     mpa,sangle,sphase,flxin,dtaur,iout,nbot,ntop,uur)

c     PURPOSE: Add atmospheric single scattering contribution
c     to the radiance array uur.  This routine is only called
c     when an aerosol phase function is specified. The beam contribution
c     is given by equation 6.28 and 6.29 of Liou.  Ignoring boundary
c     contributions, the upwelling radiance is
c
c            nz         /t_j+1                          dt
c    u(i)=sum    w  S P |       exp(-(t-tau)/mu -t/mu0) --
c           j=i   j     /t_j                            mu
c
c           where S is the extraterrestrial solar flux (w/m2/um) (i.e.,
c           what Liou calls pi*Fo), w_j is the single scattering
c           albedo at layer j, P is the phase function, and t_j is the
c           optical depth at boundary j.  The expression for
c           the downwelling radiance is
c
c           i-1         /t_j+1                          dt
c    u(i)=sum    w  S P |       exp(-(tau-t)/mu -t/mu0) --
c           j=1   j     /t_j                            mu
c
c           In these two equations mu and mu0 are the absolute values
c           of the cosine of the viewing and zenith angle, respectively,
c           and the optical depths are computed from input parameter dtaug.
c           u and tau are defined at the boundary between layers, while
c           w, and dtau are defined at the midpoint of each layer.
c
c
c            layer              boundary
c            index              index
c                  ------------   1
c              1
c                  ------------   2
c              2
c                  ------------   3
c              3
c                  ------------   4
c
c
c     on input   uur=Im               multiple scattering (from disort)
c     on output  uur=Im+Is            multiple + single scattering.
c
c

c INPUT:
c  nphi     number of azimuth angles
c  phi      azimuth angles
c  numu     number of zenith angles
c  umu      cosine of zenith angles
c  phi0     zero point of azimuth
c  amu0     solar zenith angle
c  nz       number of levels
c  dtaug    optical depth increments within layer due to aerosol, rayleigh,
c           and gas. Effects of cloud optical depth have not yet been
c           implemented
c  wreal    single scattering albedo
c  mpa      number of entries in phase function array
c  sangle   scattering angles corresponding to entries in phase function array
c  sphase   phase function array.  sphase should be normalized so that
c           integral( sphase, d omega ) = 4pi (same as Liou equation 3.63)
c  flxin    extraterrestrial solar flux density (w/m2/um) 
c  dtaur    optical depth increment within layer due to rayleigh scattering
c  uur      radiance array at computational polar angles

c
c OUTPUT:
c  uur      radiance array with atmospheric single scattering contribution

c discussion:
c
c     may want to change procedure to include spherical earth correction
      parameter (nstrms=40)
      parameter (mxly=50)
      parameter (maxulv=mxly+1)
      parameter (re=6371.2)

      integer nphi,numu,nz,mpa
      real phi(*),umu(*),dtaug(*),dtaur(*),wreal(*),sangle(*),sphase(*)
      real uur(nstrms,maxulv,nstrms)

      real tau(maxulv),expb(maxulv)
      logical debug
      data debug/.false./
      data pi/0./

ccc
c
c spherical earth correction for cosine of zenith angle
c
c      amuz(zzz)=sqrt(1.-(1.-amu0**2)*(re/(re+zzz))**2)
ccc
      !write(*,*) 'db 1'
      if(pi.eq.0.) pi=acos(-1.)

      dtor=pi/180.

      src=flxin/(4.*pi)
      rnorm=pi/(flxin*amu0)  ! conversion to planetary albedo
      sz=amu0
      sx=sqrt(1.-sz**2)

      ! to save time limit calculation to only those
      ! layers destined for output

      nzi=1
      if(iout.eq.5) then
        nz1=ntop
        nz2=ntop
      elseif(iout.eq.6) then
        nz1=nbot
        nz2=nbot
      elseif(iout.eq.20) then
        nz1=ntop
        nz2=ntop
      elseif(iout.eq.21) then 
        nz1=nbot
        nz2=nbot
      elseif(iout.eq.22) then
        nz1=1
        nz2=nz+1
      elseif(iout.eq.23) then
        nz1=ntop
        nz2=nbot
        nzi=nbot-ntop-1
      else
        write(*,*) 'Error -- beamsct '
        write(*,*) '  invalid value of iout: ', iout
        stop 
      endif
        

      ! find total optical depth and extinction factor for direct beam

      !write(*,*) 'db 2'
      tt=0.
      tau(1)=0.
      expb(1)=1.
      do iz=2,nz+1
        tt=tt+dtaug(iz-1)
        tau(iz)=tt
        expb(iz)=exp(-tt/amu0)
      enddo

      phnorm=2.

      ! loop over all angles

      do izen=1,numu
        vz=-umu(izen)  ! in this routine negative mu means up
        vza=abs(vz)
        do iphi=1,nphi

          ! find scattering angle and phase function value

          vx=sqrt(1.-vz**2)*cos((phi(iphi)-phi0)*dtor)
          smu=max(-1.,min(sz*vz+sx*vx,1.))
          sang=acos(smu)/dtor
          call locate(sangle,mpa,sang,j)
          wt=(sang-sangle(j))/(sangle(j+1)-sangle(j))
          sphaer=(1.-wt)*sphase(j)+wt*sphase(j+1)
          sphaer=2.*sphaer/phnorm  ! integral(mu,sphaer)=2
          sphray=.75*(1.+smu**2)   ! integral(mu,sphray)=2

          ! loop over all layers

          do iz=nz1,nz2,nzi

            ! upwelling

            ! write(*,'(a,3f10.5)') 'amu0,vz,phi: ',amu0,vz,phi

            if(vz.lt.0.) then
              thin=0.
              usum=0.
              do izu=iz,nz
                wtray=min(0.,dtaur(izu)/dtaug(izu))
                sph=sphaer*(1.-wtray)+sphray*wtray
                ff=sph*src*wreal(izu)
                amufac=amu0/(vza+amu0)
                exp1=exp(                earg=dtaug(izu)*(1./amu0+1./vza)
                if(earg.lt.0.0001) then
                  efac=earg
                else
                  efac=1.-exp(-earg)
                endif
                usum=usum+ff*amufac*expb(izu)*exp1*efac

                if(debug) then
                  thin=thin+ff*dtaug(izu)/vza
                 
                  if(izu.eq.iz) write(*,'(4a4,11a11)')'iz','izu','phi',
     &                 'zen','tau','e*e','efac','wtray','sph',
     &                 'wreal','u','thin'
                  
                  write(*,'(4i4,1p11e11.3)') iz,izu,iphi,izen,
     &                 tau(izu),expb(izu)*exp1,efac,wtray,sph,
     &                 wreal(izu),usum,thin
                endif

              enddo

            ! downwelling

            elseif(vz.gt.0.) then

              thin=0.
              usum=0.
              do izd=1,iz-1
                wtray=min(0.,dtaur(izu)/dtaug(izd))
                sph=(sphaer*(1.-wtray)+sphray*wtray)
                ff=sph*src*wreal(izd)
                earg=dtaug(izd)*abs(1./amu0-1./vza)
                if(earg.lt.0.0001) then
                  efac=earg
                else
                  efac=1.-exp(-earg)
                endif
                if(vza.gt.amu0) then
                  amufac=amu0/(vza-amu0)
                  exp1=exp(                  usum=usum+ff*amufac*exp1*expb(izd)*efac
                elseif(vza.lt.amu0) then 
                  amufac=amu0/(amu0-vza)
                  exp1=exp(                  usum=usum+ff*amufac*exp1*expb(izd+1)*efac
                else
                  usum=usum+ff*expb(iz)*dtaug(izd)/amu0
                endif
                if(debug) then
                  thin=thin+ff*dtaug(izd)/vza

                  if(izd.eq.1) write(*,'(4a4,11a11)')'iz0','izd',
     &                 'phi','zen','tau','e*e','efac','wtray','sph',
     &                 'wreal','u','thin'
                  write(*,'(4i4,1p11e11.3)') iz,izd,iphi,izen,
     &                 tau(izd),expb(izd)*exp1,efac,wtray,sph,
     &                 wreal(izd),usum,thin
                endif
              enddo

            else

              write(*,*) 'Error -- beamsct'
              write(*,*) ' viewing zenith angle = 90'
              write(*,*) 'iz,vz  ',iz,vz
              stop

            endif

            uur(izen,iz,iphi)=uur(izen,iz,iphi)+usum

          enddo
        enddo
      enddo

      return
      end