c file:                  taugas.f
c
c external routines:     kdistr,absint,taugas,rayleigh,gasabs
c
c internal routines:     slf296,slf260,frn296,sint,c4dta,c6dta,hno3,
c                        hertda,o2cont,o2int,o3hht0,o3hht1,o3hht2,
c                        o3int,o3uv,c8dta,cxdta,abcdta,schrun,volmix
c
c internal common        aabbcc,abc,atm,c4c8,esfblk,fh2o,h2o,o2c,o3,
c                        o3hh0,o3hh1,o3hh2,o3uvf,s260,sh2o,shur,traceg,
c                        ufmix1,ufmix2,wnlohi,absgas
c
c=======================================================================

      subroutine kdistr(nz,mxly,mxq,uu,dtauk,tauk,wtk,twgp) 1

c   PURPOSE: compute 3 component k-distribution optical depths
c            and weights.  The common block parameters in /aabbcc/
c            are set up by subroutine TAUGAS, which should be called
c            before KDISTR at each new wavelength.
c
c   input parameters:
c   ----------------
c     nz       number of active layers
c     mxly     dimensioned size of first subscript of dtauk,tauk,wtk,twgp
c     mxq      dimensioned size of first subscript of uu
c     uu       absorption integrals from ABSINT 
c
c   output parameters:
c   -----------------
c   dtauk(nz,3)   molecular optical thickness of layer n, term k 
c   tauk(nz,3)    cumulative optical depth upto layer n, term k
c   wtk(nz,3)     layer by layer weighting for k-distribution fit
c     
c
c   internal working array (dimensioned in calling routine)
c   -------------------------------------------------------
c   twgp(nz,3)  - sum of optical depth * probability by molecule
c
c***********************************************************************
c
c   new varables          11 molecules  by joseph h pierluissi
c   dpwj  probability for each molecule  fit double exponential
c   gkwj  band dependent  scaling of densitys to get k  amount
c   cps is the stored values of Pierluissi band model coefficients
c   cp1s = 10**cps
c   ibnd maps bands to molecules
c   wtk  is the  effective probability by layer
c   dtaum is defined as the sum of the optical depths by molecule
c
c************************************************************************

      parameter (nk=3)

      common /aabbcc/ aa(11),bb(11),cc(11),a(11),cps(11),ibnd(11)

      dimension dtauk(mxly,nk),tauk(mxly,nk),wtk(mxly,nk),twgp(mxly,nk)
      dimension fac(nk),gkwj(nk,11),dpwj(nk,11),cp1s(11)
      dimension taudb(11)
      dimension uu(mxq,mxly)

      data fac /1.0,0.09,0.015/

c************************************************************************

      dpc = 1./3.
      
      do 20 m = 1,11
        cp1s(m)= 10.**cps(m)
        taudb(m)=0.
 20   continue

      do  40   k = 1,nk
        do  30  mol = 1,11
          iw = ibnd(mol)
          gkwj(k,mol) = 0.
          dpwj(k,mol) = 0.
          if(iw. gt. 0) then
            gkwj(k,mol) = fac(k) * cc(mol)
            if(k .eq. 1) dpwj(k,mol)=aa(mol)
            if(k .eq. 2) dpwj(k,mol)=bb(mol)
            if(k .eq. 3) dpwj(k,mol)=1.-aa(mol)-bb(mol)
          endif
 30     continue
 40   continue
      
c   evaluate the weighted k distribution quantities for
c   water vapor, ozone and uniformly mixed gases

      taukk=0.
      do 80 k=1, nk
        do 70 n=1,nz
          tauk(n,k)=0.
          dtauk(n,k)= 0.
          nm=nz-n+1

          twgp(n,k)=0.
          do 60 mol=1,11
            ib = ibnd(mol)
            if(ib.lt.0) go to 60
            if(nm.eq.nz) then
              duu=uu(ib,nm)
            else
              duu=uu(ib,nm)-uu(ib,nm+1)
            endif
            wpth = duu*gkwj(k,mol)
            dtauk(n,k)=dtauk(n,k)+wpth*cp1s(mol)
            twgp(n,k)=twgp(n,k)+wpth*cp1s(mol)*dpwj(k,mol)
 60       continue
          
          wtk(n,k) = dpc
          if(dtauk(n,k).ne.0) wtk(n,k)=twgp(n,k)/dtauk(n,k)
          
c     effective probability by layer, wtk is based on 
c     molecular probability weighted by optical depth
          
 70     continue
 80   continue
      
      do 100 n=1,nz
        smwtk = wtk(n,1) + wtk(n,2) + wtk(n,3)
        if(smwtk .eq. .0) write(8, *) 'smwtk == 0'
        do 90 k=1, nk
          wtk(n,k) = wtk(n,k)/smwtk
          if(n.eq.1) then
            tauk(n,k) = dtauk(n,k)
          else
            tauk(n,k) = dtauk(n,k)+tauk(n-1,k)
          endif
 90     continue
 100  continue
      
      return
      end