C-----------------------------------------------------------------------------
      SUBROUTINE CSMUNK3  
C-----------------------------------------------------------------------------
C SPEED MODEL BASED ON THE 'CANONICAL' MODEL FOR A SOUND CHANNEL  
C DERIVED BY MUNK(J.ACOUSTICAL SOC. AM. 55,220-226).  
C THE FOUR PARAMETERS OF THE MODEL ARE ALLOWED TO VARY LINEARLY IN
C LONGITUDE, PH. THE USER SPECIFIES VALUES FOR THE FOUR PARAMETERS AT
C TWO LONGITUDES, PH1 AND PH2, AND THE PROGRAM INTERPOLATES LINEARLY
C TO GET THE VALUES AT OTHER LONGITUDES.
C
C CSMUNK3 ALLOWS DIFFERENT VALUES OF H AND EP ABOVE AND BELOW THE
C SOUND-CHANNEL AXIS. FIRST DERIVATIVES REMAIN CONTINUOUS.
C
C     C = CA(1+EP(ETA+EXP(-ETA)-1))
C        WHERE
C           ETA = -2(Z-ZA)/H
C           Z = R-EARTHR
      implicit double precision (a-h,o-z)
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)
      EQUIVALENCE (W(154),PH1), (W(155),CA1) ,(W(156),ZA1) 
      EQUIVALENCE (W(157),H1H) ,(W(158),EP1H)
      EQUIVALENCE (W(159),PH2), (W(160),CA2) ,(W(161),ZA2)
      EQUIVALENCE (W(162),H2H) ,(W(163),EP2H),
     1 (W(164),H1L),(W(165),EP1L),(W(166),H2L),(W(167),EP2L)
     
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))
      INTEGER CPX    ,CQTBL(10),CLTBL(10),CIRMTBL(10)
      DATA CPX/1/
      DATA CQTBL/1,11,8*0/
      DATA CLTBL/1,9*0/
      DATA CIRMTBL/1,9*0/
C*********************************************************************
      ENTRY SETSPD
C*********************************************************************
      CMX=CPX
      CALL IMOVE(CNTBL,CQTBL,10)
      CALL IMOVE(CITBL,CLTBL,10)
      CALL IMOVE(CFRMTBL,CIRMTBL,10) 
      CALL SETPSP 
      RETURN 
C*********************************************************************
      ENTRY ISPEED
C*********************************************************************
      if (w(150) .ne. 9) then
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 9 if you are using model CSMUNK3.'
      endif

      MODC(1)=7HCSMUNK3
      MODC(2)= w(152)

C    COMPUTE CONSTANT DIFFERENCES
      DPH=PH2-PH1
      DCA=CA2-CA1
      DZA=ZA2-ZA1

C    COMPUTE CONSTANT PH DERIVATIVES OF CA,ZA
      PCAPH=DCA/DPH
      PZAPH=DZA/DPH
      IF(REFC.GT.0.d0) CREF=REFC 
      CALL IPSPEED
      RETURN 
C********************************************************************** 
      ENTRY SPEED
C**********************************************************************
C     INTERPOLATE FOR ZA
      FRACT=(PH-PH1)/DPH
      ZA=ZA1+FRACT*DZA

C     SET H AND EP DEPENDING ON WHETHER RAY IS ABOVE OR BELOW ZA
      IF(R.LT.EARTHR+ZA)THEN
         H1=H1L
         H2=H2L
         EP1=EP1L
         EP2=EP2L
       ELSE
         H1=H1H
         H2=H2H
         EP1=EP1H
         EP2=EP2H
      ENDIF

C     COMPUTE PH DERIVATIVES OF H, EP
      DH=H2-H1
      DEP=EP2-EP1
      PHPH=DH/DPH
      PEPPH=DEP/DPH

C    INTERPOLATE FOR CA,H,EP
      CA=CA1+FRACT*DCA
      H=H1+FRACT*DH
      EP=EP1+FRACT*DEP

C    COMPUTE SOUND SPEED,C
      ETA=-2.d0*(R-EARTHR-ZA)/H
      EXETA1=dEXP(-ETA)-1.d0
      PRODCT=CA*(ETA+EXETA1)
      C=CA+EP*PRODCT
      CS=C*C

C    COMPUTE ALL DERIVATIVES OF C
      PCR=2.d0*CA*EP*EXETA1/H
      PCEP=PRODCT
      PCCA=1.d0+(EP/CA)*PRODCT
      PCZA=-PCR
      PCH=ETA/2.d0*PCR

C    COMPUTE ALL DERIVATIVES OF CS
      C2=C*2.d0
      PCSTH=0.d0
      PCST=0.d0
      PCSR=C2*PCR

C    CHAIN RULE TO GET PCS/PPH.
      PCSPH=C2*(PCH*PHPH+PCEP*PEPPH+PCZA*PZAPH+PCCA*PCAPH)
C
      CALL PSPEED
      RETURN
      END
