subroutine expone (nc,ncp1,imu)
c
c***********************************************************************
cccc
c     subroutine expone
c     version aug 18,1977
c
c     purpose
c
c        expone is a fortran 4 routine to calulate the iteration
c        integrals of eq 4.11 and k functions of eq 4.6
c
c     method
c        kernals of intgrals to be used in reflex and iterate
c        are calulated and stored on a disk file
c
c     calling sequence
c        call expone (ncp,ncp1)
c
c        variable   type    i/o  description
c        --------   ----    ---  ------------
c
c        nc         i*4     i    # of layers from top of
c                                atmosphere to reflecting surface
c        ncp1       i*4     o    nc+1
c
c     external references
c        dexpk1
c
c     common areas referenced
c        consts
c        eks
c        emm
c        es
c        prints
c        thkns
c
c     analysis and programming k.f. klenk sasc aug 18,1977
c
c     modifications (date , name, purpose)
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
cccc
c***********************************************************************
c
      implicit integer*4(i-n),real*8(a-h,o-z)
      real*8 ekary(202,5),extmux(202,15)
c
      include "prints.inc"
      include "kmat.inc"
      include "emm.inc"
      include "eks.inc"
      include "thkns.inc"
      include "consts.inc"
      include "es.inc"
c
c     new statement
      include "depolt.inc"
c     end of new statement
c
      do 75 i=1,ncp1
c        initialize ekary
         do 200 j=1,5
           ekary(i,j)=0.d0
  200    continue
c
c        calculate w(t)dt*attenuation factor (see eqn. 3.6,6.7)
c
      do 205 j=1,imu
      extmu(i,j)=dtsp(i)*dexp(-tt(i)/emu(j))
      extmux(i,j)=extmu(i,j)/emu(j)
  205 continue
      if (i.eq.1) go to 55
      im1=i-1
c
      do 210 j=1,5
      ekary(im1,j)=0.
  210 continue
c
      if(i.eq.2) go to 50
      im2=i-2
c
c        calculates the elements of k1 matrix (see eqn 4.11)
c        evaluate k1(1),k1(2) of eqn 4.7 for current value of
c        argument (tt(i)-tt(j))
c
      do 45 j=1,im2
      call dexpk(tt(i),tt(j))
c
c        ekary array contains k matrix*w(t)dt
c
      do 220 jk=1,5
      ekary(j,jk)=dtsp(j)*eek(jk)
  220 continue
      if(i.lt.ncp1) go to 45
c     new statements
      if(ipol.eq.0)then
          ek4(j)=(e(2)+e(4))*dtsp(j)
          ek5(j)=sq2*(e(2)-e(4))*dtsp(j)
      else
          ek4(j)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(j)
          ek5(j)=delp*(e(2)-e(4))*dtsp(j)
      endif 
c     end of new statements
   45 continue
c
      call dexpk(tt(i),tt(im1))
      do 230 k=1,5
      ekary(im1,k)=dtts(im1)*eek(k)
  230 continue
c
      if(i.lt.ncp1) go to 50
c     new statements
      if(ipol.eq.0)then
          ek4(nc)=dtsp(nc)*(e(2)+e(4))
          ek5(nc)=sq2*dtsp(nc)*(e(2)-e(4))
      else
          ek4(nc)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(nc)
          ek5(nc)=delp*(e(2)-e(4))*dtsp(nc)
      endif
c     end of new statements
   50 call dexpk1(tt(i),tt(im1))
c
      do 240 k=1,5
      ekary(im1,k)=ekary(im1,k)+dtts(i)*eek(k)
  240 continue
c
      do 250 k=1,5
      ekary(i,k)=dtts(i)*eek(k)
  250 continue
      if(i.lt.ncp1) go to 55
c     new statemnets
      if(ipol.eq.0)then
          ek4(ncp1)=dtsp(ncp1)*(odan(2)+odan(4))
          ek5(ncp1)=sq2*dtsp(ncp1)*(odan(2)-odan(4))
      else
          ek4(ncp1)=((1.0d0+2.0d0*q)*odan(2)+odan(4))*dtsp(ncp1)
          ek5(ncp1)=delp*(odan(2)-odan(4))*dtsp(ncp1)
      endif
c     end of new statements
      go to 65
   55 ip1=i+1
      call dexpk1(tt(i),tt(ip1))
c
      do 260 k=1,5
      ekary(i,k)=ekary(i,k)+dtts(ip1)*eek(k)
  260 continue
c
      do 270 k=1,5
      ekary(ip1,k)=dtts(ip1)*eek(k)
  270 continue
c
      if(i.eq.nc) go to 65
      ip2=i+2
      call dexpk(tt(i),tt(ip1))
c
      do 280 k=1,5
      ekary(ip1,k)=ekary(ip1,k)+dtts(ip2)*eek(k)
  280 continue
c
      do 60 j=ip2,ncp1
      call dexpk(tt(i),tt(j))
c
      do 290 k=1,5
      ekary(j,k)=dtsp(j)*eek(k)
  290 continue
   60 continue
c
c
   65 continue
      do 5000 ii = 1,202
         do 5001 iii = 1,5
            nek(ii,iii,i) = ekary(ii,iii)
5001     continue
5000  continue
   75 continue
c
c        optionally print calculated exponential terms
c        (zstar matrix * w(t)dt see eqn 6.3)
c
      if (jprint(3).ne.0) then
         write (33,6120) (j,ek4(j),ek5(j) ,j=1,ncp1)
         do 295  i=1,imu
            thta=dacos(emu(i))/cnvrt
            write(33,7120) thta,(j,extmux(j,i),j=1,ncp1,10)
            write(33,7130) ncp1,extmux(ncp1,i)
295      continue
      endif
      return
6120  format(1h1,t50,'debug print out for expone',///,
     1 1h ,2('  j',4x,'ek4(j)',6x,'ek5(j)',6x),//,2(1x,i3,2d12.4,3x))
7120  format(1x,//,1x,'theta=',f6.2,' *** attenuation factors/*** ',//,
     1 1x,'level',5x,'extmu/mu',/,(1x,i5,3x,d12.6))
7130  format(1x,i5,3x,d12.6)
      end