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


      subroutine taugas(wl,uu,amu0,nz,z,dtauc,dtaul,lun) 2,27
c
c     calculates optical depth due to atmospheric gases and aerosols
c
c  input:  wl        wavelength in microns
c          uu        path integral uu(mxq,nxlayr) computed by ABSINT
c          amu0      cosine of solar zenith angle (at surface)
c          nz        number of atmospheric layers
c                    nz=1  => integration to bottom level only
c          z         heights of atmospheric layers
c          lun       if non-zero print optical depth profiles to LUN
c
c  output: 
c          dtauc(k)  the optical depth of a single layer, k, due to 
c                    continuum processes (follows Beer's law).
c
c          dtaul(k)  the negative log of the mean transmission through
c                    a single layer k, due to line absorption (no
c                    continuum processes).  BE AWARE OF THIS: the
c                    negative log of the transmission through a given
c                    layer depends on the column density of absorbers 
c                    above that layer.
c
c          NOTE: the order of the depth layers are reversed in this routine,
c                uu(k,1) will be used to compute the optical depth dtauc(nz)
c                    
c
c                    continuum and line absorption due to 
c
c                    1. uniformly mixed atmospheric gases (co2,n2,o2...), 
c                    2. h2o
c                    3. ozone
c                    4. trace gases.
c
c          gasabs    gas absorption diagnostics at surface, gasabs(9),
c                    due to both line and continuum processes
c                   1    2    3   4    5   6     7      8      9   
c                  h2o  co2  o3  n2o  co  ch4  o2+n2  trace  total 
c     
cc***********************************************************************
c     
      parameter (bigexp=87.)
      parameter (mxq=63)
      parameter (grav=9.80665)                ! gravitation acc. (m s-2)
      parameter (airrho=1.273)                ! air density at stp (kg m-3)
      parameter (pzero=1013.25)               ! standard pressure (mb)
      parameter (amfac=10.*grav*airrho/pzero) ! convert scale height to airmass

      dimension uu(mxq,nz),w(mxq)
      dimension dtauc(*),dtaul(*)

      parameter (mxly=50)
      parameter (re=6371.2)

      dimension tauc(mxly),taul(mxly)

      dimension z(*)

      parameter (nta=9)

      common /gasblk/ airmass(mxly),gasabs(nta,mxly)

      common /h2o/    cph2o(3515)
      common /o3/     cpo3 ( 447)
      common /ufmix1/ cpco2(1219)
      common /ufmix2/ cpco ( 173),cpch4( 493),cpn2o( 704),cpo2 ( 382)
      common /traceg/ cpnh3( 431), cpno(  62),cpno2( 142),cpso2( 226)
      common /wnlohi/
     l     iwlh2o(15),iwlo3 ( 6),iwlco2(11),iwlco ( 4),iwlch4( 5),
     l     iwln2o(12),iwlo2 ( 7),iwlnh3( 3),iwlno ( 2),iwlno2( 4),
     l     iwlso2( 5),
     h     iwhh2o(15),iwho3 ( 6),iwhco2(11),iwhco ( 4),iwhch4( 5),
     h     iwhn2o(12),iwho2 ( 7),iwhnh3( 3),iwhno ( 2),iwhno2( 4),
     h     iwhso2( 5)
      common /aabbcc/ aa(11),bb(11),cc(11),a(11),cps(11),ibnd(11)

      dimension tau(11)
      data indh2o,indo3,indco2,indco,indch4/5*1./
      data indn2o,indo2,indnh3,indno,indno2,indso2/6*1./
      
c***********************************************************************
c non-spherical earth correction to the cosine of the solar zenith angle
c Draw a line extending from a point on the earths surface out to space
c if amu0 is the solar zenith angle at the surface then amuz(zzz) is the
c angle between the line and the local zenith direction at altitude zzz. 

c statement function:

      amuz(zzz)=sqrt(1.-(1.-amu0**2)*(re/(re+zzz))**2)

      do 10 i=1,nz
        tauc(i)=0.
        taul(i)=0.
 10   continue

c convert from wavelength in microns to wavenumber in cm-1

      iv=5*(int(10000.0/wl)/5)
      v=float(iv)
      v=10000./wl

c water continuum (v < 10000 cm-1)
c slf296 loads self-broadened water vapor continuum at 296k
c slf260 loads self-broadened water vapor continuum at 296k
c frn296 loads foreign-broadened water vapor continuum at 296k

      call slf296(v,sh2ot0)
      call slf260(v,sh2ot1)
      call frn296(v,fh2o)
      t0=296.
      t1=260.
      if(sh2ot0.gt.0.) then
        alpha2=200.**2
        xh2o=(1.-0.2333*(alpha2/((v-1050.)**2+alpha2)) )
        sh2ot0=sh2ot0*xh2o
        sh2ot1=sh2ot1*xh2o
      endif

c protect against exponential underflow at high frequency

      vtemp=v/0.6952
      if(vtemp/t1.le.bigexp) then
        xd=exp(-v/(t0*0.6952))
        radfn0=v*(1.-xd)/(1.+xd)
        xd=exp(-v/(t1*0.6952))
        radfn1=v*(1.-xd)/(1.+xd)
      else
        radfn0=v
        radfn1=v
      endif

      wfac=1.e-20
      ya=exp(-log(1.025*3.159e-8)+(2.75e-4)*v)
      yb=exp(-log(8.97e-6)+(1.300e-3)*v)
      fdg=1./(ya+yb)

c  n2 continuum absorption coefficient (2080 < v < 2740)

      call c4dta(abn2,v)

c hno3 absorption calculation (850 < v < 1735)

      call hno3(v,abno3)

c hertda computes hertzburg uv o2 absorption ( v > 36000 )  
c o2cont computes o2 continuum ( 1395 < v < 1760 )

      call hertda(abo2,v)
      call o2cont(v,sigo20,sigo2a,sigo2b)

c o2 * o2 collision-radiative process 

      call o4cont(wl,sigo4)

c diffuse ozone

c o3uv computes uv ozone for 40800-54054cm-1 (185-245nm)

c o3hht# ozone hartley band for 24370-40800 cm-1 (245-410 nm)
c   temperature dependent coefficient
c   o3hht0 (contant term), o3hht1 (linear term), o3hht2 (quadratic term)

c c8dta ozone chappius band for 13000-24200 cm-1 (413-769nm)


      doz1=0.
      doz2=0.
      doz3=0.
      if(v.gt.40800) then 
        call o3uv(v,c0)
        doz1=.269*c0
      elseif(v.gt.24370) then        
        call o3hht0(v,c0)
        call o3hht1(v,ct1)
        call o3hht2(v,ct2)
        doz1=.269*c0
        doz2=c0*ct1
        doz3=c0*ct2
      elseif(v.ge.13000. .and. v.le.24200) then
        call c8dta(abo3a,v)
        doz1=abo3a
      endif

c rayleigh scattering (included here for print out only)

      call c6dta(rayla,v)

c cxdta locates coefficient for double exponential

      call cxdta(cps(1),v,iwlh2o,iwhh2o,cph2o,indh2o)
      call cxdta(cps(2),v,iwlco2,iwhco2,cpco2,indco2)
      call cxdta(cps(3),v,iwlo3, iwho3, cpo3, indo3 )
      call cxdta(cps(4),v,iwln2o,iwhn2o,cpn2o,indn2o)
      call cxdta(cps(5),v,iwlco, iwhco, cpco, indco )
      call cxdta(cps(6),v,iwlch4,iwhch4,cpch4,indch4)
      call cxdta(cps(7),v,iwlo2, iwho2, cpo2, indo2 )
      call cxdta(cps(8),v,iwlno, iwhno, cpno, indno )
      call cxdta(cps(9),v,iwlso2,iwhso2,cpso2,indso2)
      call cxdta(cps(10),v,iwlno2,iwhno2,cpno2,indno2)
      call cxdta(cps(11),v,iwlnh3,iwhnh3,cpnh3,indnh3)

c  abcdta moves double exponential coefficients to new arrays

      call abcdta(iv)

c schrun computes uv o2 schumann-runge band model parameters

      if(v.gt.49600) call schrun(v,cps(7))

      zim=z(nz)
      do 40 i=nz,1,-1
        im=nz-i+1
        zi=z(i)
        zbar=0.5*(zi+zim)
        zim=zi
        do 20 k=1,mxq
          if (i .eq. nz) then
            w(k)=uu(k,i)/amuz(zi)
          else
            w(k)=w(k)+(uu(k,i)-uu(k,i+1))/amuz(zbar)
          endif
 20     continue

c        if(i.eq.nz) write(*,'(2x,5a11)') 'zi','zbar','w(6)','uu','amu'
c        write(*,'(i2,1p5e11.3)') i,zi,zbar,w(6),uu(6,i),amuz(zbar)

        tcunif=+sigo4*w(3)
     &         +abn2*w(4)
     &         +sigo20*(     &         +sigo2a*(     &         +sigo2b*w(2))
     &         +abo2*w(58)
        tch2o=sh2ot0*radfn0*(wfac*w(5))
     &        +((sh2ot1*radfn1)-(sh2ot0*radfn0))*(wfac*w(9))
     &        +(fh2o+fdg)*radfn0*(wfac*w(10))
        tco3=doz1*w(8)+doz2*w(59)+doz3*w(60)
        tctrc=abno3*w(11)

c  compute transmitance from double exponential band model for
c       1    2    3    4    5    6    7    8    9   10   11
c      h2o, co2,  o3, n2o,  co, ch4,  o2,  no, so2, no2, nh3

        tauc(im)=tcunif+tch2o+tco3+tctrc

c        write(*,'(a,/,1p(10i11))') 'ibnd:',(ibnd(k),k=1,11)
c        write(*,'(a,/,1p(10e11.3))') 'cps:',(cps(k),k=1,11)
c        write(*,'(a,/,1p(10e11.3))') 'w:',(w(ibnd(k)),k=1,11)

        do 30 k=1,11
          ib=ibnd(k)
          cp=cps(k)
          if(cp.gt.-20. .and. w(ib).gt.1.e-20) then
            awl=a(k)*(cp+log10(            awl=min(awl,20.)              
            tau(k)=10.**awl
          else
            tau(k)=0.
          endif
          taul(im)=taul(im)+tau(k)
 30     continue

c  gas absorption diagnostics in gasabs(9)
c    1    2    3   4    5   6     7      8      9   
c   h2o  co2  o3  n2o  co  ch4  o2+n2  trace  total 

        gasabs(1,i)=tch2o+tau(1)
        gasabs(2,i)=tau(2)
        gasabs(3,i)=tau(3)+tco3
        gasabs(4,i)=tau(4)            
        gasabs(5,i)=tau(5)         
        gasabs(6,i)=tau(6)
        gasabs(7,i)=tau(7)+tcunif
        gasabs(8,i)=tau(8)+tau(9)+tau(10)+tau(11)
        gasabs(9,i)=tauc(im)+taul(im)

        airmass(i)=amfac*w(6)

        if(lun.gt.0) then

          tauray=rayla*w(6)          
          if(i.eq.nz) write(lun,'(3x,11a11)') 
     &             'air mass','h2o','co2','o3','n2o','co',
     &             'ch4','o2+n2','trace','total','molec'

          airmass(i)=0.

          write(lun,'(i3,1p11e11.3)') i,airmass(i),
     &         tch2o+tau(1),  tau(2),  tau(3)+tco3,
     &         tau(4),  tau(5),  tau(6),  tau(7)+tcunif,
     &         tau(8)+tau(9)+tau(10)+tau(11),  tauc(im)+taul(im),
     &         tauray

        endif
 40   continue
      
      dtauc(1)=tauc(1)
      dtaul(1)=taul(1)
 9876 format(/a,1p,4(/10e11.3))
c      write(*,*) 'taugasx    amu0=',amu0
c      write(*,9876) 'tauc+taul',(tauc(i)+taul(i),i=1,nz)

      do 50 i=2,nz
        dtauc(i)=tauc(i)-tauc(i-1)
        dtaul(i)=taul(i)-taul(i-1)
 50   continue
c      write(*,9876) 'dtauc+dtaul',(dtauc(i)+dtaul(i),i=1,nz)
      return
      end