C-----------------------------------------------------------------------------
      SUBROUTINE CSMUNK1  
C-----------------------------------------------------------------------------
C 13 Mar 2008 - Corrected length in common block /PCONST/
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 = 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 
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      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),H1) ,(W(158),EP1)
      EQUIVALENCE (W(159),PH2), (W(160),CA2) ,(W(161),ZA2)
      EQUIVALENCE (W(162),H2) ,(W(163),EP2)
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
      CMX=CPX
      CALL IMOVE(CNTBL,CQTBL,10)
      CALL IMOVE(CITBL,CLTBL,10)
      CALL IMOVE(CFRMTBL,CIRMTBL,10) 
      CALL SETPSP 
      RETURN 
C
      ENTRY ISPEED
      if (w(150) .ne. 5) then
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 5 if you are using model CSMUNK1.'
      endif

      MODC(1)=7HCSMUNK1
      MODC(2)= w(152)
      DPH=PH2-PH1
      DCA=CA2-CA1
      DZA=ZA2-ZA1
      DH=H2-H1   
      DEP=EP2-EP1
C    COMPUTE PH DERIVATIVES OF CA,ZA,H,EP
      PHPH=DH/DPH 
      PEPPH=DEP/DPH
      PCAPH=DCA/DPH
      PZAPH=DZA/DPH
      IF(REFC.GT.0.d0) CREF=REFC 
      CALL IPSPEED
      RETURN 
C 
      ENTRY SPEED 
C    INTERPOLATE FOR CA,ZA,H,EP 
      FRACT=(PH-PH1)/DPH 
      CA=CA1+FRACT*DCA 
      ZA=ZA1+FRACT*DZA
      H=H1+FRACT*DH 
      EP=EP1+FRACT*DEP
C    COMPUTE SOUND SPEED,C
      ETA=-2.d0*(R-EARTHR-ZA)/H
      EXETA1=dEXP(-ETA)-1.
      PRODCT=CA*(ETA+EXETA1)
      C=CA+EP*PRODCT
      CS=C*C
C    COMPUTE 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 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       
