C-----------------------------------------------------------------------------
      SUBROUTINE CSWAV2
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,lz,lx
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
      COMMON/WW/ dum(10), MAXW,W(1000)
      EQUIVALENCE (W(1),EARTHR)
      EQUIVALENCE (CR,W(185)),(TIM,W(184))
      EQUIVALENCE (Z0,W(178)), (H,W(179)),(AMPL,W(180))
      EQUIVALENCE (VX,W(181)),(LX,W(182)),(LZ,W(183))
      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. 6.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(175) should be 6 if you are using model CSWAV2.'
      ENDIF

      MODC(3)=6HCSWAV2
      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
      DEXPO=0.d0
      EXPO=-(DZ*DZ*FWZ)
      IF(EXPO .GT. -200.d0) DEXPO=dEXP(EXPO)
C DEXPO IS THE HEIGHT-AMPLITUDE FACTOR
      TMP=PIT2*(TIM+(PID2-TH)*EARTHR/LX+Z/LZ)
      SINW=dSIN(TMP)
      COSW=dSIN(PID2-TMP)
C COSW IS THE SINUSOIDAL SPACE-TIME FACTOR
C
      F=AMPL*DEXPO*COSW
      PFR=F*(-2.d0*DZ*FWZ-SINW/COSW*PIT2/LZ)
      PFTH=F*(SINW/COSW*PIT2*EARTHR/LX)
      PFPH=0.d0
      PFT=-PIT2*VX/LX*CS*AMPL*DEXPO*SINW
C
      CS0=CS
      CSR=CR*CR
C   CS0 IS THE BACKGROUND CS PROFILE
C   CR IS A REFERENCE SOUND SPEED
  
      DENOM=1.d0-2.d0*CS0/CSR*F
      GORK=2.d0/CSR/DENOM
      CS=CS0/DENOM
      PCSR= CS*(PCSR /CS0+GORK*(F*PCSR +CS0*PFR ))
      PCSTH=CS*(PCSTH/CS0+GORK*(F*PCSTH+CS0*PFTH))
      PCSPH=CS*(PCSPH/CS0+GORK*(F*PCSPH+CS0*PFPH))
      PCST =CS*(PCST /CS0+GORK*(F*PCST +CS0*PFT ))
      IF (VX.NE.0d0) PCST= CS*(PCST /CS0+GORK*(F*PCST +CS0*PFT ))
C  CS IS THE SQUARED SOUND SPEED, AND PCS* ARE ITS SPACE/TIME DERIVATIVES
  
      RETURN
      END



