
C-----------------------------------------------------------------------------
      SUBROUTINE CTABLE
C-----------------------------------------------------------------------------
C     TABULAR TEMPERATURE PROFILE THAT MAKES A CUBIC INTERPOLATION
C     BETWEEN POINTS TO INSURE A CONTINUOUS TERMPERATURE GRADIENT 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION HPC(250),FN2C(250),ALPHA(250),TTBETA(250),GAMM(250),
     1 DELTA(250),SLOPE(250),MAT(4,5)
      DOUBLE PRECISION MAT
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      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 "CC" INSERTED HERE 
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
C     SOUND SPEED 150-174
      EQUIVALENCE (W(153),REFC)
C     COMMON DECK "B3" INSERTED HERE
      INTEGER CMX,CNTBL,CITBL,CFRMTBL,IDSC(10)
      COMMON/B3/CMX,CNTBL(10),CITBL(10),CFRMTBL(10),CGP(512)
      EQUIVALENCE (CGP,IDSC),(ANC,CGP(11)) 
      EQUIVALENCE (HPC,CGP(12)),(FN2C,CGP(262))
      integer cpx,cqtbl(10),cltbl(10),cirmtbl(10)
      DATA NOC/0/
      DATA AQC/0.0/
      DATA CPX/2/
      DATA CQTBL/1,11,512,7*0/
      DATA CLTBL/1,250,8*0/
      DATA CIRMTBL/1,2,8*0/
C
	entry setspd
	anc = aqc
	cmx = cpx
	call imove(cntbl,cqtbl,10)
	call imove(citbl,cltbl,10)
	call imove(cfrmtbl,cirmtbl,10)
	call setpsp
	return
C
      ENTRY ISPEED
C 
      CALL IPSPEED 
C
C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT
      IF(NOC.GT.0 .AND. ANC.LE.0.0) RETURN
C
C         PRINT *,(CGP(I),I=11,30)
      IF(W(150) .NE. 8.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 8 if you are using model CTABLE.'
      ENDIF
C
        ANC=AINT(ANC)
      NOC=ANC/2
      IF(ANC.NE.2*NOC .OR. NOC.LE.1) then
           PRINT*,' ERROR in CTABLE.  Invalid number of points = ',noc
           STOP   ' at entry ISPEED.'
      endif
      ANC=0.d0
C
      MODC(1)=6HCTABLE
      MODC(2)= w(152)
C           PRINT *,(HPC(I),I=1,10),(FN2C(I),I=1,10)
      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)*hpc(m)
  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.*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)
C     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 SPEED
C          PRINT *,'ENTER CTABLE'
      anc=noc
      IF(NOC.LE.0) then
           PRINT*,' ERROR in CTABLE.  Invalid number of points = ',noc
           STOP   ' at entry SPEED.'
      endif
      H=R-EARTHR 
      PCST=0.d0  
      PCSR=0.d0  
      PCSTH=0.d0 
      PCSPH=0.d0 
      IF (H.GE.HPC(1)) GO TO 12
   11 NH=2
      T=FN2C(1)+SLOPE(1)*(H-HPC(1))
        PCSR=2.d0*T*SLOPE(1)
        CS=T*T
        CALL PSPEED
C         PRINT *,'LEAVE CTABLE1 ',R,C,PCSR
      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))))
      PCSR=(TTBETA(NH)+H*(2.d0*GAMM(NH)+H*3.d0*DELTA(NH))) 
C            PRINT *,NH,ALPHA(NH),TTBETA(NH),HPC(NH),FN2C(NH)
         PCSR=2.d0*T*PCSR
         CS=T*T
         CALL PSPEED
C         PRINT *,'LEAVE CTABLE2 ',R,C,PCSR
      RETURN
   18 T=FN2C(NOC)
C
      CALL PSPEED  
         PCSR=2.d0*T*PCSR
         CS=T*T
C         PRINT *,'LEAVE CTABLE3 ',R,C,PCSR
      RETURN
         END


