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


      subroutine absint(uu,nz,z,p,t,wh,wo,lun) 1,1
c
c     calculates absorber integrals from top to bottom of atmosphere
c
c    input:  nz         number of atmospheric layers
c            z          layer altitudes
c            p          pressure in mb
c            t          temperature in kelvin 
c            wh         water vapor mass density in g/m3
c            wo         o3 mass density in g/m3
c            lun        if non-zero, write diagnostic print out to LUN
c
c    output: uu         absorption integrals
c
c
c     k      absorption process
c    ----    --------------------------------------
c     1-2    temp dep o2  1395 to 1760
c     3      o2*o2 and o2*n2 continuum (solomon et al 1998) 
c     4      n2 continuum
c     5      h2o self broadened continuum at 296 k
c     6      molecular (rayliegh) scattering
c     7      not used
c     8      ultraviolet ozone
c     9      h2o self continuum 296-260  (tx(9) total trans)
c    10      h2o foriegn continuum
c    11      hno3 (nitric acid)
c    12-16   not used
c    17-30   h2o (line absorption) (double exponential models)
c    31-35   o3  (line absorption) (double exponential models)
c    36-43   co2 (line absorption) (double exponential models)
c    44-45   co  (line absorption) (double exponential models)
c    46      ch4 (line absorption) (double exponential models)
c    47      n2o (line absorption) (double exponential models)
c    50-51   o2  (line absorption) (double exponential models)
c    52-53   nh3 (line absorption) (double exponential models)
c    54      no  (line absorption) (double exponential models)
c    55      no2 (line absorption) (double exponential models)
c    56-57   so2 (line absorption) (double exponential models)
c    58      o2  hertzberg continuum (uv)
c    59-60   temp dep o3
c    61-62   not used
c    63      o2  1395 to 1760
c
c***********************************************************************
      parameter (mxly=50)
      parameter (mxq=63)
      parameter (alosch=2.6868e19)
      parameter (xlosch=alosch*1.e5)
      parameter (pmo=2.6568e-23)
      parameter (conjoe=0.1/alosch)
      parameter (pzero=1013.25, tzero=273.15)
      parameter (con=3.3429e21)
      parameter (rhzero=273.15/296.0)
c
c conjoe = 1.e5(cm/km) * 1.e-6(ppmv) / alosch
c con    = 0.1/molecular_mass_of_h2o 
c          con/xlosch = conversion from g/m3 => density in loschmidts
c pmo    = mass of oxygen atom (g)
c
      dimension z(mxly),p(mxly),t(mxly),wh(mxly),wo(mxly)
      dimension uu(mxq,mxly),dd(mxq,mxly)
c
c**********************************************************************
c
      do 100 i=1,nz
        call volmix(     &  vfn2,vfo2,vfco2,vfch4,vfn2o,vfco,vfno2,vfso2,vfnh3,vfno,vfhno3)
        tt=t(i)
        pp=p(i)
        pss=pp/pzero
        tss=tzero/tt
        f1=(pp/pzero)/(tt/tzero)
        f2=(pp/pzero)*sqrt(tzero/tt)
        wair=alosch*f1

        rhoair=f1
        rhoh2o=con*wh(i)/xlosch     
        rhofrn=rhoair-rhoh2o
        wo2d=conjoe*wair*vfo2*pss
        vfo3=wo(i)/(3*pmo*wair)
c
c                             rhoair - air density in loschmidts
c                             rhoh2o - water vapor density in loschmidts
c                             
c
c vfo3 = o3 volume fraction in parts per million
c
        dd(1,i)=wo2d*tt
        dd(2,i)=wo2d*(tt-220.)**2
        dd(3,i)=f1**2
        dd(4,i)=1.e-6*vfn2*f1*f2
        dd(5,i)=xlosch*rhoh2o**2/rhzero
        dd(6,i)=f1
        dd(8,i)=conjoe*wair*vfo3
        dd(10,i)=xlosch*rhoh2o*rhofrn/rhzero
        dd(11,i)=f1*vfhno3*(1.e-6*1.e5)
        dd(63,i)=wo2d
c
c for h2o
c    wh(i) is in g/m3,  conh2o is in g/cm2/km
c
        conh2o=wh(i)*.1
        dd(17,i)=conh2o*pss**0.9810*tss**( 0.3324) !     0-  345 cm-1 
        dd(18,i)=conh2o*pss**1.1406*tss**(-2.6343) !   350- 1000 cm-1 
        dd(19,i)=conh2o*pss**0.9834*tss**(-2.5294) !  1005- 1640 cm-1 
        dd(20,i)=conh2o*pss**1.0443*tss**(-2.4359) !  1645- 2530 cm-1 
        dd(21,i)=conh2o*pss**0.9681*tss**(-1.9537) !  2535- 3420 cm-1 
        dd(22,i)=conh2o*pss**0.9555*tss**(-1.5378) !  3425- 4310 cm-1 
        dd(23,i)=conh2o*pss**0.9362*tss**(-1.6338) !  4315- 6150 cm-1 
        dd(24,i)=conh2o*pss**0.9233*tss**(-0.9398) !  6155- 8000 cm-1 
        dd(25,i)=conh2o*pss**0.8658*tss**(-0.1034) !  8005- 9615 cm-1 
        dd(26,i)=conh2o*pss**0.8874*tss**(-0.2576) !  9620-11540 cm-1 
        dd(27,i)=conh2o*pss**0.7982*tss**( 0.0588) ! 11545-13070 cm-1 
        dd(28,i)=conh2o*pss**0.8088*tss**( 0.2816) ! 13075-14860 cm-1 
        dd(29,i)=conh2o*pss**0.6642*tss**( 0.2764) ! 14865-16045 cm-1 
        dd(30,i)=conh2o*pss**0.6656*tss**( 0.5061) ! 16340-17860 cm-1 
c
c for o3
c    cono3 in units of (1.e5 loschmidts)=(1 atm cm/km)=(2.6868e24 cm-3)
c
        cono3 = conjoe*wair*vfo3
        dd(31,i)=cono3 *pss**0.4200*tss**( 1.3909) !     0-  200 cm-1
        dd(32,i)=cono3 *pss**0.4221*tss**( 0.7678) !   515- 1275 cm-1
        dd(33,i)=cono3 *pss**0.3739*tss**( 0.1225) !  1630- 2295 cm-1
        dd(34,i)=cono3 *pss**0.1770*tss**( 0.9827) !  2670- 2845 cm-1
        dd(35,i)=cono3 *pss**0.3921*tss**( 0.1942) !  2850- 3260 cm-1
c
c for co2
c
        conco2= conjoe*wair*vfco2
        dd(36,i)=conco2*pss**0.6705*tss**(-2.2560)
        dd(37,i)=conco2*pss**0.7038*tss**(-5.0768)
        dd(38,i)=conco2*pss**0.7258*tss**(-1.6740)
        dd(39,i)=conco2*pss**0.6982*tss**(-1.8107)
        dd(40,i)=conco2*pss**0.8867*tss**(-0.5327)
        dd(41,i)=conco2*pss**0.7883*tss**(-1.3244)
        dd(42,i)=conco2*pss**0.6899*tss**(-0.8152)
        dd(43,i)=conco2*pss**0.6035*tss**( 0.6026)
c
c for co
c
        conco = conjoe*wair*vfco
        dd(44,i)=conco *pss**0.7589*tss**( 0.6911)
        dd(45,i)=conco *pss**0.9267*tss**( 0.1716)
c
c for ch4
c
        conch4= conjoe*wair*vfch4
        dd(46,i)=conch4*pss**0.7139*tss**(-0.4185)
c
c for n2o
c
        conn2o= conjoe*wair*vfn2o
        dd(47,i)=conn2o*pss**0.3783*tss**( 0.9399)
        dd(48,i)=conn2o*pss**0.7203*tss**(-0.1836)
        dd(49,i)=conn2o*pss**0.7764*tss**( 1.1931)
c
c for o2
c
        cono2 = conjoe*wair*vfo2
        dd(50,i)=cono2 *pss**1.1879*tss**( 2.9738) !     0-  265 cm-1
        dd(51,i)=cono2 *pss**0.9353*tss**( 0.1936) !  7650-52710 cm-1
c
c for nh3
c
        connh3= conjoe*wair*vfnh3
        dd(52,i)=connh3*pss**0.8023*tss**(-0.9111)
        dd(53,i)=connh3*pss**0.6968*tss**( 0.3377)
c
c for no
c
        conno = conjoe*wair*vfno
        dd(54,i)=conno *pss**0.5265*tss**(-0.4702)
c
c for no2
c
        conno2= conjoe*wair*vfno2
        dd(55,i)=conno2*pss**0.3956*tss**(-0.0545)
c
c for so2
c
        conso2= conjoe*wair*vfso2
        dd(56,i)=conso2*pss**0.2943*tss**( 1.2316)
        dd(57,i)=conso2*pss**0.2135*tss**( 0.0733)
c
c   herzberg continuum pressure dependence calculation, shardanand 1977
c      and   yoshino et al 1988
c
        dd(58,i)=(1.+.83*f1)*cono2
c
        if(lun.ne.0) then
          if(i.eq.1) write(lun,800)
          write(lun,'(i4,0pf9.2,f9.3,f7.1,1x,1p11e9.2)')
     &       i,z(i),p(i),t(i), dd(17,i),dd(31,i),
     &       dd(36,i),dd(44,i),dd(46,i),dd(47,i),
     &       dd(50,i),dd(52,i),dd(54,i),dd(55,i),dd(56,i)
        endif

 100  continue
c
c absorption integrals
c
      do 120 j=1,nz
        do 110 k=1,mxq
          uu(k,j)=0.
 110    continue
 120  continue
c
c use scale height to get first step of integration
c
      scfac=exp(-1.)
      airwt=(2.66e-23*(vfo2+1.e6-vfn2)+2.24e-23*vfn2)*1.e-6
      sclht=1.e-5*1.38e-16*t(nz)/(airwt*980.67)
c     
      do 200 i=nz,1,-1
        if(i.eq.nz) then
          ztop=2*z(i)-z(i-1)
          ptop=p(i)**2/p(i-1)
          ttop=t(i)
        else
          ztop=z(i+1)
          ptop=p(i+1)
          ttop=t(i+1)
        endif
        dz=ztop-z(i)
        if(          tbar=.5*(ttop+t(i))
        else 
          dp=(          drho=(          tbar=dp/drho
        endif
        do 150 k=1,mxq
          den1=dd(k,i)
          if(i.eq.nz) then
            uutop=0.
            den2=dd(k,i)*scfac
          else
            uutop=uu(k,i+1)
            den2=dd(k,i+1)
          endif
          if(abs(den1-den2).le..001*den1 .or.
     &       min(den1,den2).eq.0. .or. k.eq.8)  then
            du=.5*dz*(den1+den2)
          else
            du=dz*(den1-den2)/log(den1/den2)
          endif
          if(k.eq.5) du5=du
          if(k.eq.8) du8=du
          if(k.eq.9) then
            tfac=(296.-tbar)/(296.-260.)
            tfac=max(0.,min(1.,tfac))
            
            uu(9,i)=uutop+du5*tfac
          elseif(k.eq.59) then
            uu(59,i)=uutop+.269*du8*(tbar-273.15)
          elseif(k.eq.60) then
            uu(60,i)=uutop+.269*du8*(tbar-273.15)**2
          else
            uu(k,i)=uutop+du
          endif
 150    continue
        if(lun.ne.0) then
          if(i.eq.nz) write(lun,900)
          write(lun,'(i4,0pf9.2,f9.3,f7.1,1x,1p11e9.2)')
     &       i,z(i),p(i),t(i), uu(17,i),uu(31,i),
     &       uu(36,i),uu(44,i),uu(46,i),uu(47,i),
     &       uu(50,i),uu(52,i),uu(54,i),uu(55,i),uu(56,i)
        endif

 200  continue
c
 800  format('1',/,'  atmospheric profiles',//,
     + '  (if a molecule has more than one band, then the data for the',
     + ' first band are shown.)'//
     1 '   i      z       p       t      h2o      o3       co2',
     2 '      co       ch4      n2o      o2       nh3      no ',
     3 '      no2      so2'/
     4 '         (km)    (mb)    (k) g/cm**2/km  (            ',
     5 '          atm cm/km                                   ',
     6 '                  )')
 900  format('1',/,'  absorption integrals',//,
     + '  (if a molecule has more than one band, then the data for the',
     + ' first band are shown.)'//
     1 '   i      z       p       t      h2o      o3       co2',
     2 '      co       ch4      n2o      o2       nh3      no ',
     3 '      no2      so2'/
     4 '         (km)    (mb)    (k) g/cm**2     (            ',
     5 '          atm cm                                      ',
     6 '                  )')

      return
      end