subroutine dtaus(cofx,nc,xpt,ncp1)
c
c***********************************************************************
c
cccc
c     subroutine  dtaus
c
c     purpose-
c        1. calculate the avg. optical thickness of each layer and of
c           each half layer at mid points between the layers
c        2. find the no. of layers between the atmosphere top and a
c           specified bottom (pnot)
c
c     input/output variables -
c
c        variables      type   i/o  description
c        ---------      ----   ---  -----------
c     formal parameters
c        cofx(4,487)    r*8     i   spline interpolation coeff.
c        nc             i*4     o   no. of layers between the top and
c                                   the bottom(not incl. bottom layer)
c        ncp1           i*4     o   nc+1
c     part of common 'thkns'
c         dtsp(202)     r*8     o   avg. optical thickness of each layer
c         dtts(202)     r*8     o   avg. optical thickness of each half
c                                   layer at mid points between layers
c
c     no external references
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.  Perform gravity correction when
c     computing pnotb.
c     Also use logical switch lgcorrect
c     to impliment the gravity correction to the rayleigh scattering
c     optical depth.
cccc
c***********************************************************************
c
      implicit integer*4(i-n),real*8 (a-h,o-z)
      real*8 cofx(4,487),xpt
      include "thkns.inc"
      include "contrl.inc"
      include "atmos.inc"
      include "hedout.inc"
      include "consts.inc"
      include "cgcorrect.inc"
c
c*****calculate dtsp and dtts and find nc
c
c

      if(lgcorrect)then
        pnotb = pnot * beta * (1.0d0-pscaleforg*dlog(pnot))**2     !def
      else
        pnotb = pnot * beta
      endif
      pnlog = dlog(pnotb)
      dtts(1) = 0.5*tts(2)
      dtsp(1) = dtts(1)
      do 25 i = 2, 202
      dtsp(i) = 0.5*(tts(i+1) - tts(i-1))
      dtts(i) = 0.5*(tts(i) - tts(i-1))
      if (tts(i) .lt. pnotb) go to 25
      ncp1 = i
      nc = i - 1
      go to 26
   25 continue
      ncp1 = 202
      nc = 201
c
c*****calculate tt and log(tt) of the bottom layer
c*****(located at pnot)- by spline interpolation
c
   26 do 35 i = 80, 101
      if (tsl(i) .lt. pnlog) go to 30
   28 k = i - 1
      dum1 = tsl(k+1) - pnlog
      dum2 = pnlog - tsl(k)
      dum3 = dum1*(cofx(1,k)*dum1**2 + cofx(3,k)) + dum2*(cofx(2,k)*
     1dum2**2 + cofx(4,k))
      ttl(ncp1) = dum3
      tt(ncp1) = dexp(dum3)
      tautot = tt(ncp1)
      tts(ncp1) = pnotb
      go to 40
   30 if (i .eq. 101) go to 28
   35 continue
c
c*****correct the dtts and dtsp of the bottom layer
c
   40 dtts(ncp1) = 0.5*(tts(ncp1) - tts(nc))
      dtsp(ncp1) = dtts(ncp1)
      dtsp(nc) = dtts(ncp1) + dtts(nc)
      alfaef = (tautot-pnotb)/xpt
c     --store data for archive tape
      rarray(1) = -12.0d0
      rarray(203) = -12.0d0
      do 100 i = 2,ncp1
         j = i+202
         rarray(i) = ttl(i)
         rarray(j)=dlog(tts(i))
100   continue
      rarray(405) = code
      rarray(406) = pnot
      rarray(408) = alfaef
      rarray(409) = beta
      do 300 i = 1,82
         psaray(i) = pshold(i)
300   continue
      do 200 i = 1,imuz
         j = 409 + i
         rarray(j) = thnot(i)
200   continue
      rarray(824) = wavlen
      iarray(1) = ncp1
      iarray(2) = imuz
      iarray(3) = lambda
c      if (jprint(8).eq.1) 
c     1 write(33,6100)ncp1,pnot,(i,tt(i),tts(i),dtts(i),
c     2 dtsp(i),i=1,ncp1)
c6100  format(//,6hincp1=,i3,5x,5hpnot=f5.3/1h0,2(2x,1h1,8x,
c     1 2htt,12x,3htts,10x,4hdtts,10x,4hdtsp,7x)//
c     2 (1x,2(i3,4d14.4,4x)))
      return
      end