      PROGRAM IGMV3D
C     A THREE-DIMENSIONAL TIME-DEPENDENT MHD NUMERICAL MODEL
C     BASED ON A MODIFIED LAX-WENDROFF FINITE DIFFERENCE
C     SCHEME . UNITS ARE KM, KG, SECOND AND KELVIN.
C
C     Han, S.M., S.T. Wu, and M. Dryer, A THREE-DIMENSIONAL, TIME-DEPENDENT
C        NUMERICAL MODELING OF SUPER-SONIC, SUPER-ALFVENIC MHD FLOW,
C        COMPUTERS & FLUIDS, Vol. 16, No. 1, pp 81-103, 1988.
C
C     Modified by TR Detman and Z Smith
C     Vector version of time step by TR Detman
C
C     Raw grid parameters:
C     Name   Parameter    units       value at 18 Rs in original version
C     ----   ---------    -----       ----------------------------------
C     Ro      density     Kg/Km^3     2.35E-9
C     Vr      V radial    Km/sec      250.
C     Vt      V theta     Km/sec      Vr*(Bt/Br)   (latitue dependent)
C     Vp      V phi       Km/sec      -Vr*(Bp/Br)  (latitue dependent)
C     Br      B radial    0.1 Gauss   3E-3
C     Bt      B theta     0.1 Gauss   -1.E-5*sin(2*theta)  (latitue dependent)
C     Bp      B phi       0.1 Gauss   -16.E-5*sin(theta)   (latitue dependent)
C     Tp      Temperature Kelvin      1.1E6
C
C     Notes:
C     R, gas constant, units are nT-Km/mole-K
C
C
      PARAMETER ( S2H = 1./3600. )
      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 /BLOC18/ RO0(KX,JX),VR0(KX,JX),VT0(KX,JX),VP0(KX,JX),
     &                BR0(KX,JX),BT0(KX,JX),BP0(KX,JX),TP0(KX,JX)
      COMMON /GSTUFF/ TITLE,IRUN,ISEG,DRX,DTX,DPX
      REAL*8 TCHECK,DTCHK,TSTOP,Tshock
      CHARACTER*80 TITLE
      LOGICAL GO
      DATA GO / .TRUE. /

      WRITE(6,'(/,A)') 'IGMV - VECTOR VERSION '

      CALL INITIALIZE( TCHECK,DTCHK,TSTOP,AVC,CFL)
      J = JCD
      K = KCD

C --- OUTPUT FOR DIV(B) MOVIE, ETC.
      RDT = 1.0
      CALL CFL_DT(RDT,CFL )
      CALL DBMOVIE

      KNT = 0
      WRITE (6,'(A)') ' '
      WRITE (6,'(A,I5,F13.1,F9.1)') 'STEP:', KNT, Time, dT
      CALL DIVCHK
      DO 540 I=1,4
        WRITE (6,2000) RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I),
     &                 BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I)
  540 CONTINUE

      WRITE (6,'(A)') ' '

C --- Begin Time Series output at Earth
      CALL EARTH


C *** TIME LOOP BEGINS FROM HERE. ===================================<*>
   20 CONTINUE

C --- Update inflow boundary condition
       CALL SSMLBC
       CALL SHOCKER( RDT)

C --- FORCE DIV(B) = 0 ON SHELL I=1.5 OR I=2 BY ADJUSTING BR(I=1)
      CALL DIVBH
C      CALL DIVBF

      CALL CFL_DT(RDT,CFL )
      IF ( Time+dT .GT. TSTOP ) dT = TSTOP - Time + 1.0

      KNT = KNT + 1
      WRITE (6,'(A,I5,F13.1,F9.1)') 'STEP:', KNT, Time, dT
      CALL DIVCHK

      CALL T3STEP
      CALL ARVIS(AVC )
      CALL EDGES

C --- WARNING: THIS TEMPERATURE FIX ASSUMED UNIFORM STEADY STATE
      KOUNTB = 0
      KOUNTN = 0
      CPT = TP0(KCD,JCD)/RO0(KCD,JCD)**(GAM-1.)
      TMIN = 1.0
      DO 444 I=2,IC2
	DO 444 J=1,JM
	  DO 444 K=1,KM
	    IF( TP(K,J,I) .LT. 0.0 ) THEN
	      KOUNTN = KOUNTN + 1
	    ENDIF
            CPT = TP0(K,J)/RO0(K,J)**(GAM-1.)
	    TA96 = 0.96*CPT*RO(K,J,I)**(GAM-1.)
	    IF( TP(K,J,I) .LT. TA96 ) THEN
	      KOUNTB = KOUNTB + 1
	      TMIN = AMIN1(TMIN,TP(K,J,I)*0.96/TA96)
C	      TP(K,J,I) = TA96
	    ENDIF
  444 CONTINUE
      ITIME = NINT(Time)
      WRITE (6,*) 'COUNTS:', ITIME, KOUNTB, KOUNTN, TMIN

      Time = Time + dT
C --- TIME STEP COMPLETE

C --- Record Time Series at Earth
      CALL EARTH

C --- OUTPUT ALONG RADIAL FOR MOVIE
      J = JCD
      K = KCD
      ICP = IC2 + 1
      WRITE (8) Time, ICP
      DO 447 I=1,ICP
  447   WRITE (8) RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I),
     &            BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I)

C --- OUTPUT FOR DIV(B) MOVIE
      CALL DBMOVIE

C --- TIME-HISTORY DATA OUTPUTS AT THE  SELECTED POSITIONS.
      J=JCD
      K=KCD
      DO 541 I=1,4
        WRITE (6,2000) RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I),
     &                 BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I)
  541 CONTINUE

      WRITE (6,'(A)') ' '

      IF( Time .LT. TCHECK ) GO TO 550

      CALL CHECK

      TCHECK = TCHECK + DTCHK
  550 CONTINUE

      CALL CONTROL( TSTOP)

      IF( GO .AND. ( Time .LT. TSTOP ) ) GO TO 20
      CALL CHECK

C *** DONE, ----+++*** END OF TIME LOOP ***+++-----------------

      STOP 'Run Completed'
 2000 FORMAT (1P,8E11.3)
      END
C --- ------------------------------------------------------------------
      SUBROUTINE CHECK
C --- WRITE FULL GRID DATA STRUCTURE TO BINARY FILE
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1
      COMMON /BLOC18/ RO0(KX,JX),VR0(KX,JX),VT0(KX,JX),VP0(KX,JX),
     &                BR0(KX,JX),BT0(KX,JX),BP0(KX,JX),TP0(KX,JX)
      CHARACTER TITLE*80, NAME*32, FMT*16
      COMMON /GSTUFF/ TITLE,IRUN,ISEG,DRX,DTX,DPX
      DATA FMT / '(A,I10.10)' /

      ITIME = NINT(Time)

      ISEG = ISEG + 1
      WRITE(NAME,FMT) 'grids/g3b.', ITIME
      OPEN (7,FILE=NAME,FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE (7) TITLE(1:80)
      WRITE (7) ITIME, IRUN, ISEG, IM, JM, KM, DRX, DTX, DPX
      DO 3 J=1,JM
        WRITE (7) ( RO0(K,J),VR0(K,J),VT0(K,J),VP0(K,J),
     &              BR0(K,J),BT0(K,J),BP0(K,J),TP0(K,J),K=1,KM )
    3 CONTINUE
      DO 5 I=1,IM
        DO 5 J=1,JM
          WRITE (7) ( RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I),
     &                BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I),K=1,KM )
    5 CONTINUE

      CLOSE (7)
      WRITE (6,'(2A,I11,F13.1)') 'Wrote Grid: ', NAME, ITIME, Time

      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE CFL_DT(RDT,CFL )
C     THIS ROUTINE DETERMINES A TIME STEP WHICH SATISFIES THE CFL
C     CRITERIA.
C --- Modified 5/6/91, to push the limit.
C
      PARAMETER ( S2H = 1./3600. )
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1
C
      CS2 = 2.*GAM*R*TP(KC1,JC1,IC1)
      CA2 = ( BR(KC1,JC1,IC1)**2+BT(KC1,JC1,IC1)**2+BP(KC1,JC1,IC1)**2 )
     &         /( RO(KC1,JC1,IC1)*XMU0 )
      CF = SQRT( CS2 + CA2 )
      DTMIN = DR/( ABS(VR(KC1,JC1,IC1)) + CF )
      DO 10 I=1,IC2
        DO 10 J=JC1,JC2
          DO 10 K=KC1,KC2
            TEMP = AMAX1( 0.0, TP(K,J,I) )
            CS2 = 2.*GAM*R*TEMP
            CA2 = ( BR(K,J,I)**2 +BT(K,J,I)**2 +BP(K,J,I)**2 )
     &            / ( RO(K,J,I)*XMU0 )
            CF = SQRT( CS2 + CA2 )  
            DTR =         DR/( ABS(VR(K,J,I)) + CF )
            DTT = RAD(I)*DTH/( ABS(VT(K,J,I)) + CF )
            DTP = RAD(I)*DPH/( ABS(VP(K,J,I)) + CF )
            DTMIN=AMIN1(DTMIN,DTR,DTT,DTP)
  10  CONTINUE
      DT=DTMIN*CFL*RDT

C      write(*,*) ' CFL_DT:', DT, DTMIN, RDT, Time

      RETURN
      END
C --- ------------------------------------------------------------------
C     old time step routines were here, new are in file T3step.ftn
C --- ------------------------------------------------------------------
      SUBROUTINE ARVIS(AVC )
C --- Artificial Viscosity Module
C --- Loop through shells smoothing one qty at a time
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1
      REAL RON(JX,KX,2),VRN(JX,KX,2),VTN(JX,KX,2),VPN(JX,KX,2),
     &     BRN(JX,KX,2),BTN(JX,KX,2),BPN(JX,KX,2),TPN(JX,KX,2)
C --- NEW SHELL VECTORIZABLE ARTIFICIAL VISCOSITY............

      ICNT = 1
      DO 450 I=IC1,IC2

          CALL ARVIS1(I,RO,RON(1,1,ICNT),AVC,2,1)

          CALL ARVIS1(I,VR,VRN(1,1,ICNT),AVC,2,2)
          CALL ARVIS1(I,VT,VTN(1,1,ICNT),AVC,2,2)
          CALL ARVIS1(I,VP,VPN(1,1,ICNT),AVC,2,2)

          CALL ARVIS1(I,BR,BRN(1,1,ICNT),AVC,2,1)
          CALL ARVIS1(I,BT,BTN(1,1,ICNT),AVC,1,1)
          CALL ARVIS1(I,BP,BPN(1,1,ICNT),AVC,1,1)

          CALL ARVIS1(I,TP,TPN(1,1,ICNT),AVC,2,3)

        IF( ICNT .EQ. 1 ) THEN
          ICNT = 2
        ELSE
          DO 440 J=JC1,JC2
            DO 440 K=KC1,KC2

              RO(K,J,I-1) = RON(J,K,1)
              VR(K,J,I-1) = VRN(J,K,1)/RON(J,K,1)
              VT(K,J,I-1) = VTN(J,K,1)/RON(J,K,1)
              VP(K,J,I-1) = VPN(J,K,1)/RON(J,K,1)
              BR(K,J,I-1) = BRN(J,K,1)
              BT(K,J,I-1) = BTN(J,K,1)
              BP(K,J,I-1) = BPN(J,K,1)
              TP(K,J,I-1) = TPN(J,K,1)*(GAM-1.)/(2.*R*RON(J,K,1))

              RON(J,K,1) = RON(J,K,2)
              VRN(J,K,1) = VRN(J,K,2)
              VTN(J,K,1) = VTN(J,K,2)
              VPN(J,K,1) = VPN(J,K,2)
              BRN(J,K,1) = BRN(J,K,2)
              BTN(J,K,1) = BTN(J,K,2)
              BPN(J,K,1) = BPN(J,K,2)
              TPN(J,K,1) = TPN(J,K,2)

  440     CONTINUE
        END IF

  450 CONTINUE

      I = IC2
      DO 460 J=JC1,JC2
        DO 460 K=KC1,KC2
          RO(K,J,I) = RON(J,K,1)
          VR(K,J,I) = VRN(J,K,1)/RON(J,K,1)
          VT(K,J,I) = VTN(J,K,1)/RON(J,K,1)
          VP(K,J,I) = VPN(J,K,1)/RON(J,K,1)
          BR(K,J,I) = BRN(J,K,1)
          BT(K,J,I) = BTN(J,K,1)
          BP(K,J,I) = BPN(J,K,1)
          TP(K,J,I) = TPN(J,K,1)*(GAM-1.)/(2.*R*RON(J,K,1))
  460 CONTINUE
C --- -------- END ARTIFICIAL VISCOSITY............
      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE ARVIS1(I,QP,QPN,AC,IEXP,KFLAG)
C
C     THIS ROUTINE PROVIDES ARTIFICIAL DIFFUSION OF LAPIDUS TYPE
C     TO SMOOTH COMPUTATION. A DIFFERENT COEFFICIENT CAN BE USED
C     FOR EACH FLOW VARIABLE.
C     1. CONVERT PRIMITIVE TO FLOW VARIABLE DEPENDING ON IEXP & KFLAG
C     2. SMOOTH IN 3 DIRECTIONS CONSEQUETIVELY
C     3. DO PART OF CONVERSION BACK TO PRIMITIVE: DIVIDE BY RAD(I)**IEXP
C     4. COMPLETE CONVERSION IN MAIN: DIVIDE BY SMOOTHED RO WHERE NEEDED
C
C     Do it for quantity QP on shell I, now called by ARVIS, not main

      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1
      REAL QP(KX,JX,IX), QPN(JX,KX), US(JX,KX), UT(-1:1,JX,KX)
      REAL VRX(-1:1)
      SAVE US, UT
      DATA VX0, XPN / 72644.8, -2.243243 /

C --- COMPUTE FLOW VARIABLE INTO UT FROM PRIMITIVE VARIABLE QP
      DO 20 L=-1,1
        M = I+L
        VRX(L) = VX0*(RAD(M)/RS)**XPN
        DO 20 J = 1,JM
          DO 20 K = 1,KM
            IF (KFLAG.EQ.1) THEN
              UT(L,J,K) = QP(K,J,M)*RAD(M)**IEXP
            ELSE IF (KFLAG.EQ.2) THEN
              UT(L,J,K) = QP(K,J,M)*RO(K,J,M)*RAD(M)**IEXP
            ELSE IF (KFLAG.EQ.3) THEN
              UT(L,J,K) = 2.*RO(K,J,M)*R*QP(K,J,M)/(GAM-1.)*RAD(M)**IEXP
            ENDIF
  20  CONTINUE

C --- R-DIRECTION SWEEP: SMOOTH UT(0,,) INTO US
      TEMP=AC**2*DT/DR
      DO 30 J = 1,JM
        DO 30 K = 1,KM
          DF1 = ABS( (VR(K,J,I+1)+VRX(1)) - (VR(K,J,I)+VRX(0)) )
          DF2 = ABS( (VR(K,J,I)+VRX(0)) - (VR(K,J,I-1)+VRX(-1)) )
          US(J,K) = UT(0,J,K) + TEMP*(DF1*(UT(1,J,K)-UT(0,J,K))
     1                               -DF2*(UT(0,J,K)-UT(-1,J,K)))
   30 CONTINUE

C --- THETA DIRECTION SWEEP: SMOOTH US INTO UT(0,,)
      T1 = RAD(I)
      DO 40 J = JC1,JC2
        TEMP = AC**2*DT/( DTH*T1*SIN( TH(J) ) )
        DO 40 K = KC1,KC2
          DG1 = ABS( VT(K,J+1,I) - VT(K,J,I) )
          DG2 = ABS( VT(K,J,I) - VT(K,J-1,I) )
          UT(0,J,K) = US(J,K) + TEMP*(DG1*(US(J+1,K)-US(J,K))
     1                               -DG2*(US(J,K)-US(J-1,K)))
   40 CONTINUE

C --- PHI-DIRECTION SWEEP: SMOOTH UT(0,,) INTO US
      DO 50 J = JC1,JC2
        TEMP = AC**2*DT/( DPH*T1*SIN( TH(J) ) )
        DO 50 K = KC1,KC2
          DH1 = ABS( VP(K+1,J,I) - VP(K,J,I) )
          DH2 = ABS( VP(K,J,I) - VP(K-1,J,I) )
          US(J,K) = UT(0,J,K) + TEMP*(DH1*(UT(0,J,K+1)-UT(0,J,K))
     1                               -DH2*(UT(0,J,K)-UT(0,J,K-1)))
   50 CONTINUE

C --- UNDO RAD**IEXP PART OF CONVERSION TO FLOW VARIABLES
C --- FOR LATER UPDATE OF GRID AFTER ALL PRIMITIVES ARE SMOOTHED
      DO 60 J = JC1,JC2
        DO 60 K = KC1,KC2
          QPN(J,K) = US(J,K)/RAD(I)**IEXP
   60 CONTINUE

      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE EDGES
C --- THIS ROUTINE PROVIDES BOUNDARY CONDITIONS AT FIVE SIDE BOUNDARIES
C --- BY MIXED ZERO AND FIRST ORDER EXTRAPOLATION AT EACH BOUNDARY

      INCLUDE 'igm_grid.ins'

C --- UPPER BOUNDARY, AT I=IM
      IF( IC2 .LT. IM-1 ) GO TO 15
      I=IM
      TEMP = 0.33
      DO 10 J=1,JM
        DO 10 K=1,KM
          RO(K,J,I) = RO(K,J,I-1) + TEMP*(RO(K,J,I-1) - RO(K,J,I-2))
          VR(K,J,I) = VR(K,J,I-1) + TEMP*(VR(K,J,I-1) - VR(K,J,I-2))
          VT(K,J,I) = VT(K,J,I-1) + TEMP*(VT(K,J,I-1) - VT(K,J,I-2))
          VP(K,J,I) = VP(K,J,I-1) + TEMP*(VP(K,J,I-1) - VP(K,J,I-2))
          BR(K,J,I) = BR(K,J,I-1) + TEMP*(BR(K,J,I-1) - BR(K,J,I-2))
          BT(K,J,I) = BT(K,J,I-1) + TEMP*(BT(K,J,I-1) - BT(K,J,I-2))
          BP(K,J,I) = BP(K,J,I-1) + TEMP*(BP(K,J,I-1) - BP(K,J,I-2))
          TP(K,J,I) = TP(K,J,I-1) + TEMP*(TP(K,J,I-1) - TP(K,J,I-2))
  10  CONTINUE
  15  CONTINUE

C --- NORTH BOUNDARY, AT J=1
      J=1
      TEMP = 0.33
      DO 20 I=2,IM
        DO 20 K=1,KM
          RO(K,J,I) = RO(K,J+1,I) + TEMP*(RO(K,J+1,I) - RO(K,J+2,I))
          VR(K,J,I) = VR(K,J+1,I) + TEMP*(VR(K,J+1,I) - VR(K,J+2,I))
          VT(K,J,I) = VT(K,J+1,I) + TEMP*(VT(K,J+1,I) - VT(K,J+2,I))
          VP(K,J,I) = VP(K,J+1,I) + TEMP*(VP(K,J+1,I) - VP(K,J+2,I))
          BR(K,J,I) = BR(K,J+1,I) + TEMP*(BR(K,J+1,I) - BR(K,J+2,I))
          BT(K,J,I) = BT(K,J+1,I) + TEMP*(BT(K,J+1,I) - BT(K,J+2,I))
          BP(K,J,I) = BP(K,J+1,I) + TEMP*(BP(K,J+1,I) - BP(K,J+2,I))
          TP(K,J,I) = TP(K,J+1,I) + TEMP*(TP(K,J+1,I) - TP(K,J+2,I))
  20  CONTINUE

C --- SOUTH BOUNDARY, AT J=JM
      J=JM
      TEMP = 0.33
      DO 30 I=2,IM
        DO 30 K=1,KM
          RO(K,J,I) = RO(K,J-1,I) + TEMP*(RO(K,J-1,I) - RO(K,J-2,I))
          VR(K,J,I) = VR(K,J-1,I) + TEMP*(VR(K,J-1,I) - VR(K,J-2,I))
          VT(K,J,I) = VT(K,J-1,I) + TEMP*(VT(K,J-1,I) - VT(K,J-2,I))
          VP(K,J,I) = VP(K,J-1,I) + TEMP*(VP(K,J-1,I) - VP(K,J-2,I))
          BR(K,J,I) = BR(K,J-1,I) + TEMP*(BR(K,J-1,I) - BR(K,J-2,I))
          BT(K,J,I) = BT(K,J-1,I) + TEMP*(BT(K,J-1,I) - BT(K,J-2,I))
          BP(K,J,I) = BP(K,J-1,I) + TEMP*(BP(K,J-1,I) - BP(K,J-2,I))
          TP(K,J,I) = TP(K,J-1,I) + TEMP*(TP(K,J-1,I) - TP(K,J-2,I))
  30  CONTINUE

C --- EAST BOUNDARY, AT K=1
      K=1
      TEMP = 0.33
      IF (KM.EQ.3) TEMP=0.
      DO 40 I=2,IM
        DO 40 J=1,JM
          RO(K,J,I) = RO(K+1,J,I) + TEMP*(RO(K+1,J,I) - RO(K+2,J,I))
          VR(K,J,I) = VR(K+1,J,I) + TEMP*(VR(K+1,J,I) - VR(K+2,J,I))
          VT(K,J,I) = VT(K+1,J,I) + TEMP*(VT(K+1,J,I) - VT(K+2,J,I))
          VP(K,J,I) = VP(K+1,J,I) + TEMP*(VP(K+1,J,I) - VP(K+2,J,I))
          BR(K,J,I) = BR(K+1,J,I) + TEMP*(BR(K+1,J,I) - BR(K+2,J,I))
          BT(K,J,I) = BT(K+1,J,I) + TEMP*(BT(K+1,J,I) - BT(K+2,J,I))
          BP(K,J,I) = BP(K+1,J,I) + TEMP*(BP(K+1,J,I) - BP(K+2,J,I))
          TP(K,J,I) = TP(K+1,J,I) + TEMP*(TP(K+1,J,I) - TP(K+2,J,I))
  40  CONTINUE

C --- WEST BOUNDARY, AT K=KM
      K=KM
      TEMP = 0.33
      IF (KM.EQ.3) TEMP=0.
      DO 50 I=2,IM
        DO 50 J=1,JM
          RO(K,J,I) = RO(K-1,J,I) + TEMP*(RO(K-1,J,I) - RO(K-2,J,I))
          VR(K,J,I) = VR(K-1,J,I) + TEMP*(VR(K-1,J,I) - VR(K-2,J,I))
          VT(K,J,I) = VT(K-1,J,I) + TEMP*(VT(K-1,J,I) - VT(K-2,J,I))
          VP(K,J,I) = VP(K-1,J,I) + TEMP*(VP(K-1,J,I) - VP(K-2,J,I))
          BR(K,J,I) = BR(K-1,J,I) + TEMP*(BR(K-1,J,I) - BR(K-2,J,I))
          BT(K,J,I) = BT(K-1,J,I) + TEMP*(BT(K-1,J,I) - BT(K-2,J,I))
          BP(K,J,I) = BP(K-1,J,I) + TEMP*(BP(K-1,J,I) - BP(K-2,J,I))
          TP(K,J,I) = TP(K-1,J,I) + TEMP*(TP(K-1,J,I) - TP(K-2,J,I))
  50  CONTINUE

      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE DBMOVIE
C --- WRITE MOVIE FILE CONTAINING B & DIV(B) RELATED QUANTITIES
C --- TO BE CALLED EACH TIME STEP
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC15/ Tshock, ShTha,ShPhi,ShRad,ShSpd,ShTau, JCD,KCD
      REAL*8 Tshock
      LOGICAL FIRST
      REAL BRL(IX), CRLL(IX)
      SAVE BRL, CRLL, FIRST
      DATA FIRST / .TRUE. /

      IC = IC2
      J = JCD
      K = KCD
      IF( FIRST ) THEN
        OPEN(9,FILE='divBmovi',STATUS='UNKNOWN')
        WRITE(9,*) 'DIV(B) Quantities, grid: ', IM, JM, KM
        IC = IM - 1
        DO 5 I=2,IC
          BRL(I) = BR(K,J,I)
          TEMPJP = BR(K,J+1,I)*VT(K,J+1,I)-BT(K,J+1,I)*VR(K,J+1,I)
          TEMPJP = TEMPJP*SIN(TH(J+1))
          TEMPJM = BR(K,J-1,I)*VT(K,J-1,I)-BT(K,J-1,I)*VR(K,J-1,I)
          TEMPJM = TEMPJM*SIN(TH(J-1))
          CURLJ  = (TEMPJP-TEMPJM)/(2.*DTH*RAD(I)*SIN(TH(J)))
          TEMPKP = (BR(K+1,J,I)*VP(K+1,J,I)-BP(K+1,J,I)*VR(K+1,J,I))  
          TEMPKM = (BR(K-1,J,I)*VP(K-1,J,I)-BP(K-1,J,I)*VR(K-1,J,I))
          CURLK  = (TEMPKP-TEMPKM)/(2.*DPH*RAD(I)*SIN(TH(J)))
          CRLL(I) = CURLJ + CURLK
    5   CONTINUE
        FIRST = .FALSE.
      END IF

      WRITE(9,"(F13.1,I5)") Time, IC-1
      DO 30 I=2,IC
        DIVBMX = -10.0
        DIVBMN =  10.0
        BDIVBM = 0.0
        CRLBBM = 0.0
        R2DVBM = 0.0
        R2CLBM = 0.0
        FFREEM = 1.0E+30
        RM = RAD(I-1)
        RI = RAD(I)
        RP = RAD(I+1)
        RM2 = RM**2
        RI2 = RI**2
        RP2 = RP**2
        DO 20 J=JC1,JC2
          SINTM = SIN(TH(J-1))
          RSINT = RI*SIN(TH(J))
          SINTP = SIN(TH(J+1))
          DO 10 K=KC1,KC2

C ---       COMPUTE DIV(B)
            DBR = ( RP2*BR(K,J,I+1) - RM2*BR(K,J,I-1) )/(RI2*2.*DR)
            DBT = (SINTP*BT(K,J+1,I)-SINTM*BT(K,J-1,I))/(2.*RSINT*DTH)
            DBP = (BP(K+1,J,I) - BP(K-1,J,I))/(2.*RSINT*DPH)
            DIVB = DBR + DBT + DBP
            AB = SQRT( BR(K,J,I)**2 + BT(K,J,I)**2 + BP(K,J,I)**2 )
            BDIVB = AB*DIVB

C ---       COMPUTE CURL(B)
            DSBP = ( SINTP*BP(K,J+1,I) - SINTM*BP(K,J-1,I) )/(2.*DTH)
            DBT = ( BT(K+1,J,I) - BT(K-1,J,I) )/(2.*DPH)
            CBI = ( DSBP - DBT )/RSINT
            DBR = ( BR(K+1,J,I) - BR(K-1,J,I) )/(2.*SIN(TH(J))*DPH)
            DBP = ( RP*BP(K,J,I+1) - RM*BP(K,J,I-1) )/(2.*DR)
            CBJ = ( DBR - DBP )/RI
            DBT = ( RP*BT(K,J,I+1) - RM*BT(K,J,I-1) )/(2.*DR)
            DBR = ( BR(K,J+1,I) - BR(K,J-1,I) )/(2.*DTH)
            CBK = ( DBT - DBR )/RI
            CURLB = SQRT( CBI**2 + CBJ**2 + CBK**2 )

C ---       COMPUTE CURL(B)XB
             D1 = ( BR(K+1,J,I) - BR(K-1,J,I) )/(2.*DPH)
             D2 = ( RP*BP(K,J,I+1) - RM*BP(K,J,I-1) )*SIN(TH(J))/(2.*DR)
             D3 = ( RP*BT(K,J,I+1) - RM*BT(K,J,I-1) )/(2.*DR)
             D4 = ( BR(K,J+1,I) - BR(K,J-1,I) )/(2.*DTH)
             D5 = ( BP(K,J+1,I)*SINTP - BP(K,J-1,I)*SINTM )/(2.*DTH)
             D6 = ( BT(K+1,J,I) - BT(K-1,J,I) )/(2.*DPH)
            CBI = BP(K,J,I)/RSINT*(D1 - D2) - BT(K,J,I)/RI*(D3 - D4)
            CBJ = BR(K,J,I)/RI*(D3 - D4) - BP(K,J,I)/RSINT*(D5 - D6)
            CBK = (BT(K,J,I)*(D5 - D6) - BR(K,J,I)*(D1 - D2))/RSINT
            CRLBXB = SQRT( CBI**2 + CBJ**2 + CBK**2 )

C ---       FIND MAX OF EACH & RATIO OVER COMPUTATIONAL DOMAIN
            DIVBMX = AMAX1(DIVBMX,DIVB)
            DIVBMN = AMIN1(DIVBMN,DIVB)
            BDIVBM = AMAX1(BDIVBM,ABS(BDIVB))
            CRLBBM = AMAX1(CRLBBM,ABS(CRLBXB))
            RATIO  = CRLBXB/(AB*CURLB)
            FFREEM = AMIN1(FFREEM,ABS(RATIO))

            R2DIVB = RI2*BDIVB
            R2CLBB = RI2*CRLBXB
            R2DVBM = AMAX1(R2DVBM,ABS(R2DIVB))
            R2CLBM = AMAX1(R2CLBM,ABS(R2CLBB))

   10     CONTINUE
   20   CONTINUE

C ---   COMPUTE VIOLATION OF FARADAY'S LAW ALONG JCD,KCD RADIAL
        J = JCD
        K = KCD
        TEMPJP = BR(K,J+1,I)*VT(K,J+1,I)-BT(K,J+1,I)*VR(K,J+1,I)
        TEMPJP = TEMPJP*SINTP
        TEMPJM = BR(K,J-1,I)*VT(K,J-1,I)-BT(K,J-1,I)*VR(K,J-1,I)
        TEMPJM = TEMPJM*SINTM
        CURLJ  = (TEMPJP-TEMPJM)/(2.*DTH*RSINT)
        TEMPKP = (BR(K+1,J,K)*VP(K+1,J,K)-BP(K+1,J,K)*VR(K+1,J,K))  
        TEMPKM = (BR(K-1,J,K)*VP(K-1,J,K)-BP(K-1,J,K)*VR(K-1,J,K))
        CURLK  = (TEMPKP-TEMPKM)/(2.*DPH*RSINT)
        CRL    = CURLJ + CURLK
        FDAY   = ( BR(K,J,I) - BRL(I) )/DT + ( CRLL(I) + CRL )/2.
        BRL(I) = BR(K,J,I)
        CRLL(I)= CRL

        WRITE(9,2000) DIVBMN, DIVBMX, BDIVBM, R2DVBM, FDAY,
     $                                CRLBBM, R2CLBM, FFREEM
   30 CONTINUE

      RETURN
 2000 FORMAT (1P,8E12.4)
      END
C --- ------------------------------------------------------------------
      SUBROUTINE DIVBF
C --- FORCE DIVERGENCE B TO ZERO ON SHELL I = 2
C --- BY ADJUSTMENT OF BR ON SHELL I = 1.
      INCLUDE 'igm_grid.ins'
      REAL ERR(KX,JX)

      II0 = 2
      IM1 = II0-1
      IP1 = II0+1
        RM = RAD(IM1)
        RI = RAD(II0)
        RP = RAD(IP1)
        RM2 = RM**2
        RI2 = RI**2
        RP2 = RP**2
        DO 20 J=JC1,JC2
          SINTM = SIN(TH(J-1))
          RSINT = RI*SIN(TH(J))
          SINTP = SIN(TH(J+1))
          DO 10 K=KC1,KC2
	    BTP = BT(K,J+1,II0)
	    BTM = BT(K,J-1,II0)
	    DBT = (SINTP*BTP - SINTM*BTM)/(2.*RSINT*DTH)
	    DBP = (BP(K+1,J,II0) - BP(K-1,J,II0))/(2.*RSINT*DPH)
	    BRKJ1 = (RP2*BR(K,J,IP1) + 2.*RI2*DR*(DBT+DBP))/RM2
	    ERR(K,J) = BRKJ1 - BR(K,J,IM1)
   10     CONTINUE
   20   CONTINUE

C --- DO EDGES, J=1, J=JM, K=1, K=KM BY ZERO/FIRST ORDER EXTRAPOLATION
      I = 1
      TEMP = 0.33
      DO 30 K=1,KM
        ERR(K,JM) = ERR(K,JM-1) + TEMP*(ERR(K,JM-1) - ERR(K,JM-2))
	ERR(K, 1) = ERR(K, 2)   + TEMP*(ERR(K, 2)   - ERR(K, 3))
   30 CONTINUE
      DO 40 J=1,JM
        ERR(KM,J) = ERR(KM-1,J) + TEMP*(ERR(KM-1,J) - ERR(KM-2,J))
        ERR( 1,J) = ERR( 2,J)   + TEMP*( ERR( 2,J)  - ERR( 3,J))
   40 CONTINUE

      DO 70 J=1,JM
       DO 60 K=1,KM
	BR(K,J,IM1) = BR(K,J,IM1) - 0.85*ERR(K,J)
   60  CONTINUE
   70 CONTINUE

C      ITIME = NINT(Time)
C      OPEN (7,FILE='Monopoly.bin',FORM='UNFORMATTED',ACCESS='APPEND')
C      WRITE (7) ITIME, ERR
C      CLOSE (7)

      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE DIVBH
C --- FORCE DIVERGENCE B TO ZERO ON SHELL I = 1 1/2
C --- BY ADJUSTMENT OF BR ON SHELL I = 1.
      INCLUDE 'igm_grid.ins'

      I = 2
      M = 1
        RM = RAD(M)
        RI = RAD(I)
        RC = ( RM + RI )/2.
        RM2 = RM**2
        RI2 = RI**2
        RC2 = RC**2
        DO 20 J=JC1,JC2
          SINTM = SIN(TH(J-1))
          SINTP = SIN(TH(J+1))
          RISINT = RI*SIN(TH(J))
          RMSINT = RM*SIN(TH(J))
          DO 10 K=KC1,KC2
            DBTI = (SINTP*BT(K,J+1,I)-SINTM*BT(K,J-1,I))/(2.*RISINT*DTH)
            DBTM = (SINTP*BT(K,J+1,M)-SINTM*BT(K,J-1,M))/(2.*RMSINT*DTH)
            DBT  = ( DBTI + DBTM )/2.
            DBPI = (BP(K+1,J,I) - BP(K-1,J,I))/(2.*RISINT*DPH)
            DBPM = (BP(K+1,J,M) - BP(K-1,J,M))/(2.*RMSINT*DPH)
            DBP  = ( DBPI + DBPM )/2.
            BR(K,J,M) = RI2/RM2*BR(K,J,I) + RC2/RM2*DR*( DBT + DBP )
   10     CONTINUE
   20   CONTINUE

C --- DO EDGES, J=1, J=JM, K=1, K=KM BY ZERO ORDER EXTRAPOLATION
      I = 1
      DO 30 K=1,KM
        BR(K,JM,I) = BR(K,JM-1,I) 
        BR(K,1,I)  = BR(K,2,I)
   30 CONTINUE
      DO 40 J=1,JM
        BR(KM,J,I) = BR(KM-1,J,I)
        BR(1,J,I)  = BR(2,J,I)
   40 CONTINUE
      
      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE DIVCHK
C --- CHECK DIVERGENCE B.
C --- COMPUTE MAX, MIN, & RMS OF SUM(dB)/SUM|dB|,
C ---         SUM | dR*DIV(B)/B |
C ---         SUM |  R*DIV(B)/B |
C ---         RATIO OF ABOVE
      INCLUDE 'igm_grid.ins'
      COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1

      SUMD = 0.0
      SUMR = 0.0
      KNT  = 0
      RRSUM = 0.0
      FRMAX = 0.0
      FRA = 0.0

      DO 30 I=2,IM-1
        RM = RAD(I-1)
        RI = RAD(I)
        RP = RAD(I+1)
        RM2 = RM**2
        RI2 = RI**2
        RP2 = RP**2
        DO 20 J=JC1,JC2
          SINTM = SIN(TH(J-1))
          RSINT = RI*SIN(TH(J))
          SINTP = SIN(TH(J+1))
          RSINT = RI*SIN(TH(J))
          DO 10 K=KC1,KC2

C ---       COMPUTE DIV(B) MEASURES
	    BTP = BT(K,J+1,I)
	    BTM = BT(K,J-1,I)
            DBR = ( RP2*BR(K,J,I+1) - RM2*BR(K,J,I-1) )/(2.*RI2*DR)
            DBT = (SINTP*BTP-SINTM*BTM)/(2.*RSINT*DTH)
            DBP = (BP(K+1,J,I) - BP(K-1,J,I))/(2.*RSINT*DPH)
            DIVB = DBR + DBT + DBP

	    SUMA = ABS(DBR) + ABS(DBT) + ABS(DBP)
	    RSUM = DIVB/SUMA
	    RRSUM = RRSUM + RSUM**2
            KNT = KNT + 1

            AB = SQRT( BR(K,J,I)**2 + BT(K,J,I)**2 + BP(K,J,I)**2 )
            SUMD = SUMD + ABS( DR*DIVB/AB )
            SUMR = SUMR + ABS( RI*DIVB/AB )

C ---       COMPUTE NONDIMENSIONAL FORCE RATIO W/ GRAVITY
            FR = RI2*ABS( AB*DIVB )/( RO(K,J,I)*GMS*XMU0 )
            FRA = FRA + FR
            FRMAX = AMAX1(FRMAX,FR)

   10     CONTINUE
   20   CONTINUE
   30 CONTINUE

      SUMD = SUMD/KNT
      SUMR = SUMR/KNT
      BARR = SUMR/SUMD
      RRMS = SQRT( RRSUM/KNT )
      FRA  = FRA/KNT

      ITIME = NINT(Time)
      WRITE (6,900) ' DIVCHK: ', ITIME, SUMD, SUMR, BARR, RRMS
     &                                , FRA, FRMAX
  900 FORMAT(A,I11,6G11.4)

      RETURN
      END
C --- ------------------------------------------------------------------
      SUBROUTINE CONTROL( TSTOP)
C --- Special I/O routine to permit outside control of execution
      INCLUDE 'igm_grid.ins'
      PARAMETER ( S2H = 1./3600. )
      CHARACTER SIGNAL*2
      LOGICAL FIRST
      REAL*8 TSTOP, TSTOPX
      SAVE FIRST
      DATA FIRST / .TRUE. /

      ITIME = NINT(Time)
      IF( FIRST ) THEN
        OPEN(3,FILE='control.igm',STATUS='UNKNOWN')
        WRITE(3,9) 'GO', ITIME, dT, TSTOP
        CLOSE(3)
        FIRST = .FALSE.
      END IF

      OPEN(3,FILE='control.igm',STATUS='OLD',ERR=80)
      READ(3,9,ERR=80,END=80) SIGNAL, ITIMEX, dTX, TSTOPX
      REWIND 3
      WRITE(3,9) 'GO', ITIME, dT, TSTOP, (TSTOP-ITIME)/dT
      CLOSE(3)
      TSTOP = TSTOPX
      IF( SIGNAL .NE. 'GO' ) TSTOP = Time
      RETURN

   80 CONTINUE
      TSTOP = Time
      RETURN
    9 FORMAT(A,I11,F8.2,F13.1,F8.1)
      END
