C------------------------------------------------------------------------------
      SUBROUTINE TTABLE
C------------------------------------------------------------------------------
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C------------------------------------------------------------------------------
C     TABULAR TEMPERATURE PROFILE THAT MAKES A CUBIC INTERPOLATION
C     BETWEEN POINTS TO INSURE A CONTINUOUS TEMPERATURE GRADIENT

C     Compile with options /c (compile only) and /4Yb (error tracing)
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION MAT
      DIMENSION HPC(250),FN2C(250),ALPHA(250),TTBETA(250),GAMM(250),
     1 DELTA(250),SLOPE(250),MAT(4,5)
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "RKAM" INSERTED HERE
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
C     COMMON DECK "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH
C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/ dum(10),MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
C
C     TEMPERATURE    200-224
      EQUIVALENCE (W(200),TMODEL),(W(201),TFORM),(W(202),TID)
C
C     COMMON DECK "B5" INSERTED HERE
      INTEGER TMX,TNTBL,TITBL,TFRMTBL,IDST(10)
      COMMON/B5/TMX,TNTBL(10),TITBL(10),TFRMTBL(10),TGP(262)
      EQUIVALENCE (TGP,IDST),(ANT,TGP(11))
      EQUIVALENCE (AN,TGP(11)),(HPC,TGP(12)),(FN2C,TGP(262))
      integer tpx,tqtbl(10),tltbl(10),tirmtbl(10)
      DATA RECOGT,NOC/6.d0,0/
      DATA AQT/0.d0/
      DATA TPX/2/
      DATA TQTBL/1,11,512,7*0/
      DATA TLTBL/1,250,8*0/
      DATA TIRMTBL/1,2,8*0/
C
      entry settmp
      ant = aqt
      tmx = tpx
      call imove(tntbl,tqtbl,10)
      call imove(titbl,tltbl,10)
      call imove(tfrmtbl,tirmtbl,10)
      call setptp
      return
C
      ENTRY ITEMP
C
      CALL IPTEMP
C
C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT
      IF(NOC.GT.0 .AND. ANT.LE.0.d0) RETURN
C
      IF(idint(RECOGT) .NE. idint(TMODEL))
     1    CALL RERROR('TEMP   ','WRNG MODEL',RECOGT)
C
      ANT=dINT(ANT)
      NOC=ANT/2
      IF(ANT.NE.2*NOC .OR. NOC.LE.1)
     1     CALL RERROR('TTABLE','BAD NUMBER',ANT)
      ANT=0.d0
C
      MODT(1)=6HTTABLE
      MODT(2)=TID
C
C
      SLOPE(1)=(FN2C(2)-FN2C(1))/(HPC(2)-HPC(1))
      SLOPE(NOC)=0.d0
      NMAX=1
      DO 6 I=2,NOC
      IF (FN2C(I).GT.FN2C(NMAX)) NMAX=I
      IF (I.EQ.NOC) GO TO 4
      DO 3 J=1,3
      M=I+J-2
      MAT(J,1)=1.d0
      MAT(J,2)=HPC(M)
      MAT(J,3)=HPC(M)**2
  3   MAT(J,4)=FN2C(M)
      CALL GAUSEL (MAT,4,3,4,NRANK)
      IF (NRANK.LT.3) GO TO 60
      SLOPE(I)=MAT(2,4)+2.d0*MAT(3,4)*HPC(I)
    4 DO 5 J=1,2
      M=I+J-2
      MAT(J,1)=1.d0
      MAT(J,2)=HPC(M)
      MAT(J,3)=HPC(M)*hpc(m)
      MAT(J,4)=HPC(M)*hpc(m)*hpc(m)
      MAT(J,5)=FN2C(M)
      L=J+2
      MAT(L,1)=0.d0
      MAT(L,2)=1.d0
      MAT(L,3)=2.d0*HPC(M)
      MAT(L,4)=3.d0*HPC(M)*hpc(m)
    5 MAT(L,5)=SLOPE(M)
      CALL GAUSEL (MAT,4,4,5,NRANK)
      IF (NRANK.LT.4) GO TO 60
      ALPHA(I)=MAT(1,5)
      TTBETA(I)=MAT(2,5)
      GAMM(I)=MAT(3,5)
    6 DELTA(I)=MAT(4,5)
      HMAX=HPC(NMAX)
      NH=2
      RETURN
C
   60 PRINT 6000, I,HPC(I)
 6000 FORMAT(' THE',I4,'TH POINT IN THE TEMPERATURE PROFILE HAS'
     1,' THE HEIGHT',F8.2,' KM, WHICH IS THE SAME AS ANOTHER POINT.')
      stop 'error in tabular temperature profile'
C
      ENTRY TEMP
C
      ant=noc
      IF(NOC.LE.0)
     1    CALL RERROR('TTABLE','BAD N VALUE',ant)
      H=R-EARTHR
      PTT=0.d0
      PTR=0.d0
      PTTH=0.d0
      PTPH=0.d0
      IF (H.GE.HPC(1)) GO TO 12
   11 NH=2
      T=FN2C(1)+SLOPE(1)*(H-HPC(1))
      PTR=SLOPE(1)
      RETURN
   12 IF (H.GE.HPC(NOC)) GO TO 18
      NSTEP=1
      IF (H.LT.HPC(NH-1)) NSTEP=-1
   15 IF (HPC(NH-1).LE.H.AND.H.LT.HPC(NH)) GO TO 16
      NH=NH+NSTEP
      GO TO 15
   16 T=(ALPHA(NH)+H*(TTBETA(NH)+H*(GAMM(NH)+H*DELTA(NH))))
      PTR=(TTBETA(NH)+H*(2.d0*GAMM(NH)+H*3.d0*DELTA(NH)))
      RETURN
   18 T=FN2C(NOC)
C
      CALL PTEMP
      RETURN
      END
