C-----------------------------------------------------------------------------
      SUBROUTINE CSSPOK2
C-----------------------------------------------------------------------------
C     SOUND SPEED SQUARED IS MODELED AS A FUNCTION OF ANGLE 'ALPHA'    
C     ABOUT A HORIZONTAL LINE AT SPECIFIED HEIGHT AND LATITUDE. DEPENDE
C     ON 'ALPHA' IS AS A SEQUENCE OF LINEAR SEGMENTS JOINED BY HYPERBOL
C     TANGENTS.                                                        
      implicit double precision (a-h,o-z)
      DOUBLE PRECISION ALC(20),Z(19),BIN(19),B(19),DL(19),LAM0
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 "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))
C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH 
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA  
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10 
      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1))
C     SOUND SPEED 150-174      
      EQUIVALENCE (W(153),REFC)
      EQUIVALENCE (Z0,CGP(12)),(C0,CGP(32)),(DL0,CGP(52))
      EQUIVALENCE (Z,CGP(13)),(BIN,CGP(33)),(DL,CGP(53)) 
      EQUIVALENCE (REFH,W(154)), (LAM0,W(155))           
      INTEGER CPX    ,CQTBL(10) ,CLTBL(10),CIRMTBL(10)   
      DATA N /0/
      DATA PFST,PFSTH,PFSPH/3*0.0/
      DATA AQC/0.0/
      DATA CPX/2/  
      DATA CQTBL/1,11,72,7*0/ 
      DATA CLTBL/1,20,8*0/ 
      DATA CIRMTBL/1,2,8*0/
C 
C 
      ENTRY SETSPD
      PCST=PFST
      PCSTH=PFSTH
      PCSPH=PFSPH
      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
      IF(REFC.GT.0) CREF=REFC   
      CALL IPSPEED 
C     IF HAD PREVIOUS CALL BUT NOTHING THIS TIME, EXIT NOW
C     RETAINING PREVIOUS TABULAR DATA COUNT 
      IF(N.GT.0 .AND. ANC.EQ.0.d0) RETURN 
      if (w(150) .ne. 4) then
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 4 if you are using model CSSPOK2.'
      endif

      MODC(1)=8HCSSPOK2 
      MODC(2)= w(152)
      N=ANC/3           
      IF(ANC.NE.3*N.OR.N.LE.0) then
           PRINT*,' ERROR in CSSPOK2: Invalid number of points = ',N
           STOP   ' at entry ISPEED.'
      endif
      N=N-2
      ANC=0.d0
      TH0=PID2-LAM0 
      R0=EARTHR+REFH 
C     CONVERT 'C' ARRAY INPUT(OVERLAYS 'B' ARRAY) TO 'B' ARRAY
      ZM1=Z0 
      CS0=C0*C0
      CSM1=CS0 
      NP1=N+1  
      DO 10 I=1,NP1
         ZR=Z(I)   
         ALC(I)=ALCOSH((ZR-Z0) / DL(I)) 
         CS=BIN(I)*bin(i)
         B(I)=(CS-CSM1)/(ZR-ZM1)
         ZM1=ZR 
10       CSM1=CS
      RETURN    
C 
      ENTRY SPEED 
      anc=n
      IF(N.LE.0) then
           PRINT*,' ERROR in CSSPOK2.  Invalid number of points: ',anc
           STOP   ' at entry SPEED.'
      endif
      X=(R-R0)/(EARTHR*(TH0-TH)) 
      ZR=dATAN2(R-R0,EARTHR*(TH0-TH))
      SUM=0.d0 
      PCSZ=B(1)
      DO 1 I = 1, N
      SAV=0.5d0*(B(I+1)-B(I))  
      PCSZ= PCSZ+ SAV * (1.d0+DTANH ((ZR-Z(I)) /DL(I))) 
   1  SUM = SUM+DL(I) * SAV *(ALCOSH((ZR-Z(I))/DL(I))-ALC(I))
      CS = CS0+SUM + 0.5d0*(B(1) + B(N + 1)) * (ZR-Z0)
      IF(CS.LT.0.d0) CALL PMDSTOP 
      F=X/(1.d0+X*X)
      PZR=F/(R-R0)  
      PZTH=-F/(TH-TH0)
      PCSR=PCSZ*PZR   
      PCSTH=PCSZ*PZTH 
      CALL PSPEED     
      RETURN          
      END             


