subroutine dexpk1(a,b)
c
c***********************************************************************
c
cccc
c     subroutine dexpk1
c
c     purpose-
c        calculate the elements of k1 matrix averaged over the interval
c        a-b .  (dave's eq. 4.11)
c
c     method-
c        the averages are calculated using eq(5.1.26) in abramowitz and
c        segun.
c
c     calling sequence -  call dexpk1(a,b)
c
c
c        variables     type   i/o  description
c        ---------     ----   ---  -----------
c        a             r*8     i   abs(a-b) is the argument of the
c        b             r*8     i   exponential  integrals
c
c     result of the subroutine is returned in common block  'es'
c        e(6)         r*8      o   exponential  integrals
c        eek(6)      r*8       o   elements of the k-matrix
c
c     external references -   dexpi
c
c     modifications-  by  k.klenk  8/19/77
c
c        functions k1(1) and k1(2) are also calculated and are returned
c        in eek(4) and eek(5) respectively (dave's eqn. 4.7)
c
c     last modified by zia ahmad 9/10/93
c     purpose: to include the effect of molecular anisotropy
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.
c
c***********************************************************************
c
cccc
c
      implicit integer*4(i-n),real*8 (a-h,o-z)
      real *8 de(6)
      include "consts.inc"
      include "es.inc"
c
c     new statemnet
      include "depolt.inc"
c     end of new statement
c
c
c*****same as dexpk- calculate expon. integrals
c
      dt = dabs(a-b)
      de(1) = - dexpi(-dt)
      dum = dexp(-dt)
      do 70 j = 2, 6
   70 de(j) = odan(j)*(dum-dt*de(j-1))
      do 80 i=1,6
   80 e(i) = de(i)
c
c*****calculate elements of k1 matrix (azimuth independent term)
c     by averaging over dt
c     also calculate k1(1) and k1(2) functions for azimuth
c     dependent terms
c     new statements added here
      if(ipol.eq.0)then
          eek(1) = 0.375 *(1.2d+00-de(2)-de(6))/dt
          eek(2) = c38sq2*(2.0d+00/15.d+00-de(4)+de(6))/dt
          eek(3)= 0.75d+00*(8.d+00/15.d+00-de(2)+2.d+00*de(4)-de(6))/dt
          eek(4)=(c1415-e(2)-e(4)+2*e(6))*0.375d0/dt
          eek(5)=(c2815-e(2)-2*e(4)-e(6))*0.1875d0/dt
      else
          eek(1)=q12s*((sdp**2+q**2)*(1.0d0-e(2))+
     1           2.0d0*q*(1.0d0/3.0d0-e(4))+
     2           (1.0d0/5.0d0-e(6)) )/dt
          eek(2)=q12s*delp*(q*(1.0d0-e(2))+
     1           (1.0d0-q)*(1.0d0/3.0d0-e(4))-
     2           (1.0d0/5.0d0-e(6)))/dt
          eek(3)=q12s*delp**2*((1.0d0-e(2))-
     1           2.0d0*(1.0d0/3.0d0-e(4))+
     2           (1.0d0/5.0d0-e(6)))/dt
          eek(4)=q2*0.375d0*(c1415-e(2)-e(4)+2.0d0*e(6))/dt
          eek(5)=q2*0.1875d0*(c2815-e(2)-2.0d0*e(4)-e(6))/dt
      endif
c     end of new statements 
      return
      end