C --- ------------------------------------------------------------------
      SUBROUTINE SHOCKER( RDT)
C ---
C     MHD shock jump conditions are used to introduce boundary
C     perturbation on a portion of the lower boundary at RAD(1)
C     Tshock - POSIX Time at beginning of shock input
C     ShTau  - duration of shock plateau
C     ShLat  - Latitude of shock center on disk, degrees
C     ShLng  - Longitude of shock center on disk, degrees
C     ShRad  - Shock input Radius, input degrees convert to radians
C     SDT  - duration of ramp-up + shock plateau, should be > ramp.
C     SDTE - T at end of shock disturbance, should be > 2.*ramp
C     RAMP - length of linear rise and fall segments of disturbance
C     TFCS - tfc during linear start ramp from t=0 to t=ramp.
C     TFCE - tfc during linear end ramp from t=sdt to t=sdte=std+ramp.
C     RCIR - radius of 'input' region on lower boundary
C ---
C     ShTha - Theta Shock center, radians
C     ShPhi - Phi Shock center, radians
C     JCD   - Theta grid point closest to shock center
C     KCD   - Phi grid point closest to shock center

      PARAMETER ( OMEGA = 2*3.141593/(27.2753*24*3600) )
      PARAMETER ( RAMP = 3600. )
      PARAMETER ( D2R = 3.141593/180. )
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1
      COMMON /BLOC15/ Tshock, ShTha,ShPhi,ShRad,ShSpd,ShTau, JCD,KCD
      COMMON /BLOC16/ SJRo,SJVr,SJVt,SJVp,SJBr,SJBt,SJBp,SJTp
      COMMON/BLOC18/RO0(JX,KX),VR0(JX,KX),VT0(JX,KX),VP0(JX,KX),
     &              BR0(JX,KX),BT0(JX,KX),BP0(JX,KX),TP0(JX,KX)
      REAL*8 Tshock
      CHARACTER*40 FMT1, FMT2
      DATA FMT1 / '(A,F12.0,F9.1,F12.0,3F9.3)' /
      DATA FMT2 / '(A,F8.1,2F8.4,2F8.3)' /
      LOGICAL FIRST
      DATA FIRST / .TRUE. /
      SAVE FIRST

      IF ( Time .LT. Tshock ) THEN     ! Not yet
        RDT = 1.0
        RETURN
      ENDIF

      T = Time - Tshock
      SDT = ShTau + RAMP
      SDTE = SDT + RAMP
      IF( T .GT. SDTE+3600. ) THEN     ! shock is over
        RDT = 1.0
        RETURN
      ENDIF

      IF ( .NOT. FIRST ) GO TO 40
C --- Compute shock jumps; use current Vr and Tp near shock center
      VrBk = Vr(KCD,JCD,1)     ! Background Speed
      TpBk = Tp(KCD,JCD,1)     ! Background Temperature
      RR = 10.*0.00826527*TpBk/(ShSpd - VrBk)**2
      RR = AMAX1(4./(RR+1.),1.)

      SJRo = RR
      SJVr = (ShSpd*(RR-1.) + VrBk)/(RR*VrBk)
      SJVt = 1.
      SJVp = 1.
      SJBr = 1.
      SJBt = 1.
      SJBp = 1.
      SJTp = (4.*RR - 1.)/(RR*(4. - RR))
      write(6,FMT1) '  Start shock:', time, VrBk,TpBk, SJRo,SJVr,SJTp
      FIRST = .FALSE.

   40 CONTINUE
      R1 = RAD(1)
      IF( T .LT. RAMP ) TFC = T/RAMP
      IF( T .GE. RAMP .AND. T .LE. SDT ) TFC = 1.0
      IF( T .GT. SDT .AND. T .LT. SDTE ) TFC = ( SDTE - T )/RAMP
      IF( T .GE. SDTE ) TFC = 0.0

      PhiC = ShPhi + OMEGA*T
      RCIR = R1*SIN(ShTha)*ShRad
      DO 5 J=1,JM
         DO 5 K=1,KM
            A1 = R1*SIN(ShTha)*(PHI(K)-PhiC)
            A2 = R1*(ShTha-TH(J))
            DSTCR = SQRT( A1**2 + A2**2 )
            A = DSTCR/RCIR
            IF ( A .GE. 1.0 ) A = 1.0
            CFSP = COS(A*PI/2.)
            Frac = TFC*CFSP

	    FRo       = (1.+(SJRo-1.)*Frac)
            Ro(K,J,1) = Ro(K,J,1)*FRo
	    FVr       = (1.+(SJVr-1.)*Frac)
            Vr(K,J,1) = Vr(K,J,1)*FVr
            Vt(K,J,1) = Vt(K,J,1)*(1.+(SJVt-1.)*Frac)
            Vp(K,J,1) = Vp(K,J,1)*(1.+(SJVp-1.)*Frac)
            Br(K,J,1) = Br(K,J,1)*(1.+(SJBr-1.)*Frac)
            Bt(K,J,1) = Bt(K,J,1)*(1.+(SJBt-1.)*Frac)
            Bp(K,J,1) = Bp(K,J,1)*(1.+(SJBp-1.)*Frac)
            Tp(K,J,1) = Tp(K,J,1)*(1.+(SJTp-1.)*Frac)
	    if ( j .eq. JCD .and. k .eq. KCD )
     &         write(6,FMT2) '  Shocker:', T, CFSP, TFC, FRo,FVr
    5 CONTINUE

C --- Mod to control DT during up ramp and plateau
      IF( T .LT. SDT ) THEN
        K = KCD
        J = JCD
        CMCS2 = 2.*GAM*R*RO(K,J,1)*TP(K,J,1)
        SHCS2 = 2.*GAM*R*RO0(J,K)*TP0(J,K)*SJRO*SJTP
        CMB2 = BR(K,J,1)**2 + BT(K,J,1)**2 + BP(K,J,1)**2
        CMCA2 = CMB2/( RO(K,J,1)*XMU0 )
        SHB2 = (SJBR*BR0(J,K))**2+(SJBT*BT0(J,K))**2+(SJBP*BP0(J,K))**2
        SHCA2 = SHB2/( SJRO*RO0(J,K)*XMU0 )
        CMCF = SQRT( CMCS2 + CMCA2 )
        SHCF = SQRT( SHCS2 + SHCA2 )
        CMDT = DR/( ABS(VR(K,J,1)) + CMCF )
        SHDT = DR/( ABS(SJVR*VR0(J,K)) + SHCF )
        RDT = AMIN1( SHDT/CMDT, 1.0 )
      ELSE
        RDT = 1.0
      ENDIF

      RETURN
      END