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

      subroutine depthscl(kdist,nz,k,npass,nstr,wl,amu0,z,dtls, 1,3
     &     dtcs,dtcv,dtlv,dtaur,dtaua,dtauab,taucld,dtauk,wtk,wcld,waer,
     &     gcld,gaer,dtaug,wreal,pmom,wt,imom,imoma,lun)
c
c routine:   depthscl
c
c purpose:   compute optical depths and other layer-by-layer info
c
c input:
c  nz       number of atmospheric layers
c  k        term in k-distribution expansion
c  npass    total number of k-distribution terms in expansion
c  wl       wavelength (microns) 
c  amu0     cosine of solar zenith
c  z        height of layer (km)
c  dtcs     optical depth increment along slant path due to gaseous continua
c  dtls     o.d.i. along slant path due to gaseous lines (non beers)
c  dtcv     o.d.i. along vertical due to gaseous continua
c  dtlv     o.d.i. along vertical due to gaseous lines (non beers)
c  dtaua    o.d.i. along vertical due to all aerosols
c  dtauab   o.d.i. along vertical due to boundary layer aerosols
c  taucld   o.d.i. along vertical due to clouds
c  dtauk    o.d.i. of all three k-distibution terms
c  wtk      probabilities of k-distribution terms
c  wcld     single scatter albedo due to clouds
c  waer     single scatter albedo due to aerosols
c  gcld     asymmetry factor for clouds
c  gaer     asymmetry factor for aerosols
c  pmom     moments of the the CLOUD phase function (set in taucloud)
c  imom     0 use Henyey-Greenstein phase function for cloud scattering
c           1 use cloud scattering phase function read from pmom.dat
c  imoma    0 use Henyey-Greenstein phase function for b.l. aerosol scattering
c           1 use aerosol scattering phase function read from INPUT
c  lun      logical unit for diagnostic output, no output for lun=0

c
c   output:
c  dtaug    optical depth increment along vertical
c  wreal    effective single scatter albedo
c  pmom     moments of TOTAL phase function
c  wt       probabilities for the k'th k-distribution pass 


      parameter (nstrms=40)
      parameter (mxly=50)
      parameter (mxq=63)
      parameter (re=6371.2, pi=3.1415926536)
      parameter (prec=1.0d-7)

      dimension z(*),dtls(*),dtcs(*),dtcv(*),dtlv(*),
     &     dtaur(*),dtaua(*),dtauab(*),taucld(*),dtauk(mxly,3),
     &     wtk(mxly,3),wcld(*),waer(*),
     &     gcld(*),gaer(*),wreal(*),pmom(0:nstrms,*)

      dimension pmcld(0:nstrms),pmaer(0:nstrms),
     &          greal(mxly),dtaug(mxly),taug(mxly)
c
c
      tautot=0.
      tauray=0.
      tauaer=0.
      tcloud=0.
      tauls=0.      
      taulv=0.
      tauc=0.
      tauk=0.
      tauscat=0.
      tgline=0.
      tgcont=0.
      taug1=0.
      taug2=0.
      taug3=0.       

      if(npass.gt.1) then
        wnorm=0.
        w1=0.
        w2=0.
        w3=0.
        do 10 j=1,nz
          w1=wtk(j,1)*dtlv(j)+w1
          w2=wtk(j,2)*dtlv(j)+w2
          w3=wtk(j,3)*dtlv(j)+w3
 10     continue
        wnorm=w1+w2+w3
        if(wnorm.eq.0) then
          w1=1.
          w2=0.
          w3=0.
        else
          w1=w1/wnorm
          w2=w2/wnorm
          w3=w3/wnorm
        endif
        if(k.eq.1) wt=w1
        if(k.eq.2) wt=w2
        if(k.eq.3) wt=w3
      else
        wt=1.0
      endif

      if(lun.ne.0) then
        if(npass.gt.1.and.k.eq.1) then
          write(lun,'(a/1p(10e11.3))') 'dtauk(j,1)',(dtauk(j,1),j=1,nz)
          write(lun,'(a/1p(10e11.3))') 'dtauk(j,2)',(dtauk(j,2),j=1,nz)
          write(lun,'(a/1p(10e11.3))') 'dtauk(j,3)',(dtauk(j,3),j=1,nz)
          write(lun,'(a/1p(10e11.3))') 'wtk(j,1)',(wtk(j,1),j=1,nz)
          write(lun,'(a/1p(10e11.3))') 'wtk(j,2)',(wtk(j,2),j=1,nz)
          write(lun,'(a/1p(10e11.3))') 'wtk(j,3)',(wtk(j,3),j=1,nz)
          write(lun,'(/a,1p3e11.3/)') 'w1,w2,w3: ',w1,w2,w3
        endif
        write(lun,'(a,i3,a,f10.4)') 'k =',k,'    wl =',wl
        write(lun,'(2x,a8,8a11)')
     &       'z','dtau','dtaua','dtaur','dtcld','dtauc','dtaul',
     &       'w','g'
      endif

      do 30 j=1,nz
        dtauc=dtaur(j)+dtaua(j)+taucld(j)+dtcv(j)
        tauc=tauc+dtauc
        dtsct=dtaur(j)+wcld(j)*taucld(j)+waer(j)*dtaua(j)
        tauscat=tauscat+dtsct
        tauls=tauls+dtls(j)
        taulv=taulv+dtlv(j)
        ultrans=exp(-tauls)

c NOTE: gas continua absorption is grouped with other continuous
c       opacities,  this assumes dtcs=dtcv/amu


c normalize direct beam to slant path transmission by
c adjusting k-distributions optical depths.  
c kdist=1 no renormalization
c kdist=2 renormalize (clear sky situation)
c kdist=3 ramp correction factor to one as tauscat goes > 1

        if(npass.gt.1) then 
          taugx1=taug1+dtauk(j,1)
          taugx2=taug2+dtauk(j,2)
          taugx3=taug3+dtauk(j,3)

          if(tauls.lt.12. .and. kdist.ge.2) then
            call taucor(w1,w2,w3,taugx1,taugx2,taugx3,amu0,tauls,cf)
            if(kdist.eq.3.and.tauscat.gt.1) then
              ramp=exp(-tauscat)/exp(-1.)
              cf=cf*ramp+(1.-ramp)
            endif
            taugx1=taugx1*cf
            taugx2=taugx2*cf
            taugx3=taugx3*cf
          endif
          
          if(k.eq.1) dtgasl=(taugx1-taug1)
          if(k.eq.2) dtgasl=(taugx2-taug2)
          if(k.eq.3) dtgasl=(taugx3-taug3)
          dtgasl=max(dtgasl,0.)

          taug1=taugx1
          taug2=taugx2
          taug3=taugx3
          transk=w1*exp(-taug1/amu0)
     &          +w2*exp(-taug2/amu0)
     &          +w3*exp(-taug3/amu0)
        else
          if(taulv.gt..001) then
            afac=amu0*tauls/taulv
          else 
            afac=1.
          endif
          if(tauscat.gt.1.) then
            ramp=exp(-tauscat)/exp(-1.)
            afac=afac*ramp+(1.-ramp)
          endif
          dtgasl=afac*dtlv(j)
          transk=ultrans
        endif
        
        dtaug(j)=dtgasl+dtauc
        tautot=tautot+dtgasl+dtauc
        taug(j)=tautot
        wreal(j)=0.
        if(dtaug(j).ne.0.) wreal(j)=dtsct/dtaug(j)

c NOTE: WREAL is limited to values below (1.0 - prec) to prevent a 
c divide by zero error in DISORT when the gas absorption is very low.

        wreal(j)=min(wreal(j),1.-prec)

        if(dtsct.ne.0.) greal(j)=
     &       (gcld(j)*wcld(j)*taucld(j)+gaer(j)*waer(j)*dtaua(j))/dtsct
         
        call getmom(3,gaer(j),nstrms,pmaer)

        if(imoma.eq.1) call phaerw(wl,dtaua(j),dtauab(j),pmaer)
        if(imom.ne.1) call getmom(3,gcld(j),nstr,pmcld)

        do 20 i=0,nstrms
          if(imom.eq.1) pmcld(i)=pmom(i,j)
          pmom(i,j)=0.
          if(dtsct.ne.0.) pmom(i,j)=(pmcld(i)*wcld(j)*taucld(j)+
     &                               pmaer(i)*waer(j)*dtaua(j))/dtsct
 20     continue

        pmom(0,j)=1.0

c add in contribution from Rayleigh scattering

        if(dtsct.ne.0.) pmom(2,j)=pmom(2,j)+0.1*dtaur(j)/dtsct

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

c print diagnostics

        if(lun.ne.0) then 
           tauaer=tauaer+dtaua(j)
           tauray=tauray+dtaur(j)
           tcloud=tcloud+taucld(j)
           tgcont=tgcont+dtcv(j)
           tgline=tgline+dtgasl
           
           write(lun,'(i2,f8.2,1p8e11.3)')
     &         j,z(nz-j+1),dtaug(j),dtaua(j),dtaur(j),taucld(j),
     &          dtcv(j),dtgasl,wreal(j),greal(j)
        endif
        
 30   continue

      if(lun.ne.0.) then 
        write(lun,'(10x,8("  ---------"))') 
        write(lun,'(2x,a8,1p6e11.3)') "totals: ",
     &       tauaer+tauray+tcloud+tgcont+tgline,
     &       tauaer,tauray,tcloud,tgcont,tgline
        
        if(imom.ne.0.or.imoma.ne.0) then
          write(lun,*)
          write(lun,'(a,/,1p(10e11.3))') 'pmom:',(pmom(i,nz),i=0,nstrms)
          write(lun,'(a,/,1p(10e11.3))') 'pmaer:',(pmaer(i),i=0,nstrms)
        endif
      endif
c  compute optical depth correction to user angle slant paths

      return
      end