c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/opthik.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine opthik(p,t,tl,cofx,tmp0,tmpcof)
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
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 input via atm_101
c        x(101)     r*8     i   cumulative ozone thicknesses
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
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn.
c
c last mod. Oct. 14, 1998...def
c purpose: add the standard profile option (use temp. dist. w/o weighting).
c
c last mod. Nov. 23, 1998...def
c purpose: check to make sure that tt decreases with altitude.
c
c last mod. Apr. 13, 1999...def
c purpose: call getgasprf just to load the values into totgas.
c
c     last mod. May 3, 1999...def
c     purpose: add option to compute addition absorbing gas profile using the
c     internal, implicit neutral density profile and a constant mixing ratio.
c     Currently it is set to do only o2 or o4 absorption and is controled by
c     the lo2abs and lo4abs logicals read in the ENV file.
c
c last mod. Nov 2, 1999...def
c purpose: use .not. instead of .eqv. in logical testing.
c
c last mod. Apr. 14, 2000...def
c purpose: add switch for using weighted temperatures with user input profile
c
cccc
c***********************************************************************
c
      implicit none
c passed
      real*8 p(101),t(101),tl(101),cofx(4,487)
      real*8 tmpcof(2),tmp0
c commons
      include 'switches.cmn'
      include 'atmos.cmn'
      include 'consts.cmn'
      include 'thkns.cmn'
      include 'cgas2.cmn'
      include 'atm_101.cmn'
c local
      integer igas,i,jb,kb,j
      real*8 taother,elg,dum2a,betal,plog,gr,grl,cnsta,cnstb
      real*8 adum,dum1,dum2,dum3,ttsl,deltmp,alpha,ta,ts
      real*8 interp_gc
      external interp_gc
c
      if (jprint(5) .eq. 1)write(33,*)'opthik. debugging output------'
      if (jprint(5) .eq. 1)write(33,*)'i,x_101(i),p(i),t(i),tsl(i),',
     &'tl(i),t_101(i),t(i)-ts(i)'
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
      if((ngas.gt.1).and.(ildgas.eq.0))then
c
c Load other absorbing gas profiles and absorbtion coefficients.
c Then calc the cumlative from top of atmo, pressure=1.0e-05 atm
c to each altitude.  Then convert to pressure scale with 1 atm at 
c the surface.
c
        do igas=1,ngas-1
           call getgasprf(igas)
        enddo
        ildgas=1
      endif
      betal = dlog(beta)
      gr=1.0d0
      grl=0.0d0
      ta=0.0d0
      taother=0.0d0
      do 10 i = 1, n_101
         if(lv7tab)then
           plog = dum2a*float(i-n_101)
           p(i) = dexp(plog)
         else
           plog=p_101lg(i)
           p(i)=p_101(i)
         endif
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
c     the earth.
         if(gc_type.gt.0)then
           if(gc_type.eq.1)then
             gr=(1.0d0-pscaleforg*dlog(p(i)))**2
           else
             gr=gc_factor(i)
           endif
           grl=dlog(gr)
         endif
c
         deltmp=t_101(i) - tmp0
         alpha=alpha0+tmpcof(1)*deltmp+tmpcof(2)*deltmp**2
         if(lwgttmp)then
           ta = alpha*x_101(i)
         else
           if(i.eq.1)then
             ta = alpha*x_101(i)
           else
             ta = ta+alpha*(x_101(i)-x_101(i-1))
           endif
         endif
         if(ngas.gt.1)then
c
c add the other gases
c
            taother=0.0d0
            do igas=1,ngas-1
               taother=taother+totgas(i,igas)*alphagas(igas)
            enddo
         endif
         ts = beta*p(i)*gr
         t(i) = ta + ts + taother
         tl(i) = dlog(t(i))
         tsl(i) = betal + plog + grl
         if (jprint(5) .eq. 1)
     1write(33,101) i,x_101(i),p(i),t(i),tsl(i),tl(i),t_101(i),
     2t(i)-exp(tsl(i))
c         if(i.gt.1)then
c           if(t(i).lt.t(i-1))then
c      write(6,*)'opthik. t decreases'
c      write(6,101) i-1,x_101(i-1),p(i-1),t(i-1),tsl(i-1),tl(i-1),
c     &t_101(i-1)
c      write(6,101) i,x_101(i),p(i),t(i),tsl(i),tl(i),t_101(i)
c           endif
c         endif
10    continue
      if(ngas.gt.1)then
        do igas=1,ngas-1
           write(6,1000)igas+1,alphagas(igas)*totgas(n_101,igas)
1000  format(' opthik. Absorption vertical optical depth ',
     &'for gas ',i2,'=',1pe10.3)
        enddo
      endif
c
c*****fit a spline between log(ts) and log(t)
c
      call splset (tsl,tl,cofx,n_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
      pressure(1)=0.0d0
      gr=1.0d0
      grl=0.0d0
      do 20 i = 2, 202
         adum = float(202-i)
         pressure(i)=exp(cnsta*adum*(1. + cnstb*adum))
         if(.not.lv7tab)then
           if(gc_type.gt.0)then
             if(gc_type.eq.1)then
               gr=(1.0d0-pscaleforg*dlog(pressure(i)))**2
             else
               gr=interp_gc(pressure(i))
             endif
             grl=dlog(gr)
           endif
         endif
         ttsl = betal + cnsta*adum*(1. + cnstb*adum) + grl
         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, n_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)) + 
     1dum2*(cofx(2,kb)*dum2**2 + cofx(4,kb))
            ttl(i) = dum3
            tt(i) = dexp(dum3)
            if(i.gt.2)then
c check to see if tt decreases.  If so, then use linear interpolation
              if(tt(i).lt.tt(i-1)) then
                 ttl(i)=
     &tl(j)*dum2/(tsl(j)-tsl(kb))+tl(kb)*dum1/(tsl(j)-tsl(kb))
              endif
            endif
            go to 20
   15    continue
   20 continue
      return
101   format(i4,1P7e15.4)
      end