
C-----------------------------------------------------------------------------
      SUBROUTINE CSMUNK2  
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 
C     COMMON DECK "CONST" INSERTED HERE   
      implicit double precision (a-h,o-z)
      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),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
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. 6) then
           PRINT*,' Model check number mismatch: '
           STOP ' W(150) should be 6 if you are using model CSMUNK2.'
      endif

      MODC(1)=7HCSMUNK2 
      MODC(2)= w(152)
      IF(REFC.GT.0) CREF=REFC
      DPH=PH2-PH1
      CALL IPSPEED
      RETURN 
C 
      ENTRY SPEED 
C    COMPUTE SOUND SPEED,C1, AT (PH1,Z)  
      ETA1=-2.d0*(R-EARTHR-ZA1)/H1 
      EXETA1=dEXP(-ETA1)-1.d0
      PROD1=CA1*(ETA1+EXETA1)
      C1=CA1+EP1*PROD1       
C    COMPUTE SOUND SPEED C2 AT (PH2,Z)
      ETA2=-2.d0*(R-EARTHR-ZA2)/H2 
      EXETA2=dEXP(-ETA2)-1.d0 
      PROD2=CA2*(ETA2+EXETA2) 
      C2=CA2+EP2*PROD2
C    COMPUTE VERTICAL GRADIENTS AT PH1,PH2
      PCR1=2.d0*CA1*EP1/H1*EXETA1 
      PCR2=2.d0*CA2*EP2/H2*EXETA2 
C    DIFFERENCES IN C AND ITS VERTICAL GRADIENT 
      DPCR=PCR2-PCR1
      DC = C2-C1 
C    INTERPOLATE FOR C AT (PH,Z) 
      FRACT=(PH-PH1)/DPH
      C=C1+FRACT*DC 
      CS=C*C
C    DERIVATIVES OF C AND CS 
      PCPH=DC/DPH
      PCR=PCR1+FRACT*DPCR
      PCSR=2.d0*C*PCR  
      PCSPH=2.d0*C*PCPH
      PCSTH=0.d0  
      PCST=0.d0  
      CALL PSPEED
      RETURN  
      END    


