subroutine opthik(x,p,t,tl,cofx,tmp0,tmpcof,tmp101)
c
c***********************************************************************
cccc
c     subroutine opthik
c
c     purpose-
c        divide the atmosphere into layers and calculate their
c        absorption,scattering and total optical depths measured from
c        top of the atmosphere.
c
c     method-
c        uses spline interpolation to calculate optical thickness values
c        on finer grid points than those obtained from the input data
c
c     calling sequence-  call opthik(x,p,t,tl,cofx)
c
c        variable   type   i/o  description
c        --------   ----   ---  -----------
c
c        x(101)     r*8     i   cumulative ozone thicknesses
c
c        p(101)     r*8     o   pressure at each layer (atmospheres)
c        t(101)     r*8     o   total optical depth of each layer
c        tl(101)    r*8     o   log(t)
c        cofx(4,487) r*8    o   spline interpolation coeff.
c
c     this subroutine also returns(thru' the common 'thkns' )
c
c        tsl(101)    r*8    o  log(scattering optical depth)
c        tts(202)    r*8    o  interpolated scattering opt. depth
c        tt(202)     r*8    o  interpolated total opt. depth
c        ttl(202)    r*8    o  log(tt)
c
c
c     external references-
c        splset
c     last modified......11/4/94   zia ahmad
c     purpose............to correct the raleigh optical thickness
c            for 1/r**2 effect of gravity
c     last modified......03/10/95  dave flittner
c     purpose............to set radius of earth equal to value in consts
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     to impliment the gravity correction to the rayleigh scattering
c     optical depth.
c
cccc
c***********************************************************************
c
      implicit integer*4(i-n),real*8 (a-h,o-z)
      real *8 p(101),t(101),tl(101),x(101),cofx(4,487),tmpcof(2),
     & tmp0,tmp101(101)
c
      include "prints.inc"
      include "atmos.inc"
      include "consts.inc"
      include "thkns.inc"
      include "cgcorrect.inc"
c
c*****divide the atmosphere into 101 layers by a linear interpolation
c*****of log(p)            (pmax=1, pmin=1.0e-5)
c*****store p and calculate log(ts),t and log(t)
c
      elg = dlog(10.d+00)
      dum2a = 0.05*elg
      betal = dlog(beta)
      do 10 i = 1, 101
      plog = dum2a*float(i-101)
      p(i) = dexp(plog)
c
c     determine the height above the ground(p=1.0) using a constant
c     scale height of H=6.95 and the g(h).  Assume the radius of earth
c     as: r0=r (in commmon block consts) and then compute the ratio: g0/g(h)
c     the value of H is set in iniclz and is divided by the radius of   !def
c     the earth.                                                        !def
      if(lgcorrect)then                                        !def
        gr=(1.0d0-pscaleforg*dlog(p(i)))**2                      !def
        grl=dlog(gr)
      else                                                    !def
        gr=1.0d0                                              !def
        grl=0.0d0                                             !def
      endif                                                   !def
c
      deltmp=tmp101(i) - tmp0
      alpha=alpha0+tmpcof(1)*deltmp+tmpcof(2)*deltmp**2
      ta = alpha*x(i)
C      ts = beta*p(i)
      ts = beta*p(i)*gr
      t(i) = ta + ts
      tl(i) = dlog(t(i))
C      tsl(i) = betal + plog 
      tsl(i) = betal + plog + grl
      if (jprint(5) .eq. 1)
     1 write(33,101) i,x(i),p(i),t(i),tsl(i),tl(i),tmp101(i)
10    continue
c
c*****fit a spline between log(ts) and log(t)
c
      call splset (tsl,tl,cofx,101)
c
c
c*****increase the no. of layers to 202 by preferentially adding
c*****layers at the atmosphere bottom using a quadratic function.
c
      cnsta = -3.4219e-3
      cnstb = 6.7056e-2
      tt(1) = 0.
      tts(1) = 0.
      jb = 2
      do 20 i = 2, 202
      adum = float(202-i)
      ttsl = betal + cnsta*adum*(1. + cnstb*adum)
      tts(i) = dexp(ttsl)
c
c*****search for two layers in the old atmos. stradlling a given
c*****layer in the new atmos. and find tt by spline interpolation
c
      do 15 j = jb, 101
      if (ttsl .gt. tsl(j)) go to 15
      dum1 = tsl(j) - ttsl
      kb = j - 1
      dum2 = ttsl - tsl(kb)
      dum3 = dum1*(cofx(1,kb)*dum1**2 + cofx(3,kb)) + dum2*(cofx(2,kb)*
     1dum2**2 + cofx(4,kb))
      ttl(i) = dum3
      tt(i) = dexp(dum3)
      go to 20
   15 continue
   20 continue
      return
101   format(i4,6d15.4)
      end