c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/dtaus.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      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.
c
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn.
c
c     last mod. Oct. 19, 1998...def
c     purpose: add label for debugging output.
c
cccc
c***********************************************************************
c
      implicit none
      include 'parameter.inc'
c passed in
      real*8 cofx(4,487),xpt
c commons
      include 'thkns.cmn'
      include 'contrl.cmn'
      include 'atmos.cmn'
      include 'consts.cmn'
      include 'switches.cmn'
c passed out
      integer*4 nc,ncp1
c
      real*8 interp_gc
      external interp_gc
c local
      integer*4 i,k
      real*8 dum1,dum2,dum3,pnlog,pnotb
c
c*****calculate dtsp and dtts and find nc
c
c

      if(gc_type.gt.0)then
        if(gc_type.eq.1)then
          pnotb = pnot * beta * (1.0d0-pscaleforg*dlog(pnot))**2     !def
        else
          pnotb = pnot * beta * interp_gc(pnot)
        endif
      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
      if (jprint(8).eq.1) write(33,6000)
      if (jprint(8).eq.1) 
     1 write(33,6100)ncp1,pnot,(i,tt(i),tts(i),dtts(i),
     2 dtsp(i),i=1,ncp1)
6000  format('dtaus. debugging output-----')
6100  format(//,6hincp1=,i3,5x,5hpnot=,f5.3/1h0,2(2x,1h1,8x,
     1 2htt,12x,3htts,10x,4hdtts,10x,4hdtsp,7x)//
     2 (1x,2(i3,4d14.4,4x)))
      return
      end