C-----------------------------------------------------------------------------
      SUBROUTINE CSWAVE
C-----------------------------------------------------------------------------
      implicit double precision (a-h,o-z)
C     SOUND SPEED PERTURBATION MODEL
C     WITH SINUSOIDAL VARIATION VERTICALLY, HORIZONTALLY AND WITH TIME
C     AND GAUSSIAN AMPLITUDE DEPENDENCE ON DEPTH.
C     WAVE TRAVELS FROM THE NORTH POLE TOWARD THE SOUTH POLE.
C     SPECIFY LATITUDE INSTEAD OF CO-LATITUDE.
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 "B4" INSERTED HERE
      INTEGER DCMX,DCNTBL,DCITBL,DCFRMTB,IDSDC(10)
      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(10)
      EQUIVALENCE (DCGP,IDSDC)
C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC,lx,lz
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (W(1),EARTHR)
      EQUIVALENCE (Z0,W(178)), (H,W(179)),(AMPL,W(180))
      EQUIVALENCE (VX,W(181)),(LX,W(182)),(LZ,W(183)),(TIM,W(184))
      INTEGER DFMX,DFNTBL(10),DFITBL(10),DFFRMTB(10)
      DATA DFMX/1/
      DATA DFNTBL/1,11,8*0/
      DATA DFITBL/1,9*0/
      DATA DFFRMTB/1,9*0/
C
C
      ENTRY SETPSP
      DCMX=DFMX
      CALL IMOVE(DCNTBL,DFNTBL,10)
      CALL IMOVE(DCITBL,DFITBL,10)
      CALL IMOVE(DCFRMTB,DFFRMTB,10)
      RETURN
C
      ENTRY IPSPEED
      IF(W(175) .NE. 5.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(175) should be 5 if you are using model CSWAVE.'
      ENDIF

      MODC(3)=6HCSWAVE
      MODC(4)= w(177)
C IF GAUSSIAN WIDTH IN HEIGHT = 0, SET HEIGHT DEPENDENCE = 0
      FWZ=0.d0
      IF (H.NE.0.d0) FWZ=1.d0/H/H
      RETURN
C
      ENTRY PSPEED
      IF(AMPL.EQ.0.d0) RETURN
      Z=R-EARTHR
      DZ=Z-Z0
C  Z IS HEIGHT OF RAY ABOVE EARTH'S SURFACE
C  DZ IS HEIGHT OF RAY ABOVE HEIGHT OF WAVE AMPLITUDE MAXIMUM
C
      DEXPO=0.d0
      EXPO=-(DZ*DZ*FWZ)
      IF(EXPO .GT. -200.d0) DEXPO=dEXP(EXPO)
C DEXPO IS THE HEIGHT-AMPLITUDE FACTOR
C
      TMP=PIT2*(TIM+(PID2-TH)*EARTHR/LX+Z/LZ)
      SINW=dSIN(TMP)
      COSW=dSIN(PID2-TMP)
C COSW IS THE SINUSOIDAL SPACE-TIME FACTOR
      FACTOR=1.d0+AMPL*DEXPO*COSW
      PCSR=PCSR*FACTOR-CS*AMPL*DEXPO*(2.d0/H/H*DZ*COSW+PIT2/LZ*SINW)
      PCSTH=PCSTH*FACTOR+CS*AMPL*PIT2*EARTHR/LX*SINW*DEXPO
      PCSPH=PCSPH*FACTOR
      PCST=0.d0
      IF (VX.NE.0.d0) PCST=-PIT2*VX/LX*CS*AMPL*DEXPO*SINW
      CS=CS*FACTOR
C CS IS THE SQUARED SOUND SPEED, AND PCS ARE ITS SPACE/TIME DERIVATIVES
C
      RETURN
      END



