C-----------------------------------------------------------------------------
      SUBROUTINE CBLOB2
C-----------------------------------------------------------------------------
C     SOUND SPEED PERTURBATION MODEL
C     MULTIPLICATIVE PERTURBATION WITH EXPONENTIALLY DECAYING
C     EFFECT IN ALL THREE DIRECTIONS. GIVE LATITUDE 
C     INSTEAD OF CO-LATITUDE.
      implicit double precision (a-h,o-z)
      COMMON/PCONST/CREF,RGAS,GAMMA 
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
      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(20)
      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(10)
      EQUIVALENCE (DCGP,IDSDC) 
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     DELTA SOUND SPEED 175-199
      EQUIVALENCE (C0,W(178)), (Z0,W(179)),(CBTH0,W(180))
      EQUIVALENCE (PH0,W(181)),(WZ,W(182)),(WTH,W(183)),(WPH,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
      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. 2.0) then
           PRINT*,' Model check number mismatch: '
           STOP   ' W(175) should be 2 if you are using model CBLOB2.'
      endif

      MODC(3)=8HCBLOB2  
      MODC(4)= w(177)  ! ID number for a particular model's parameters
      FWZ=0.d0
      FWTH=0.d0
      FWPH=0.d0
      TH0= PID2-CBTH0
      IF(WZ.NE.0.d0) FWZ=2.d0/WZ/WZ
      IF(WTH.NE.0.d0) FWTH=2.d0/WTH/WTH
      IF(WPH.NE.0.d0) FWPH=2.d0/WPH/WPH
      RETURN
C
      ENTRY PSPEED
      IF(C0.EQ.0.0) RETURN
      DZ=R-EARTHR-Z0
      DTH=TH-TH0
      DPH=PH-PH0
      DEXPO=0.0 
      EXPO=-0.5*(DZ*DZ*FWZ+DTH*DTH*FWTH+DPH*DPH*FWPH)
      IF(EXPO .GT. -200.d0) DEXPO=C0*dEXP(EXPO)
      DEL=1.d0+DEXPO 
      PCSR=PCSR*DEL-CS*DEXPO*FWZ*DZ 
      PCSTH=PCSTH*DEL-CS*DEXPO*FWTH*DTH
      PCSPH=PCSPH*DEL-CS*DEXPO*FWPH*DPH
      CS=CS*DEL
      RETURN
      END 
    
 
