=================== beginning of mhd2.f ================================
C      FILE NAME  /JSSS/RECTG81
C      COMPUTER SIMULATION OF THE INTERACTION BETWEEN THE SOLAR
C      WIND AND THE EARTH'S MAGNETOSPHERE BY USING TWO
C      DIMENSIONAL GLOBAL MAGNETOHYDRODYNAMIC (MHD) MODEL
C      IN THE TWO STEP LAX-WENDROFF METHOD
C      CARTESIAN COORDINATE   FINITE RESISTIVITY
C      1986.05.16  BY TATSUKI OGINO
       PARAMETER (NR=40,NZ= 20,N1=NR+2,N2=N1*(NZ+2))
       PARAMETER (NB=6,NBB=7,N3=N2*NB,N4=N2*NBB,N5=N2*8)
       PARAMETER (LAST=128,LLL=6,IBIG=0,IICH=  32)
       PARAMETER (IIQ0=32,IIP0= 32,IIS0= 32,THR=16.0,TAM=0.50)
       PARAMETER (RMU0=0.005,EUD0=0.00,RO01=5.0E-4,PR01=3.56E-8)
       PARAMETER (NRP=  0,ARU=100.0,PMU0=1.0,DFU0=1.0)
       PARAMETER (EAT0=0.010,RAT1=0.2,RATT=0.2)
       PARAMETER (NIN=0,ITAP=9)
       PARAMETER (K2=N1*2,K3=N1*3,K4=N1*4,K5=N1*5,K6=N1*6,K7=N1*7)
       PARAMETER (KK2=K3*2,KK3=K3*3,KK4=K3*4,KK5=K3*5,KK6=K3*6)
       PARAMETER (N11=N1+1,NR1=NR+1,NR2=NR+2,NZ1=NZ+1,NZ2=NZ+2)
       PARAMETER (M2=N2*2,M3=N2*3,M4=N2*4,M5=N2*5)
       PARAMETER (NRZ=NR*NZ,MN1=N1*NZ1)
       DIMENSION F(N3),FF(KK6),U(KK6),P(K7),ZMAX(LLL),ZMIN(LLL),
     1           PP(N2),CJ(10),CP(10),FBB(6),HHX(N1)
C      M  1-RO  2-VR  3-VZ  4-P  5-BR  6-BZ  7-JY  8-JYZ
C
C      CJ      R1   DR  DZ  B0  ALP  BET   AID AR1 GAM GRA
       DATA CJ/10.0,0.3,0.5,1.0,10.0, 1.20,0.7,0.8,0.9,1.0/
C
C      CP      R0   Z0   RA  VO0  P0   GRA  07  VSW  IBR IBZ
C      DATA CP/600.,150.,16.,6.81,2.68,1.35,7.0,.044,9.0,0.0/
       DATA CP/202.,102.,16.,0.00,2.68,1.35,7.0,.044,0.0,-1.5/
C
C      F(N3)   MAIN VARIABLES OF 2-DIMENSIONAL MHD SIMULATION
C      FF(KK6) KEEP THE ORIGINAL QUANTITIES IN FIRST STEP
C      U(KK6)  VARIABLES IN FIRTST STEP
C      P(K7)   VARIABLES IN INTERMEDIATE STAGE
C      PP(N2)  CONSTANT MODIFICATION CURRENT FOR DIPOLE FIELD
C      CJ(10)  VARIABLES TO DETERMINE INITIAL CONDITION
C      CP(10)  VARIABLES TO DETERMINE INITIAL CONDITION
C      FBB(6)  COEFFICIENT TO DEFINE BOUNDARY CONDITION
C      HHX(N1) DISTANCE IN ARRAY
C      NR     GRID NUMBER IN X-DIRECTION WITHOUT BOUNDARY
C      NZ     GRID NUMBER IN Z-DIRECTION WITHOUT BOUNDARY
C      NB     COMPONENTS OF VECTOR (RO,VX,VZ,P,BX,BZ)
C      NBB    COMPONENTS OF VECTOR (RO,VX,VZ,P,BX,BZ,JY)
C      LAST   LAST STEP OF SIMULATION
C      IIQ0   SAMPLING INTERVAL
C      IIP0   SAMPLING INTERVAL
C      IIS0   SAMPLING INTERVAL
C      THR    A MAGUNIFICATION COEFFICIENT FOR TIME STEP
C      TAM    CONTROL FACTOR FOR TIME INTEGRATION
C      RMU0   COEFFICIENT OF VISCOSITY
C      FUD0   COEFFICIENT OF FRICTION
C      RO01   DENSITY OF SOLAR WIND
C      PR01   PLASMA PRESSURE OF SOLAR WIND
C      PMU0   FACTOR FOR PLASMA PRESSURE DIFFUSION
C      DFU0   FACTOR FOR PLASMA DENSITY DIFFUSION
C      EAT0   COEFFICIENT OF RESISTIVITY
C      RAT1   FACTOR TO DETERMINE MINIMUM PLASMA DENSITY
C      RRTT   FACTOR TO DETERMINE MINIMUM PLASMA DENSITY
C      NIN    CONTROL NUMBER FOR FIRST RUN OR FOLLOWING RUN
C      ITAP   CONTROL NUMBER TO SKIP OLD FILES
C      CJ(1)  R1
C      CJ(2)  DR  RATT  FACTOR TO DETERMINE MINIMUM PLASMA DENSITY
C      CJ(3)  DZ
C      CJ(4)  BO  MAGNITUDE OF DIPOLE MAGNETIC FIELD
C      CJ(5)  ALP FACTOR TO DETERMINE INITIAL SOLAR WIND POSITION
C      CJ(6)  BET FACTOR TO DETERMINE INITIAL SOLAR WIND POSITION
C      CJ(7)  AID
C      CJ(8)  AR1 RADIUS OF INNER BOUNDARY
C      CJ(9)  GAM RATIO OF SPECIFIC HEAT
C      CJ(10) GRA GRAVITATIONAL ACCELARETION
C      CP(1)  RO  DISTANCE IN X-DIRECTION
C      CP(2)  ZO  DISTANCE IN Z-DIRECTION
C      CP(3)  RA  DISTANCE OF IONOSPHERIC INNER BOUNDARY
C      CP(4)  VO0 ---- ROTATIONAL VELOCITY
C      CP(5)  PO  PLASMA PRESSURE AT IONOSPHERE
C      CP(6)  GRA GRAVITATIONAL ACCELERATION
C      CP(7)  07  ----
C      CP(8)  VSW VELOCITY OF SOLAR WIND
C      CP(9)  IBR ---- IMF IN X-DIRECTION
C      CP(10) IBZ IMF IN Z-DIRECTION
C      GAM    RATIO OF SPECIFIC HEAT
C      GRA    GRAVITATIONAL ACCELERATION
C      PO0    PLASMA PRESSURE OF IONOSPHERE
C      BIS    Z-COMPONENT OF INTERPLANETARY MAGNETIC FIELD
C      VSW    SOLAR WIND VELOCITY
C      AR1    CRITICAL RADIUS TO INNER IONOSPHERE
C      HR     GRID INTERVAL IN X-DIRECTION
C      HZ     GRID INTERVAL IN Z-DIRECTION
C      T      TIME STEP OF SECOND STEP OF SIMULATION
C      T1     TIME STEP OF FIRST STEP
C      RMU    VISICOUS COEFFICIENT OF VELOCITY
C      PMU    DIFFUSION COEFFICIENT OF PLASMA PRESSURE
C      DFU    DIFFUSION COEFFICIENT OF PLASMA DENSITY
C      EUD    FRICTION COEFFICIENT OF VELOCITY
C      RO02   MINIMUM PLASMA DENSITY
C      PR02   MINIMUM PLASMA PRESSURE
C
       REWIND 12
C
C      INTRODUCTION OF UPSTREAM UNIFORM IMF BZ COMPONENT
C      FOR NORTHWARD, ZERO AND SOUTHWARD CASES
C
       DO 300 II1=1,3
       CP(10)=1.5*FLOAT(2-II1)
C
C      DETERMINATION OF INITIAL PARAMETERS IN MHD MODEL
       PI=3.1415926
       GAM=2.0
       GRA=CP(6)*1.0E-6
       PO0=CP(5)*1.0E-7
       BIS=CP(10)*1.0E-4
C      CP(10)=0.0
C      CP(10)=CP(10)
       CJ(2)=RATT
       VSW=CP(8)
       AR1=CP(3)
       AR2=AR1*AR1
       HR=CP(1)/FLOAT(NR1)
       HZ=CP(2)/FLOAT(NZ1)
       T=0.5*HR*THR
       T1=TAM*T
       CJ(8)=AR1
       CJ(9)=GAM
       CJ(10)=GRA
C
       DR1=0.5*T1/HR
       DZ1=0.5*T1/HZ
       DR2=0.5*T/HR
       DZ2=0.5*T/HZ
       DR3=T1/(HR*HR)
       DZ3=T1/(HZ*HZ)
       DR4=T/(HR*HR)
       DZ4=T/(HZ*HZ)
C
       RMU=RMU0*RO01
       PMU=RMU0*PMU0
       DFU=RMU0*DFU0
       EUD=EUD0
       RO02=RAT1*RO01
       PR02=RAT1*PR01
C
       WRITE (6,12) LAST,NR,NZ,N1,N2,N3,N4,RMU,EUD,
     1 HR,HZ,T,T1,RO01,PR01,GRA,DR2,DZ2,DR4,DZ4,EAT0,
     2 (CP(I),I=1,10),(CJ(J),J=1,10)
   12  FORMAT(1H ,5X,7I10/(1H ,5X,1P8E15.5))
C
       DO 20 I=1,N3
   20  F(I)=0.0
       DO 221 I=1,N1
  221  HHX(I)=HR*FLOAT(I)
       DO 222 I=1,NB
  222  FBB(I)=1.0
       FBB(3)=-1.0
       FBB(5)=-1.0
C
C      INITIAL CONDITION BEGINNING
c      CALL EQUIP4(NR,NZ,NRP,RO01,PR01,CJ,CP,F,PP)
       CALL EQUIP5(NR,NZ,NRP,RO01,PR01,CJ,CP,F,PP)
       GO TO 882
       WRITE(6,212) (F(I),I=1,N3)
  212  FORMAT(1H ,2X,1P10E12.4)
  882  CONTINUE
C
       JXG=N2
       JYG=NZ2
       JXN=NR2
       LANK=20
       IP=0
       LB=1
       LL=6
       I1=II
C      CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F)
       DO 130 J=1,2
       LB=1+3*(J-1)
       LL=3
C      CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F)
  130  CONTINUE
C
C      INITIAL CONDITION END
C
       IF(ITAP.EQ.999) GO TO 400
       DO 410 I=1,ITAP
C      READ(13) F
  410  CONTINUE
  400  CONTINUE
C      TWO DIMENSIONAL CARTESIAN MODEL
       TIME=0.0
       IIQ=0
       IIP=0
       IIS=0
C
C
C      MAIN TEMPORAL ITARATION (ROOP DO 100) BEGINS
C      IN THE TWO STEP LAX-WENDROFF METHOD
C
       DO 100 II=1,LAST
       IIQ=IIQ+1
       IIP=IIP+1
       IIS=IIS+1
C
C      BOUNDARY CONDITION BEGINS
C      BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2
       X1=0.5*HR*FLOAT(NR1-2+2*NRP)
       DO 30 M=1,NB
       I3=N1*2
       L1=N2*(M-1)
       X=(HR*FLOAT(NR1)-HHX(2))/X1
       X=AMIN1(X,1.0)
       X=1.0+X/3.0
       XX=1.0-X
       I1=2+L1
       I2=I1+MN1
       P(2)=X*F(I2-N1)+XX*F(I2-I3)
       U(2)=F(I1+N1)*FBB(M)
       DO 31 I=3,NR1
       X=(HR*FLOAT(NR1)-HHX(I))/X1
       X=AMIN1(X,1.0)
       X=1.0+X/3.0
       XX=1.0-X
       I1=I+L1
       I2=I1+MN1
       P(I)=X*F(I2-N1-1)+XX*F(I2-I3-2)
   31  U(I)=F(I1+N1)*FBB(M)
       DO 32 I=2,NR1
       I2=I+L1+MN1
       F(I2)=P(I)
   32  F(I2-MN1)=U(I)
C      BOUNDARY CONDITION AT NR=1
C      FREE BOUNDARY AT NR=NR2
       I1=N2*(M-1)
       DO 33 J=1,NZ2
       I1=I1+N1
   33  P(J)=F(I1-1)
       I1=N2*(M-1)
       DO 30 J=1,NZ2
       I1=I1+N1
       F(I1)=P(J)
   30  CONTINUE
       DO 34 I=1,N2
       F(I+M3)=AMAX1(F(I+M3),PR02)
       F(I)=AMAX1(F(I),RO02)
   34  CONTINUE
C
       IF(II.NE.1) GO TO 136
       IF(NIN.EQ.0) WRITE(12) F
       JXG=N2
       JYG=NZ2
       JXN=NR2
       LANK=20
       IP=0
       LB=1
       LL=6
       I1=II
C      CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F)
       DO 134 J=1,2
       LB=1+3*(J-1)
       LL=3
C      CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F)
  134  CONTINUE
  136  CONTINUE
C
C      BOUNDARY CONDITION ENDS
C
C      J=1 STEP CALCULATION BEGINS
C      STEP OF J=1
       J=1
C
C      CALCULATION OF CURRENT FROM MAGNETIC FIELD
C      CURRENT
       L1=N1*(J-1)
       DO 381 I=1,NR1
       I1=I+L1
       I5=I1+M4
       I6=I1+M5
       J7=I+K6
       P(J7)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ
     1        -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR)
     2        -PP(I1)
  381  CONTINUE
C
C      INTERPOLATION OF P(J1) AND U(J2)
       DO 382 M=1,NB
       L4=L1+N2*(M-1)
       L2=N1*(M-1)
       L3=K3*(M-1)
       DO 382 I=1,NR1
       I1=I+L4
       J1=I+L2
       J2=I+L3
       P(J1)=0.25*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11))
       U(J2)=P(J1)
  382  CONTINUE
C
C      ZEROTH STEP CALCULATION BEGINS
C      ZEROTH STEP
       XX4=0.5*HR*FLOAT(2*NRP-NR1-1)
       Z=0.5*HZ*FLOAT(2*J-2)
       DO 501 I=1,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       RO2=AMAX1(RO2,AR2)
       I1=I+L1
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       J5=J4+N1
       J6=J5+N1
       J7=J6+N1
       JJ1=I
       X2=RO2/AR2-1.0
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C
       U(JJ1+KK5)=U(JJ1+KK5)
     1       +DR1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11)
     2       +F(I3+1)*F(I5+1)-F(I2+1)*F(I6+1)
     3       -F(I3+N1)*F(I5+N1)+F(I2+N1)*F(I6+N1)
     4       -F(I3)*F(I5)+F(I2)*F(I6))
       U(JJ1+KK4)=U(JJ1+KK4)
     1       -DZ1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11)
     2       -F(I3+1)*F(I5+1)+F(I2+1)*F(I6+1)
     3       +F(I3+N1)*F(I5+N1)-F(I2+N1)*F(I6+N1)
     4       -F(I3)*F(I5)+F(I2)*F(I6))
       X3=DR1*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2))
     1    +DZ1*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3))
       U(JJ1+KK3)=U(JJ1+KK3)-P(J4)*GAM*X3
     1       -DR1*P(J2)*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))
     2       -DZ1*P(J3)*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))
       U(JJ1+KK2)=U(JJ1+KK2)+T1*((-P(J7)*P(J5))/P(J1)
     1       -EUD*P(J3)-GRA*Z/RO2)
     2       -DR1*P(J2)*(F(I3+N11)+F(I3+1)-F(I3+N1)-F(I3))
     3       -DZ1*P(J3)*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3))
     4       -DZ1*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))/P(J1)
       U(JJ1+K3)=U(JJ1+K3)+T1*((P(J7)*P(J6))/P(J1)
     1       -EUD*P(J2)-GRA*R/RO2)
     2       -DR1*P(J2)*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2))
     3       -DZ1*P(J3)*(F(I2+N11)-F(I2+1)+F(I2+N1)-F(I2))
     4       -DR1*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))/P(J1)
       U(JJ1)=U(JJ1)
     1       -DR1*(F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1)
     2       -F(I2+N1)*F(I1+N1)-F(I2)*F(I1))
     3       -DZ1*(F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1)
     4       +F(I3+N1)*F(I1+N1)-F(I3)*F(I1))
  501  CONTINUE
C
C      SMOOTH CONNECTION WITH INNER BOUNDARY
       DO 502 I=1,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       J5=J4+N1
       J6=J5+N1
       J7=J6+N1
       JJ1=I
       X2=RO2/AR2-1.0
       X2=AMAX1(X2,0.0)
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C
       U(JJ1+KK5)=U(JJ1+KK5)*X2+P(J6)*(1.0-X2)
       U(JJ1+KK4)=U(JJ1+KK4)*X2+P(J5)*(1.0-X2)
       U(JJ1+KK3)=U(JJ1+KK3)*X2+P(J4)*(1.0-X2)
       U(JJ1+KK2)=U(JJ1+KK2)*X2+P(J3)*(1.0-X2)
       U(JJ1+K3)=U(JJ1+K3)*X2+P(J2)*(1.0-X2)
       U(JJ1)=U(JJ1)*X2+P(J1)*(1.0-X2)
       U(JJ1)=AMAX1(U(JJ1),RO02)
  502  CONTINUE
C
C      PRESERVE F(I1) WITH FF(J1)
       DO 621 M=1,NB
       L4=L1+N2*(M-1)
       L3=K3*(M-1)
       DO 621 I=1,N1
       I1=I+L4
       J1=I+L3
       FF(J1)=F(I1)
  621  CONTINUE
C
C      J=J STEP CALCULATION BEGINS
C      STEP OF J=J
       DO 90 J=2,NZ1
C
C      CALCULATION OF CURRENT FROM MAGNETIC FIELD
C      CURRENT
       L1=N1*(J-1)
       DO 38 I=1,NR1
       I1=I+L1
       I5=I1+M4
       I6=I1+M5
       J7=I+K6
       P(J7)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ
     1        -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR)
     2        -PP(I1)
   38  CONTINUE
C
C      INTERPOLATION OF P(J1) AND U(J2)
       DO 39 M=1,NB
       L4=L1+N2*(M-1)
       L2=N1*(M-1)
       L3=N1+K3*(M-1)
       DO 39 I=1,NR1
       I1=I+L4
       J1=I+L2
       J2=I+L3
       P(J1)=0.25*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11))
       U(J2)=P(J1)
   39  CONTINUE
C
C      FIRST STEP CALCULATION BEGINS
C      FIRST STEP
       XX4=0.5*HR*FLOAT(2*NRP-NR1-1)
       Z=0.5*HZ*FLOAT(2*J-2)
       DO 50 I=1,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       RO2=AMAX1(RO2,AR2)
       I1=I+L1
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       J5=J4+N1
       J6=J5+N1
       J7=J6+N1
       JJ1=I+N1
       X2=RO2/AR2-1.0
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C
       U(JJ1+KK5)=U(JJ1+KK5)
     1       +DR1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11)
     2       +F(I3+1)*F(I5+1)-F(I2+1)*F(I6+1)
     3       -F(I3+N1)*F(I5+N1)+F(I2+N1)*F(I6+N1)
     4       -F(I3)*F(I5)+F(I2)*F(I6))
       U(JJ1+KK4)=U(JJ1+KK4)
     1       -DZ1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11)
     2       -F(I3+1)*F(I5+1)+F(I2+1)*F(I6+1)
     3       +F(I3+N1)*F(I5+N1)-F(I2+N1)*F(I6+N1)
     4       -F(I3)*F(I5)+F(I2)*F(I6))
       X3=DR1*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2))
     1    +DZ1*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3))
       U(JJ1+KK3)=U(JJ1+KK3)-P(J4)*GAM*X3
     1       -DR1*P(J2)*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))
     2       -DZ1*P(J3)*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))
       U(JJ1+KK2)=U(JJ1+KK2)+T1*((-P(J7)*P(J5))/P(J1)
     1       -EUD*P(J3)-GRA*Z/RO2)
     2       -DR1*P(J2)*(F(I3+N11)+F(I3+1)-F(I3+N1)-F(I3))
     3       -DZ1*P(J3)*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3))
     4       -DZ1*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))/P(J1)
       U(JJ1+K3)=U(JJ1+K3)+T1*((P(J7)*P(J6))/P(J1)
     1       -EUD*P(J2)-GRA*R/RO2)
     2       -DR1*P(J2)*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2))
     3       -DZ1*P(J3)*(F(I2+N11)-F(I2+1)+F(I2+N1)-F(I2))
     4       -DR1*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))/P(J1)
       U(JJ1)=U(JJ1)
     1       -DR1*(F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1)
     2       -F(I2+N1)*F(I1+N1)-F(I2)*F(I1))
     3       -DZ1*(F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1)
     4       +F(I3+N1)*F(I1+N1)-F(I3)*F(I1))
   50  CONTINUE
C
C      SMOOTH CONNECTION WITH INNER BOUNDARY
       DO 51 I=1,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       J5=J4+N1
       J6=J5+N1
       J7=J6+N1
       JJ1=I+N1
       X2=RO2/AR2-1.0
       X2=AMAX1(X2,0.0)
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C
       U(JJ1+KK5)=U(JJ1+KK5)*X2+P(J6)*(1.0-X2)
       U(JJ1+KK4)=U(JJ1+KK4)*X2+P(J5)*(1.0-X2)
       U(JJ1+KK3)=U(JJ1+KK3)*X2+P(J4)*(1.0-X2)
       U(JJ1+KK2)=U(JJ1+KK2)*X2+P(J3)*(1.0-X2)
       U(JJ1+K3)=U(JJ1+K3)*X2+P(J2)*(1.0-X2)
       U(JJ1)=U(JJ1)*X2+P(J1)*(1.0-X2)
       U(JJ1)=AMAX1(U(JJ1),RO02)
   51  CONTINUE
C
C
C      SECOND STEP CALCULATION BEBINING
C      SECOND STEP
C
C      PRESERVE F(I1) WITH FF(J1)
       DO 62 M=1,NB
       L4=L1+N2*(M-1)
       L3=N1+K3*(M-1)
       DO 62 I=1,N1
       I1=I+L4
       J1=I+L3
       FF(J1+N1)=F(I1+N1)
       FF(J1)=F(I1)
   62  CONTINUE
C
C      CURRENT CALCULATION FROM MAGNETIC FIELD
C      CURRENT
       L1=N1*(J-1)
       DO 68 I=2,NR1
       I1=I+L1
       I5=I+N1+KK4
       I6=I5+K3
       J7=I+K6
       P(J7)=0.5*((U(I5)-U(I5-N1)+U(I5-1)-U(I5-N11))/HZ
     1        -(U(I6)+U(I6-N1)-U(I6-1)-U(I6-N11))/HR)
     2        -0.25*(PP(I1)+PP(I1-1)+PP(I1-N1)+PP(I1-N11))
   68  CONTINUE
C
C      INTERPOLATION OF P(J1)
       DO 69 M=1,NB
       L2=N1+K3*(M-1)
       L3=N1*(M-1)
       DO 69 I=2,NR1
       I1=I+L2
       J1=I+L3
       P(J1)=0.25*(U(I1)+U(I1-1)+U(I1-N1)+U(I1-N11))
   69  CONTINUE
C
C      SECOND STEP CALCULATION BEGINS
C      DIFFUSION CALCULATION IS DONE IN SECOND STEP ONLY
C      SECOND STEP
       XX4=0.5*HR*FLOAT(2*NRP-NR1-2)
       Z=0.5*HZ*FLOAT(2*J-3)
       DO 80 I=2,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       RO2=AMAX1(RO2,AR2)
       I1=I+N1
       I2=I1+K3
       I3=I2+K3
       I4=I3+K3
       I5=I4+K3
       I6=I5+K3
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       J5=J4+N1
       J6=J5+N1
       J7=J6+N1
       JJ1=I+L1
       X1=FF(I4)/(FF(I1)*PO0)
       X2=RO2/AR2-1.0
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C      CLASSICAL RESISTIVITY  EAT=EAT0*X2/SQRT(X1*X1*X1)
       EAT=EAT0*X2/SQRT(X1*X1*X1)
C      CONSTANT RESISTIVITY   EAT=EAT0*X2
C      EAT=EAT0*X2
C
       F(JJ1+M5)=F(JJ1+M5)
     1       +DR2*(U(I3)*U(I5)-U(I2)*U(I6)
     2       +U(I3-N1)*U(I5-N1)-U(I2-N1)*U(I6-N1)
     3       -U(I3-1)*U(I5-1)+U(I2-1)*U(I6-1)
     4       -U(I3-N11)*U(I5-N11)+U(I2-N11)*U(I6-N11))
     5       +EAT*(DR4*(FF(I6+1)-2.0*FF(I6)+FF(I6-1))
     6       +DZ4*(FF(I6+N1)-2.0*FF(I6)+FF(I6-N1)))
       F(JJ1+M4)=F(JJ1+M4)
     1       -DZ2*(U(I3)*U(I5)-U(I2)*U(I6)
     2       -U(I3-N1)*U(I5-N1)+U(I2-N1)*U(I6-N1)
     3       +U(I3-1)*U(I5-1)-U(I2-1)*U(I6-1)
     4       -U(I3-N11)*U(I5-N11)+U(I2-N11)*U(I6-N11))
     5       +EAT*(DR4*(FF(I5+1)-2.0*FF(I5)+FF(I5-1))
     6       +DZ4*(FF(I5+N1)-2.0*FF(I5)+FF(I5-N1)))
       X3=DR2*(U(I2)+U(I2-N1)-U(I2-1)-U(I2-N11))
     2    +DZ2*(U(I3)-U(I3-N1)+U(I3-1)-U(I3-N11))
       F(JJ1+M3)=F(JJ1+M3)-P(J4)*GAM*X3
     1       -DR2*P(J2)*(U(I4)+U(I4-N1)-U(I4-1)-U(I4-N11))
     2       -DZ2*P(J3)*(U(I4)-U(I4-N1)+U(I4-1)-U(I4-N11))
     3       +PMU*(DR4*(FF(I4+1)-2.0*FF(I4)+FF(I4-1))
     4       +DZ4*(FF(I4+N1)-2.0*FF(I4)+FF(I4-N1)))
       F(JJ1+M2)=F(JJ1+M2)+T*((-P(J7)*P(J5))/P(J1)
     1       -EUD*P(J3)-GRA*Z/RO2)
     2       -DR2*P(J2)*(U(I3)+U(I3-N1)-U(I3-1)-U(I3-N11))
     3       -DZ2*P(J3)*(U(I3)-U(I3-N1)+U(I3-1)-U(I3-N11))
     4       -DZ2*(U(I4)-U(I4-N1)+U(I4-1)-U(I4-N11))/P(J1)
     5       +RMU*(DR4*(FF(I3+1)-2.0*FF(I3)+FF(I3-1))
     6       +DZ4*(FF(I3+N1)-2.0*FF(I3)+FF(I3-N1)))/FF(I1)
       F(JJ1+N2)=F(JJ1+N2)+T*((P(J7)*P(J6))/P(J1)
     1       -EUD*P(J2)-GRA*R/RO2)
     2       -DR2*P(J2)*(U(I2)+U(I2-N1)-U(I2-1)-U(I2-N11))
     3       -DZ2*P(J3)*(U(I2)-U(I2-N1)+U(I2-1)-U(I2-N11))
     4       -DR2*(U(I4)+U(I4-N1)-U(I4-1)-U(I4-N11))/P(J1)
     5       +RMU*(DR4*(FF(I2+1)-2.0*FF(I2)+FF(I2-1))
     6       +DZ4*(FF(I2+N1)-2.0*FF(I2)+FF(I2-N1)))/FF(I1)
       F(JJ1)=F(JJ1)
     1       -DR2*(U(I2)*U(I1)+U(I2-N1)*U(I1-N1)
     2       -U(I2-1)*U(I1-1)-U(I2-N11)*U(I1-N11))
     3       -DZ2*(U(I3)*U(I1)-U(I3-N1)*U(I1-N1)
     4       +U(I3-1)*U(I1-1)-U(I3-N11)*U(I1-N11))
     5       +DFU*(DR4*(FF(I1+1)-2.0*FF(I1)+FF(I1-1))
     6       +DZ4*(FF(I1+N1)-2.0*FF(I1)+FF(I1-N1)))
   80  CONTINUE
C
C      SMOOTH CONNECTION WITH INNER IONOSPHERIC BOUNDARY
       DO 81 I=2,NR1
       R=HHX(I)+XX4
       RO2=R*R+Z*Z
       I1=I+N1
       I2=I1+K3
       I3=I2+K3
       I4=I3+K3
       I5=I4+K3
       I6=I5+K3
       JJ1=I+L1
       X2=RO2/AR2-1.0
       X2=AMAX1(X2,0.0)
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
C
       F(JJ1+M5)=F(JJ1+M5)*X2+FF(I6)*(1.0-X2)
       F(JJ1+M4)=F(JJ1+M4)*X2+FF(I5)*(1.0-X2)
       F(JJ1+M3)=F(JJ1+M3)*X2+FF(I4)*(1.0-X2)
       F(JJ1+M2)=F(JJ1+M2)*X2+FF(I3)*(1.0-X2)
       F(JJ1+N2)=F(JJ1+N2)*X2+FF(I2)*(1.0-X2)
       F(JJ1)=F(JJ1)*X2+FF(I1)*(1.0-X2)
   81  CONTINUE
C
C      PRESERVE OF F(JJ1) WITH P(I)
       DO 82 I=2,NR1
       JJ1=I+L1
       JJ2=JJ1+N2
       JJ3=JJ2+N2
       JJ4=JJ3+N2
       P(I+K3)=F(JJ4)
       P(I+K2)=F(JJ3)
       P(I+N1)=F(JJ2)
       P(I)=F(JJ1)
   82  CONTINUE
C
C      DIAGNOSTICS OF MAXIMUM VELOCITY
       X3=0.0
       DO 83 I=2,NR1
       J1=I
       J2=J1+N1
       J3=J2+N1
       J4=J3+N1
       JJ1=I+L1
       X3=AMAX1(X3,P(J2)*P(J2)+P(J3)*P(J3))
       F(JJ1+M3)=AMAX1(P(J4),PR02)
       F(JJ1)=AMAX1(P(J1),RO02)
   83  CONTINUE
       IF(X3.LT.1.0) GO TO 84
       WRITE(12) F
       WRITE(6,663) II,LAST,NR,NZ,X1,X2,X3
  663  FORMAT(1H ,5X,4I10,1P3E15.5)
       GO TO 300
   84  CONTINUE
C
C      MOVEMENT OF GRID POINTS IN U(I1) AND FF(I1)
       DO 90 M=1,NB
       L2=K3*(M-1)
       DO 90 I=1,N1
       I1=I+L2
       U(I1)=U(I1+N1)
       FF(I1)=FF(I1+N1)
   90  CONTINUE
C
C
       TIME=TIME+T
C
C      FINISH OF A TIME STEP INTEGRATION IN TWO STEP
C      LAX-WENDROFF METHOD
C
       IF(IIQ.NE.IIQ0) GO TO 401
       IIQ=0
C
C      UPSTREAM BOUNDARY CONDITION WITH SMOOTH TRANSITION
       X1=2.0
       BIS1=0.0
       IF(II.GE.IICH) BIS1=BIS
       NR11=NR1-2*NRP
       DO 102 J=1,NZ2
       L1=N1*(J-1)
       DO 102 I=1,NR2
       X=(2.0*HHX(I)/HR-FLOAT(NR11+2))/FLOAT(NR11)
       X=(X1+1.0)*X+X1
       X=AMAX1(-X,0.0)
       X=X*X
       I1=I+L1
       F(I1+M5)=(1.0-X)*F(I1+M5)+BIS1*X
       F(I1+M4)=(1.0-X)*F(I1+M4)
       F(I1+M3)=(1.0-X)*F(I1+M3)+PR01*X
       F(I1+M2)=(1.0-X)*F(I1+M2)
       F(I1+N2)=(1.0-X)*F(I1+N2)+VSW*X
       F(I1)=(1.0-X)*F(I1)+RO01*X
  102  CONTINUE
C
  401  IF(IIP.NE.IIP0) GO TO 100
       IIP=0
C
C      DIAGNOSTICS PARTS
C
       JXG=N2
       JYG=NZ2
       JXN=NR2
       LANK=20
       IP=0
       LB=1
       LL=6
       I1=II
C      CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F)
       DO 120 J=1,2
       LB=1+3*(J-1)
       LL=3
C      CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F)
  120  CONTINUE
C
C 501  IF(IIS.NE.IIS0) GO TO 100
       IF(IIS.NE.IIS0) GO TO 100
       IIS=0
C      SIMULATION ARGUMENTS ARE STORED AS SAMPLING
       WRITE(12) F
  100  CONTINUE
  300  CONTINUE
    9  CONTINUE
       REWIND 12
       STOP
       END
       SUBROUTINE EQUIP4(NR,NZ,NRP,RO01,PR01,CJ,CP,F,P)
C      INITIAL CONDITION OF 2-DIMENSIONAL DIPOLE MAGNETIC FIELD
       DIMENSION P(1),F(1),CJ(1),CP(1)
C      CJ 1-R1 2-DR 3-DZ 4-B0 5-RM 6-A2 7-AID 8-AR1 9-GAM 10-GRA
C      CP 1-R0 2-Z0 3-RA 4-VO0 5-P0 6-GRA 7-07 8-VSW 9-09 10-10
C
       AID=CJ(7)
       AR1=CJ(8)
       GAM=CJ(9)
       GRA=CJ(10)
       RATT=CJ(2)
       P0=(GAM-1.0)*GRA/GAM
       BIS=CP(10)*1.0E-4
       X1=1.0/(GAM-1.0)
       X2=GAM*X1
       B0=CJ(4)
       RO02=RATT*RO01
       PR02=PR01
C
       IC0=4
       ICD=6
       PI=3.1415926
       NR1=NR+1
       NR2=NR+2
       N1=NR2
       N11=N1+1
       NB=8
       NZ1=NZ+1
       NZ2=NZ+2
       N2=N1*NZ2
       HR=CP(1)/FLOAT(NR1)
       HZ=CP(2)/FLOAT(NZ1)
       VSW=CP(8)
C
       DO 10 J=1,NZ2
       VR1=VSW
       DO 10 I=1,NR2
       R=0.5*HR*FLOAT(2*I-NR2-1+2*NRP)
       Z=0.5*HZ*FLOAT(2*J-3)
       RO2=R*R+Z*Z
       RO1=SQRT(RO2)
       RO3=RO1*RO2
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       F(I1)=1.0/RO2
       IF(F(I1).LT.RO02) F(I1)=RO02
       X=FLOAT(IC0+ICD-I)/FLOAT(ICD)
       IF(X.GE.1.0) X=1.0
       IF(X.LE.0.0) X=0.0
       F(I2)=X*CP(8)
       F(I4)=1.0E-7*CP(5)/RO1
       IF(F(I4).LT.PR02) F(I4)=PR02
       F(I5)=-2.0*B0*R*Z/(RO2*RO2)
       F(I6)=B0*(R*R-Z*Z)/(RO2*RO2)
       X3=SQRT(F(I5)*F(I5)+F(I6)*F(I6))
       IF(VSW.LT.1.0E-6) GO TO 10
       X=(RO01-CJ(1)*(X3/VSW)**2)/F(I1)
       IF(X.LE.0.0) X=0.0
       X=VSW*SQRT(X)
       VR1=AMIN1(VR1,X)
       F(I2)=VR1
       XX=(1.0-RATT)*VR1/VSW
       F(I1)=F(I1)+RO01*XX
       F(I4)=F(I4)+PR01*XX
       F(I6)=F(I6)+BIS
   10  CONTINUE
C
C      CURRENT
       DO 28 J=1,NZ1
       DO 28 I=1,NR1
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       P(I1)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ
     1        -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR)
   28  CONTINUE
       RETURN
       END
       SUBROUTINE EQUIP5(NR,NZ,NRP,RO01,PR01,CJ,CP,F,P)
C      INITIAL CONDITION OF SUPERIMPOSED DIPOLE MAGNETIC FIELD
C      AS THE CHAPMAN-FERRARO MODEL
       DIMENSION P(1),F(1),CJ(1),CP(1)
C      CJ 1-R1 2-DR 3-DZ 4-B0 5-RM 6-A2 7-AID 8-AR1 9-GAM 10-GRA
C      CP 1-R0 2-Z0 3-RA 4-VO0 5-P0 6-GRA 7-07 8-VSW 9-09 10-10
C
       AID=CJ(7)
       AR1=CJ(8)
       GAM=CJ(9)
       GRA=CJ(10)
       RATT=CJ(2)
       P0=(GAM-1.0)*GRA/GAM
       BIS=CP(10)*1.0E-4
       X1=1.0/(GAM-1.0)
       X2=GAM*X1
       B0=CJ(4)
       RO02=RATT*RO01
       PR02=PR01
C
       IC0=4
       ICD=6
       PI=3.1415926
       NR1=NR+1
       NR2=NR+2
       N1=NR2
       N11=N1+1
       NB=8
       NZ1=NZ+1
       NZ2=NZ+2
       N2=N1*NZ2
       HR=CP(1)/FLOAT(NR1)
       HZ=CP(2)/FLOAT(NZ1)
       VSW=CP(8)
       ALP=CJ(5)
       BET=CJ(6)
       ALP1=1.0/SQRT(ALP)
       XM=(2.0*ALP*B0*B0/(RO01*VSW*VSW))**(1.0/4.0)
       XM=-XM
C
       DO 10 J=1,NZ2
       VR1=VSW
       DO 10 I=1,NR2
       R=0.5*HR*FLOAT(2*I-NR2-1+2*NRP)
       Z=0.5*HZ*FLOAT(2*J-3)
       RO2=R*R+Z*Z
       RO1=SQRT(RO2)
       RO3=RO1*RO2
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       F(I1)=1.0/RO2
       IF(F(I1).LT.RO02) F(I1)=RO02
       X=FLOAT(IC0+ICD-I)/FLOAT(ICD)
       IF(X.GE.1.0) X=1.0
       IF(X.LE.0.0) X=0.0
       F(I2)=X*CP(8)
       F(I4)=1.0E-7*CP(5)/RO1
       IF(F(I4).LT.PR02) F(I4)=PR02
       F(I5)=-2.0*B0*R*Z/(RO2*RO2)
       F(I6)=B0*(R*R-Z*Z)/(RO2*RO2)
   10  CONTINUE
C
C      CURRENT
       DO 28 J=1,NZ1
       DO 28 I=1,NR1
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       P(I1)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ
     1        -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR)
   28  CONTINUE
C
       DO 40 J=1,NZ2
       VR1=VSW
       DO 40 I=1,NR2
       R=0.5*HR*FLOAT(2*I-NR2-1+2*NRP)
       Z=0.5*HZ*FLOAT(2*J-3)
       R1M=R-2.0*XM
       RO2=R*R+Z*Z
       RO1=SQRT(RO2)
       RO3=RO1*RO2
       RO2M=R1M*R1M+Z*Z
       RO1M=SQRT(RO2M)
       RO3M=RO1M*RO2M
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       F(I1)=1.0/RO2
       IF(F(I1).LT.RO02) F(I1)=RO02
       X=FLOAT(IC0+ICD-I)/FLOAT(ICD)
       IF(X.GE.1.0) X=1.0
       IF(X.LE.0.0) X=0.0
       F(I2)=X*CP(8)
       F(I4)=1.0E-7*CP(5)/RO1
       IF(F(I4).LT.PR02) F(I4)=PR02
       BX1=-2.0*B0*R*Z/(RO2*RO2)
       BZ1=B0*(R*R-Z*Z)/(RO2*RO2)
       BX2=-2.0*B0*R1M*Z/(RO2M*RO2M)
       BZ2=B0*(R1M*R1M-Z*Z)/(RO2M*RO2M)
       F(I5)=BX1+BX2
       F(I6)=BZ1+BZ2
       IF(R.LT.XM) F(I5)=0.0
       IF(R.LT.XM) F(I6)=0.0
       Y1=ALP1+(1.0-ALP1)*(R-XM)/(BET-1.0)/XM
       IF(R.GT.XM) Y1=0.0
       Y1=AMAX1(Y1,0.0)
       Y1=AMIN1(Y1,1.0)
       VR1=VSW*Y1
C
       PX1=SQRT(F(I5)*F(I5)+F(I5)*F(I5))
       IF(ABS(VSW).LT.1.0E-8) GO TO 40
C      Y1=(RO01-CJ(1)*(PX1/VSW)**2)/RO01
C      IF(Y1.LE.0.0) Y1=0.0
C      Y1=VSW*SQRT(Y1)
C      VR1=AMIN1(VR1,Y1)
       F(I2)=VR1
       XX=VR1/VSW
       F(I4)=F(I4)+(PR01-PR02)*XX
       F(I6)=F(I6)+BIS*Y1
   40  CONTINUE
C
       RETURN
       END
=================== end of mhd2.f ================================

===================== beginning of mhd2graph.f =========================
C      FILE NAME  /JSSS/RDSPF801
       PARAMETER (X0=1.5,Y0=2.0,XL0=6.5,YL=6.5)
       PARAMETER (DXL=0.5,DYL=1.0)
       PARAMETER (NR= 40,NZ= 20)
       PARAMETER (NR1=NR+1,NR2=NR+2,NZ1=NZ+1,NZ2=NZ+2)
       PARAMETER (N1=NR2,N2=N1*NZ2,NB=6,N3=N2*NB)
       PARAMETER (NZZ=NZ*2,N11=N1+1)
       PARAMETER (NZZ1=NZZ+1,NZZ2=NZZ+2,NN2=N1*NZZ2,NN3=NN2*2)
       PARAMETER (NRG=20,NZG=20,ICU=1,KP=0,LANK=10)
       PARAMETER (LAST= 5,IIQ0=2,RRL=202.0,ZZL=102.0,MM=2,MOD=2)
       PARAMETER (NXP= 0,ARU=10.0,AR1=16.0,LAN1=90,LAN2=10)
       PARAMETER (IAR=1,URMIN=0.01,BIS=-0.0E-4,IXC=21,MDFU=1)
       DIMENSION U(NN3),F(N3),PP(N2),P(N2)
       DIMENSION IZM(9),ZMIN(9),ZMAX(9)
       DATA IZM/2,2,2,2,2,2,2,0,0/
       DATA ZMIN/-0.3,0.0E-8,0.E-4,  0.0,-0.044,-1.E-3,-1.E-4,8.,9./
       DATA ZMAX/ 0.3,1.0E-6,2.E-3,1.E-3, 0.044, 1.E-3, 1.E-4,8.,9./
       REWIND(12)
C
       CALL PLOTS
C
       DO 300 JJJ=1,3
C
       AR2=AR1*AR1
       PI=3.1415926
       HR=RRL/FLOAT(NR1)
       HZ=ZZL/FLOAT(NZ1)
       HR2=0.5*HR
       HZ2=0.5*HZ
       HRG=FLOAT(NR1)/FLOAT(NRG+1)
       HZG=FLOAT(2*NZ+1)/FLOAT(NZG+1)
       HRR=RRL/FLOAT(NRG+1)
       HZZ=ZZL/FLOAT(NZG+1)
       NRZG=NRG*NZG
       XEP=0.5*(FLOAT(NR1-2*NXP))/FLOAT(NR1)
       IIQ=IIQ0-1
       JJ=0
C      READ(12) F
C
       DO 100 II=1,LAST
       IIQ=IIQ+1
       READ(12,END=9) F
       IF(IIQ.NE.IIQ0) GO TO 100
       IIQ=0
       IF(II.LT.1) GO TO 100
C
       DO 10 II1=1,MOD
C      BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2
       DO 11 M=1,NB
       DO 11 I=2,NR1
       I1=I+N2*(M-1)
       I2=I1+N1*NZ1
       I3=N1*2
       F(I1)=F(I1+N1)
       IF(M.EQ.3) F(I1)=-F(I1)
       IF(M.EQ.5) F(I1)=-F(I1)
   11  F(I2)=(4.0*F(I2-N1)-F(I2-I3))/3.0
C      BOUNDARY CONDITION AT NR=1
       DO 13 J=1,NZ2
C      FREE BOUNDARY AT NR=NR2
       DO 13 M=1,NB
       I1=NR2+N1*(J-1)+N2*(M-1)
       F(I1)=(4.0*F(I1-1)-F(I1-2))/3.0
   13  CONTINUE
       IF(II.EQ.MOD) GO TO 10
       DO 15 M=1,NB
       DO 17 I=1,N2
       I1=I+N2*(M-1)
       U(I)=F(I1)
   17  CONTINUE
       DO 15 J=2,NZ1
       DO 15 I=2,NR1
       I1=I+N1*(J-1)
       I2=I1+N2*(M-1)
       F(I2)=0.125*(4.0*U(I1)+U(I1+1)+U(I1-1)+U(I1+N1)+U(I1-N1))
   15  CONTINUE
   10  CONTINUE
C
CC*
CC*
C
       DO 20 M=1,6
       DO 22 I=1,N1
       I1=I+N1+N2*(M-1)
       I2=I+N1*(M-1)
       U(I2)=F(I1)
       IF(M.EQ.2) U(I2)=-U(I2)
       IF(M.EQ.5) U(I2)=-U(I2)
       IF(M.EQ.3) U(I2)=0.0
       IF(M.EQ.5) U(I2)=0.0
   22  CONTINUE
       DO 24 K=2,NZ2
       I1=IXC+N1*(K-1)+N2*(M-1)
       I2=NZ2+K-2+N1*(M+5)
       I3=NZ2-K+1+N1*(M+5)
       U(I2)=0.5*(F(I1)+F(I1+1))
       U(I3)=0.5*(F(I1)+F(I1+1))
       IF(M.EQ.2) U(I2)=-U(I2)
       IF(M.EQ.2) U(I3)=-U(I3)
       IF(M.EQ.3) U(I3)=-U(I3)
       IF(M.EQ.5) U(I2)=-U(I2)
   24  CONTINUE
   20  CONTINUE
C
       DO 36 J=1,6
       DO 36 K=1,MDFU
       DO 37 I=1,N1
       I2=I+N1*(J-1)
       P(I)=U(I2)
   37  CONTINUE
       DO 36 I=2,NR1
       I2=I+N1*(J-1)
       U(I2)=0.25*(P(I-1)+2.0*P(I)+P(I+1))
   36  CONTINUE
C
       DO 46 J=1,6
       DO 46 K=1,MDFU
       DO 47 I=1,N1
       I2=I+N1*(J-1)+N1*6
       P(I)=U(I2)
   47  CONTINUE
       DO 46 I=2,NR1
       I2=I+N1*(J-1)+N1*6
       U(I2)=0.25*(P(I-1)+2.0*P(I)+P(I+1))
   46  CONTINUE
C
C       CALL XYOPEN(0,0.,0.,32.5,32.5,0)
C       CALL PLOTS('NAME$',32.5)
       CALL FACTOR(1.00)
       I1=II
       CALL DATA(1.5,0.5,LAST,I1,NXP,NR)
       IAA=1
       XL=1.5*XL0
C       1 DENSITY - RO
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0
       YB=Y0
       IB=1
       LL=1
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       2 PRESSURE-P
       IAA=2
       CALL NEWPEN(IPEN)
       IM=2
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0
       YB=Y0
       IB=4
       LL=1
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       3 VELOCITY - VX VY VZ
       IAA=1
       CALL NEWPEN(IPEN)
       IM=5
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0
       YB=Y0+YL+DYL
       IB=2
       LL=2
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       4 MAGNETIC FIELD - BX BY BZ
       CALL NEWPEN(IPEN)
       IM=6
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0
       YB=Y0+2.0*(YL+DYL)
       IB=5
       LL=2
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C
       IAA=1
       XL=1.5*XL0
C       1 DENSITY - RO
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0
       IB=1+6
       LL=1
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       2 PRESSURE-P
       IAA=2
       CALL NEWPEN(IPEN)
       IM=2
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0
       IB=4+6
       LL=1
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       3 VELOCITY - VX VY VZ
       IAA=1
       CALL NEWPEN(IPEN)
       IM=5
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0+YL+DYL
       IB=2+6
       LL=2
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       4 MAGNETIC FIELD - BX BY BZ
       CALL NEWPEN(IPEN)
       IM=6
       KP1=IZM(IM)
       VMIN=1.5*ZMIN(IM)
       VMAX=1.5*ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0+2.0*(YL+DYL)
       IB=5+6
       LL=2
       CALL GRAP7G(IAA,IB,LL,NR,NR,XB,YB,XL,YL,U,VMIN,VMAX,KP1)
C       CALL PLOTE
       CALL CHART
C
CC*
CC*
C
       X1=0.0
       DO 26 I=1,N1
       J1=I+N1*NZZ1
       I1=I+N1*NZ1+N2*5
       X2=0.5*HR*F(I1)
       IF(I.EQ.1) U(J1)=X1+X2
       IF(I.GE.2) U(J1)=U(J1-1)+X1+X2
       X1=X2
   26  CONTINUE
C
       DO 30 I=1,N1
       X2=0.0
       J1=I+N1*(NZZ2-1)
       J2=I
       I1=I+N1*(NZ2-1)+N2*4
       X1=-0.5*HZ*F(I1)
       U(J2)=U(J1)
       DO 30 J=2,NZ1
       J1=I+N1*(NZZ2-J)
       J2=I+N1*(J-1)
       I1=I+N1*(NZ2-J)+N2*4
       X2=-0.5*HZ*F(I1)
       U(J1)=U(J1+N1)-X1-X2
       U(J2)=U(J1)
       X1=X2
   30  CONTINUE
C
       I=N1/2-NXP
       J=NZZ2/2
       I1=I+N1*(J-1)
       X1=0.25*(U(I1)+U(I1+1)+U(I1+N1)+U(I1+N11))
       DO 32 I=1,NN2
       U(I)=U(I)-X1
   32  CONTINUE
C
C       CALL XYOPEN(0,0.,0.,32.5,32.5,0)
C       CALL PLOTS('NAME$',32.5)
       CALL FACTOR(1.00)
       I1=II
       CALL DATA(1.5,0.5,LAST,I1,NXP,NR)
       XL=XL0
C
C       1 FLUX
       IM=1
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0
       YB=Y0
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LAN1,VO)
C       2 PRESSURE-P
       K=4
       DO 121 J=2,NZ2
       DO 121 I=1,NR2
       I1=I+N1*(J-1)+N2*(K-1)
       J1=I+N1*(J+NZ-1)
       J2=I+N1*(NZ2-J)
       U(J1)=F(I1)
       U(J2)=U(J1)
  121  CONTINUE
       IM=2
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LAN2,VO)
C       3 DENSITY-RO
       K=1
       DO 131 J=2,NZ2
       DO 131 I=1,NR2
       I1=I+N1*(J-1)+N2*(K-1)
       J1=I+N1*(J+NZ-1)
       J2=I+N1*(NZ2-J)
       U(J1)=F(I1)
       U(J2)=U(J1)
  131  CONTINUE
       IM=3
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+2.0*(XL+DXL)
       YB=Y0
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LANK,VO)
C       4 MOD-B
       K=5
       DO 141 J=2,NZ2
       DO 141 I=1,NR2
       I1=I+N1*(J-1)+N2*(K-1)
       J1=I+N1*(J+NZ-1)
       J2=I+N1*(NZ2-J)
       U(J1)=SQRT(F(I1)*F(I1)+F(I1+N2)*F(I1+N2))
       U(J2)=U(J1)
  141  CONTINUE
       IM=4
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0
       YB=Y0+YL+DYL
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LANK,VO)
C       5 VR,VZ
       IM=5
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0+1.0*(YL+DYL)
       NZG1=NZG/2+1
       DO 40 J=NZG1,NZG
       DO 40 I=1,NRG
       R1=FLOAT(I)*HRG
       Z1=0.5*FLOAT(2*J-NZG-1)*HZG
       IR1=IFIX(R1)+1
       IZ1=IFIX(Z1)+1
       R1=R1-FLOAT(IR1-1)
       Z1=Z1-FLOAT(IZ1-1)
       I1=IR1+NR2*(IZ1-1)+N2
       I2=I1+N2
       J1=I+NRG*(J-1)
       J2=J1+NRZG
       J3=I+NRG*(NZG-J)
       J4=J3+NRZG
       U(J1)=(1.0-R1)*(1.0-Z1)*F(I1)+R1*Z1*F(I1+N1+1)
     1       +(1.0-R1)*Z1*F(I1+N1)+R1*(1.0-Z1)*F(I1+1)
       U(J2)=(1.0-R1)*(1.0-Z1)*F(I2)+R1*Z1*F(I2+N1+1)
     1       +(1.0-R1)*Z1*F(I2+N1)+R1*(1.0-Z1)*F(I2+1)
       U(J3)=U(J1)
       U(J4)=-U(J2)
   40  CONTINUE
       CALL GRPH1M(NRG,NZG,NXP,XB,YB,XL,YL,HRR,HZZ,AR2,U,VMIN,VMAX,
     1             URMIN,KP1,IAR,XEP)
C       6 BR,BZ
       IM=6
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+2.0*(XL+DXL)
       YB=Y0+YL+DYL
       NZG1=NZG/2+1
       DO 50 J=NZG1,NZG
       DO 50 I=1,NRG
       R1=FLOAT(I)*HRG
       Z1=0.5*FLOAT(2*J-NZG-1)*HZG
       IR1=IFIX(R1)+1
       IZ1=IFIX(Z1)+1
       R1=R1-FLOAT(IR1-1)
       Z1=Z1-FLOAT(IZ1-1)
       I1=IR1+NR2*(IZ1-1)+N2*4
       I2=I1+N2
       J1=I+NRG*(J-1)
       J2=J1+NRZG
       J3=I+NRG*(NZG-J)
       J4=J3+NRZG
       U(J1)=(1.0-R1)*(1.0-Z1)*F(I1)+R1*Z1*F(I1+N1+1)
     1       +(1.0-R1)*Z1*F(I1+N1)+R1*(1.0-Z1)*F(I1+1)
       U(J2)=(1.0-R1)*(1.0-Z1)*F(I2)+R1*Z1*F(I2+N1+1)
     1       +(1.0-R1)*Z1*F(I2+N1)+R1*(1.0-Z1)*F(I2+1)
       U(J3)=-U(J1)
       U(J4)=U(J2)
   50  CONTINUE
       CALL GRPH1M(NRG,NZG,NXP,XB,YB,XL,YL,HRR,HZZ,AR2,U,VMIN,VMAX,
     1             URMIN,KP1,IAR,XEP)
C
C      CURRENT
       DO 38 J=2,NZ1
       DO 38 I=2,NR1
       I1=I+N1*(J-1)
       I2=I1+N2
       I3=I2+N2
       I4=I3+N2
       I5=I4+N2
       I6=I5+N2
       J1=I+N1*(J+NZ-1)
       J2=I+N1*(NZ2-J)
       J3=J1+NN2
       J4=J2+NN2
       R=0.5*HR*FLOAT(2*I-NR1-2+2*NXP)
       Z=0.5*HZ*FLOAT(2*J-3)
       RO2=R*R+Z*Z
       X2=RO2/AR2-1.0
       IF(X2.LT.0.0) X2=0.0
       X2=ARU*X2*X2
       X2=X2/(X2+1.0)
       U(J1)=0.125*((F(I5+N11)+2.0*F(I5+N1)+F(I5+NR1)
     1        -F(I5-NR1)-2.0*F(I5-N1)-F(I5-N11))/HZ
     2        -(F(I6+N11)+2.0*F(I6+1)+F(I6-NR1)
     3        -F(I6+NR1)-2.0*F(I6-1)-F(I6-N11))/HR)
       U(J3)=0.125*((F(I2+N11)+2.0*F(I2+N1)+F(I2+NR1)
     1        -F(I2-NR1)-2.0*F(I2-N1)-F(I2-N11))/HZ
     2        -(F(I3+N11)+2.0*F(I3+1)+F(I3-NR1)
     3        -F(I3+NR1)-2.0*F(I3-1)-F(I3-N11))/HR)
       U(J1)=U(J1)*X2
       U(J3)=U(J3)*X2
       IF(II.EQ.1) PP(I1)=U(J1)
       IF(II.GT.1) U(J1)=U(J1)-PP(I1)
       U(J2)=U(J1)
       U(J4)=-U(J3)
   38  CONTINUE
C
C      BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2
       DO 33 K=1,2
       DO 31 I=2,NR1
       I1=I+NN2*(K-1)
       I2=I1+N1*NZZ1
       I3=N1*2
       U(I1)=(4.0*U(I1+N1)-U(I1+I3))/3.0
       U(I2)=(4.0*U(I2-N1)-U(I2-I3))/3.0
   31  CONTINUE
C      BOUNDARY CONDITION AT NR=1
C      FREE BOUNDARY AT NR=NR2
       DO 33 J=1,NZZ2
       I1=1+N1*(J-1)+NN2*(K-1)
       I2=NR1+I1
       U(I1)=(4.0*U(I1+1)-U(I1+2))/3.0
       U(I2)=(4.0*U(I2-1)-U(I2-2))/3.0
   33  CONTINUE
C
C       7 CURRENT-JY
       IM=7
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0
       YB=Y0+2.0*(YL+DYL)
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LANK,VO)
C       8 VORTICITY ROT(V)-Y
       DO 181 I=1,NN2
       U(I)=U(I+NN2)
  181  CONTINUE
       IM=8
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+XL+DXL
       YB=Y0+2.0*(YL+DYL)
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LANK,VO)
C       9 ELECTRIC FIELD  EY
       K=2
       DO 191 J=2,NZ2
       DO 191 I=1,NR2
       I1=I+N1*(J-1)+N2*(K-1)
       J1=I+N1*(J+NZ-1)
       J2=I+N1*(NZ2-J)
       U(J1)=F(I1+N2)*F(I1+3*N2)-F(I1)*F(I1+4*N2)
       U(J2)=U(J1)
  191  CONTINUE
       IM=9
       KP1=IZM(IM)
       VMIN=ZMIN(IM)
       VMAX=ZMAX(IM)
       XB=X0+2.0*(XL+DXL)
       YB=Y0+2.0*(YL+DYL)
       CALL GRAP4M(NR,NZZ,NXP,XB,YB,XL,YL,HR2,HZ2,AR2,U,VMIN,VMAX,
     1             KP1,LANK,VO)
C       CALL PLOTE
       CALL CHART
  100  CONTINUE
  300  CONTINUE
    9  CONTINUE
C       CALL XYCLOS
       CALL PLOT(0.,0.,999) 
       STOP
       END

c-------------------- suba.f --------------------------------------
C      /JSSS/SUBM2A
       SUBROUTINE PHASS(III,NR,NZ,F,U)
       DIMENSION F(1),U(1)
       NX=2*NR+3
       NY=2*NZ+3
       NR1=NR+1
       NR2=NR+2
       NZ1=NZ+1
       NZ2=NZ+2
       N1=NR2*NZ2
       N2=NX*NY
       NX3=NX*3
       DO 10 I=1,N2
       U(I)=0.0
   10  CONTINUE
       DO 20 J=1,NZ2
       DO 20 I=1,NR2
       I1=2*I-1+NX*(2*J-2)
       I2=I+NR2*(J-1)+N1*(III-1)
       U(I1)=F(I2)
   20  CONTINUE
       DO 30 J=1,NZ2
       J1=2*J-1
       I1=2
       K=I1+NX*(J1-1)
       U(K)=(3.0*U(K-1)+6.0*U(K+1)-U(K+3))/8.0
       I1=2*NR1
       K=I1+NX*(J1-1)
       U(K)=(-U(K-3)+6.0*U(K-1)+3.0*U(K+1))/8.0
       DO 30 I=2,NR
       I1=2*I
       K=I1+NX*(J1-1)
       U(K)=(-U(K-3)+9.0*U(K-1)+9.0*U(K+1)-U(K+3))/16.0
   30  CONTINUE
       DO 40 I=1,NX
       J1=2
       K=I+NX*(J1-1)
       U(K)=(3.0*U(K-NX)+6.0*U(K+NX)-U(K+NX3))/8.0
       J1=2*NZ1
       K=I+NX*(J1-1)
       U(K)=(-U(K-NX3)+6.0*U(K-NX)+3.0*U(K+NX))/8.0
       DO 40 J=2,NZ
       J1=2*J
       K=I+NX*(J1-1)
       U(K)=(-U(K-NX3)+9.0*U(K-NX)+9.0*U(K+NX)-U(K+NX3))/16.0
   40  CONTINUE
       RETURN
       END
       SUBROUTINE DIFFU(M,NX,NY,U,V)
       DIMENSION U(1),V(1)
C
       IF(M.EQ.99) GO TO 99
       ALP=1.0/8.0
       NX1=NX-1
       NY1=NY-1
       N1=NX*NY
       DO 100 II=1,M
       DO 12 I=1,N1
       V(I)=U(I)
   12  CONTINUE
       DO 14 J=2,NY1
       DO 14 I=2,NX1
       I1=I+NX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+V(I1-NX)+V(I1+NX))
   14  CONTINUE
       NY2=NX*(NY-2)
C
       DO 16 I=2,NX1
       I1=I
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1+NX))
       I1=I+NX*(NY-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(V(I1-1)+V(I1+1)+2.0*V(I1-NX))
   16  CONTINUE
       DO 18 J=2,NY1
       I1=1+NX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1+1)+V(I1-NX)+V(I1+NX))
       I1=NX+NX*(J-1)
       U(I1)=(1.0-4.0*ALP)*V(I1)+ALP*(2.0*V(I1-1)+V(I1-NX)+V(I1+NX))
   18  CONTINUE
C
  100  CONTINUE
   99  CONTINUE
       RETURN
       END
       SUBROUTINE DATA(X,Y,LY,LM,LD,LH)
C      X,Y ------- DATE NO HYOGI ITI
       CHARACTER DATE*16
       WRITE(DATE,100) LY,LM,LD,LH
       CALL SYMBOL(X,Y,0.35,DATE,0.0,16)
  100  FORMAT(I4,I4,I4,I4)
       RETURN
       END
       SUBROUTINE GAROHD(X0,Y0,X1,Y1,B,CR)
       C=CR*B
       EP0=1.0E-10
       EP1=1.0E-3
       EP2=1.0E3
       IF(ABS(X1-X0).LT.EP0) GO TO 2
       IF(ABS(Y1-Y0).LT.EP0) GO TO 3
       AA=(Y1-Y0)/(X1-X0)
       IF(ABS(AA).GT.EP2) GO TO 2
       IF(ABS(AA).LT.EP1) GO TO 3
       AA2=AA*AA
       A1=B/SQRT(1.0+AA2)
       A2=C/SQRT(1.0+1.0/AA2)
C
       XH=X1-A1
       IF(X1.LT.X0) XH=X1+A1
       X2=XH+A2
       X3=XH-A2
       YH=Y1+AA*(XH-X1)
       Y2=YH-(X2-XH)/AA
       Y3=YH-(X3-XH)/AA
       GO TO 1
    2  CONTINUE
       X2=X1-C
       X3=X1+C
       Y2=Y1-B
       IF(Y1.LT.Y0) Y2=Y1+B
       Y3=Y2
       GO TO 1
    3  CONTINUE
       Y2=Y1-C
       Y3=Y1+C
       X2=X1-B
       IF(X1.LT.X0) X2=X1+B
       X3=X2
    1  CONTINUE
C
       CALL PLOT(X0,Y0,3)
       CALL PLOT(X1,Y1,2)
       CALL PLOT(X2,Y2,3)
       CALL PLOT(X1,Y1,2)
       CALL PLOT(X3,Y3,2)
       RETURN
       END
       SUBROUTINE GCIRCL(X,Y,R)
       DIMENSION SC(10)
       DATA SC/0.0,0.174,0.342,0.5,0.643,0.766,0.866,0.94,0.985,1.0/
       CALL PLOT(X+R,Y,3)
       DO 10 I=1,9
       X1=X+SC(10-I)*R
       Y1=Y+SC(I+1)*R
       CALL PLOT(X1,Y1,2)
   10  CONTINUE
       DO 20 I=10,18
       X1=X-SC(I-8)*R
       Y1=Y+SC(19-I)*R
       CALL PLOT(X1,Y1,2)
   20  CONTINUE
       DO 30 I=19,27
       X1=X-SC(28-I)*R
       Y1=Y-SC(I-17)*R
       CALL PLOT(X1,Y1,2)
   30  CONTINUE
       DO 40 I=28,36
       X1=X+SC(I-26)*R
       Y1=Y-SC(37-I)*R
       CALL PLOT(X1,Y1,2)
   40  CONTINUE
       RETURN
       END
       SUBROUTINE GRECT(X,Y,A,CR)
       ACR=A*CR
       X1=X-A
       Y1=Y-ACR
       X2=X+A
       Y2=Y+ACR
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X1,Y2,2)
       CALL PLOT(X2,Y2,2)
       CALL PLOT(X2,Y1,2)
       CALL PLOT(X1,Y1,2)
       RETURN
       END
       SUBROUTINE GRPH1M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,
     &                   VMAX,URMI0,KP,IAR,XEP)
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
       NXY=NX*NY
       NL=NXY*2
       AD=1.0
       DX=XL/FLOAT(NX+1)
       DY=YL/FLOAT(NY+1)
       AR=DX*0.40
       CR=0.40
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 8 J=1,NY
       DO 8 I=1,NX
       Y=0.5*HY*FLOAT(2*J-NY-1)
       X=0.5*HX*FLOAT(2*I-NX-1+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 8
       I1=I+NX*(J-1)
       I2=I1+NXY
       X1=SQRT(U(I1)*U(I1)+U(I2)*U(I2))
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
    8  CONTINUE
       IF(KP.EQ.0) GO TO 126
       IF(KP.EQ.1) GO TO 122
       IF(KP.EQ.-1) GO TO 124
       UMIN=VMIN
       UMAX=VMAX
       GO TO 126
  122  UMAX=VMAX
       GO TO 126
  124  UMIN=VMIN
  126  CONTINUE
C
       UMIN=-UMIN
       UMAX=AMAX1(UMAX,UMIN)
       UMIN=-UMAX
       URMI=URMI0*UMAX
C
       DO 18 I=1,NXY
       I2=I+NXY
       R=SQRT(U(I)*U(I)+U(I2)*U(I2))
       IF(R.GT.UMAX) U(I)=U(I)*UMAX/R
       IF(R.GT.UMAX) U(I2)=U(I2)*UMAX/R
   18  CONTINUE
C      
       WRITE(M3,100) UMIN
       WRITE(M4,100) UMAX
  100  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
       X=X0+XEP*XL
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       Y=Y0+0.5*YL
       CALL PLOT(X,Y,3)
       X=X0+XL
       CALL PLOT(X,Y,2)
C
       NY1=NY/2
       DO 10 J=1,NY1
       Y1=Y0+DY*FLOAT(2*J-1)
       DO 12 I=1,NX
       X1=X0+DX*FLOAT(I)
       I1=I+NX*(2*J-2)
       I2=I1+NXY
       R=SQRT(U(I1)*U(I1)+U(I2)*U(I2))
       IF(R.LT.URMI) GO TO 12
       UX=AD*DX*U(I1)/UMAX
       UY=AD*DY*U(I2)/UMAX
       X2=X1+UX
       Y2=Y1+UY
       XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY))
       AR1=AR
       IF(IAR.EQ.1) AR1=AR*XX1*1.5
       CALL GAROHD(X1,Y1,X2,Y2,AR1,CR)
   12  CONTINUE
       Y1=Y1+DY
       DO 14 I=1,NX
       X1=X0+DX*FLOAT(NX+1-I)
       I1=NX+1-I+NX*(2*J-1)
       I2=I1+NXY
       R=SQRT(U(I1)*U(I1)+U(I2)*U(I2))
       IF(R.LT.URMI) GO TO 14
       UX=AD*DX*U(I1)/UMAX
       UY=AD*DY*U(I2)/UMAX
       X2=X1+UX
       Y2=Y1+UY
       XX1=SQRT((UX*UX+UY*UY)/(DX*DX+DY*DY))
       AR1=AR
       IF(IAR.EQ.1) AR1=AR*XX1*1.5
       CALL GAROHD(X1,Y1,X2,Y2,AR1,CR)
   14  CONTINUE
   10  CONTINUE
       RETURN
C
C
C
       ENTRY GRPH2M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,KP)
       NXY=NX*NY
       AD=0.5
       DX=XL/FLOAT(NX+1)
       DY=YL/FLOAT(NY+1)
       DR=AMIN1(DX,DY)
       DRR=DR*1.0E-2
       DRR=DR*2.0E-2
       DRR=AMAX1(DRR,5.0E-3)
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 20 J=1,NY
       DO 20 I=1,NX
       Y=0.5*HY*FLOAT(2*J-NY-1)
       X=0.5*HX*FLOAT(2*I-NX-1+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 20
       I1=I+NX*(J-1)
       X1=U(I1)
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
   20  CONTINUE
       IF(KP.EQ.0) GO TO 26
       IF(KP.EQ.1) GO TO 22
       IF(KP.EQ.-1) GO TO 24
       UMIN=VMIN
       UMAX=VMAX
       GO TO 26
   22  UMAX=VMAX
       GO TO 26
   24  UMIN=VMIN
   26  CONTINUE
C
       DO 28 I=1,NXY
       IF(U(I).GT.UMAX) U(I)=UMAX
       IF(U(I).LT.UMIN) U(I)=UMIN
   28  CONTINUE
C
       WRITE(M3,200) UMIN
       WRITE(M4,200) UMAX
  200  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,M2)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
C
       DI=0.0
       NY1=NY/2
       DO 30 J=1,NY1
       Y1=Y0+DY*FLOAT(2*J-1)
       DO 32 I=1,NX
       X1=X0+DX*FLOAT(I)
       I1=I+NX*(2*J-2)
       RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN)
       IF(RAD.LT.DRR) GO TO 32
       X2=X1
       Y2=Y1
       CALL GCIRCL(X2,Y2,RAD)
   32  CONTINUE
C
       Y1=Y1+DY
       DO 34 I=1,NX
       X1=X0+DX*FLOAT(NX+1-I)
       I1=NX+1-I+NX*(2*J-1)
       RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN)
       IF(RAD.LT.DRR) GO TO 34
       X2=X1
       Y2=Y1
       CALL GCIRCL(X2,Y2,RAD)
   34  CONTINUE
   30  CONTINUE
       RETURN
C
C
C
       ENTRY GRPH3M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,KP)
       NXY=NX*NY
       AD=0.5
       DX=XL/FLOAT(NX+1)
       DY=YL/FLOAT(NY+1)
       RR=DY/DX
       DR=AMIN1(DX,DY)
       DRR=DR*1.0E-2
       DRR=DR*2.0E-2
       DRR=AMAX1(DRR,5.0E-3)
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 40 J=1,NY
       DO 40 I=1,NX
       Y=0.5*HY*FLOAT(2*J-NY-1)
       X=0.5*HX*FLOAT(2*I-NX-1+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 40
       I1=I+NX*(J-1)
       X1=U(I1)
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
   40  CONTINUE
       IF(KP.EQ.0) GO TO 46
       IF(KP.EQ.1) GO TO 42
       IF(KP.EQ.-1) GO TO 44
       UMIN=VMIN
       UMAX=VMAX
       GO TO 46
   42  UMAX=VMAX
       GO TO 46
   44  UMIN=VMIN
   46  CONTINUE
C
       DO 48 I=1,NXY
       IF(U(I).GT.UMAX) U(I)=UMAX
       IF(U(I).LT.UMIN) U(I)=UMIN
   48  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
  300  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
C
       DI=0.0
       NY1=NY/2
       DO 50 J=1,NY1
       Y1=Y0+DY*FLOAT(2*J-1)
       DO 52 I=1,NX
       X1=X0+DX*FLOAT(I)
       I1=I+NX*(2*J-2)
       RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN)
       IF(RAD.LT.DRR) GO TO 52
       X2=X1
       Y2=Y1
       CALL GRECT(X2,Y2,RAD,RR)
   52  CONTINUE
C
       Y1=Y1+DY
       DO 54 I=1,NX
       X1=X0+DX*FLOAT(NX+1-I)
       I1=NX+1-I+NX*(2*J-1)
       RAD=AD*DR*(U(I1)-UMIN)/(UMAX-UMIN)
       IF(RAD.LT.DRR) GO TO 54
       X2=X1
       Y2=Y1
       CALL GRECT(X2,Y2,RAD,RR)
   54  CONTINUE
   50  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP3A(NR,NZ,MB,NX,NY,NXP,X0,Y0,XL,YL,
     &                  AR2,U,VMIN,VMAX,URMI0,KP,IAR,LAS1)
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
C
       NR1=NR+1
       NR2=NR+2
       NZ1=NZ+1
       NZ2=NZ+2
       N1=NR2
       N2=NR2*NZ2
       HX=FLOAT(NR+1)/FLOAT(NX+1)
       HY=FLOAT(NZ+1)/FLOAT(NY+1)
       RMIN=0.0
       RMAX=FLOAT(NR1)
       ZMIN=0.0
       ZMAX=FLOAT(NZ+1)
C
       NXY=NX*NY
       NL=NXY*2
       AD=3.0*HX/FLOAT(LAS1)
       DX=XL/FLOAT(NR+1)
       DY=YL/FLOAT(NZ+1)
       AR=DX*HX*0.60
       CR=0.40
C
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 8 J=1,NZ2
       DO 8 I=1,NR2
       Y=0.5*DY*FLOAT(2*J-NZ2-1)
       X=0.5*DX*FLOAT(2*I-NR2-1+2*NXP)
       R=X*X+Y*Y
C      IF(R.LT.AR2) GO TO 8
       I1=I+N1*(J-1)
       I2=I1+N2
       X1=SQRT(U(I1)*U(I1)+U(I2)*U(I2))
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
    8  CONTINUE
       IF(KP.EQ.0) GO TO 126
       IF(KP.EQ.1) GO TO 122
       IF(KP.EQ.-1) GO TO 124
       UMIN=VMIN
       UMAX=VMAX
       GO TO 126
  122  UMAX=VMAX
       GO TO 126
  124  UMIN=VMIN
  126  CONTINUE
C
       UMIN=-UMIN
       UMAX=AMAX1(UMAX,UMIN)
       UMIN=-UMAX
       URMI=URMI0*UMAX
C
       DO 18 I=1,N2
       I1=I+N2*(MB-1)
       I2=I1+N2
       R=SQRT(U(I1)*U(I1)+U(I2)*U(I2))
       IF(R.GT.UMAX) U(I1)=U(I1)*UMAX/R
       IF(R.GT.UMAX) U(I2)=U(I2)*UMAX/R
   18  CONTINUE
C      
       WRITE(M3,100) UMIN
       WRITE(M4,100) UMAX
  100  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
       X=X0+0.5*XL*FLOAT(NR1-2*NXP)/FLOAT(NR1)
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       Y=Y0+0.5*YL
       CALL PLOT(X,Y,3)
       X=X0+XL
       CALL PLOT(X,Y,2)
C
       NY1=NY/2+1
       DO 10 J=1,NY
       DO 12 I=1,NX
       R1=FLOAT(I)*HX
       Z1=FLOAT(J)*HY
       X1=X0+DX*R1
       Y1=Y0+DY*Z1
       X2=X1
       Y2=Y1
       RR=0
C
       DO 20 III=1,LAS1
       X1=X2
       Y1=Y2
       IR1=IFIX(R1)+1
       IZ1=IFIX(Z1)+1
       I1=IR1+NR2*(IZ1-1)
       I2=I1+N2
       R=R1-IFIX(R1)
       Z=Z1-IFIX(Z1)
       UX=R*Z*U(I1+N1+1)+(1.0-R)*Z*U(I1+N1)
     &   +R*(1.0-Z)*U(I1+1)+(1.0-R)*(1.0-Z)*U(I1)
       UY=R*Z*U(I2+N1+1)+(1.0-Z)*R*U(I2+N1)
     &   +R*(1.0-Z)*U(I2+1)+(1.0-R)*(1.0-Z)*U(I2)
       R=SQRT(UX*UX+UY*UY)
       IF(R.LT.0.01*URMI) GO TO 12
       IF(R.GT.UMAX) UX=UX*UMAX/R
       IF(R.GT.UMAX) UY=UY*UMAX/R
       UX=AD*UX/UMAX
       UY=AD*UY/UMAX
       IF(R.LT.UMAX) RR=RR+R/UMAX
       IF(R.GE.UMAX) RR=RR+1.0
       R2=R1+UX
       Z2=Z1+UY
       X2=X0+DX*R2
       Y2=Y0+DY*Z2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
       R1=R2
       Z1=Z2
       IF(R.LT.URMI) GO TO 22
       IF(R1.LE.RMIN.OR.R1.GE.RMAX) GO TO 22
       IF(Z1.LE.ZMIN.OR.Z1.GE.ZMAX) GO TO 22
   20  CONTINUE
   22  CONTINUE
       XX1=RR/FLOAT(LAS1)
       AR1=AR
       IF(IAR.EQ.1) AR1=AR*XX1*1.5
       CALL GAROHD(X1,Y1,X2,Y2,AR1,CR)
   12  CONTINUE
   10  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP4M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,
     &                   KP,LANK,V0)
C
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
C
       NXY=NX*NY
       NX1=NX+1
       NY1=NY+1
       NX2=NX+2
       NY2=NY+2
       N1=NX2
       N2=N1*NY2
C
       DX=XL/FLOAT(NX1)
       DY=YL/FLOAT(NY1)
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 10 J=1,NY2
       DO 10 I=1,NX2
       Y=0.5*HY*FLOAT(2*J-NY-3)
       X=0.5*HX*FLOAT(2*I-NX-3+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 10
       I1=I+NX2*(J-1)
       X1=U(I1)
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
   10  CONTINUE
       IF(KP.EQ.0) GO TO 16
       IF(KP.EQ.1) GO TO 12
       IF(KP.EQ.-1) GO TO 14
       IF(KP.EQ.5) GO TO 15
       UMIN=VMIN
       UMAX=VMAX
       GO TO 16
   12  UMAX=VMAX
       GO TO 16
   14  UMIN=VMIN
       GO TO 16
   15  CONTINUE
       X=ABS(UMAX-V0)
       Y=ABS(V0-UMIN)
       X=AMAX1(X,Y)
       UMAX=V0+X
       UMIN=V0-X
   16  CONTINUE
C
       DO 18 I=1,N2
       IF(U(I).GT.UMAX) U(I)=UMAX
       IF(U(I).LT.UMIN) U(I)=UMIN
   18  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
  300  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
       X=X0
       Y=Y0+0.5*YL
       CALL PLOT(X,Y,3)
       X=X0+XL
       CALL PLOT(X,Y,2)
       X=X0+0.5*XL*(FLOAT(NX1-2*NXP))/FLOAT(NX1)
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
C
       ULD=FLOAT(LANK+1)/(UMAX-UMIN)
       EP1=1.0E-8
C
       DO 100 J=1,NY1
       DO 100 I=1,NX1
C
       I1=I+N1*(J-1)
       I2=I1+1
       I3=I1+1+N1
       I4=I1+N1
C
       U1=(U(I1)-UMIN)*ULD
       U2=(U(I2)-UMIN)*ULD
       U3=(U(I3)-UMIN)*ULD
       U4=(U(I4)-UMIN)*ULD
C
       K1=IFIX(U1)
       K2=IFIX(U2)
       K3=IFIX(U3)
       K4=IFIX(U4)
C
       J1=IABS(K2-K1)
       J2=IABS(K3-K2)
       J3=IABS(K4-K3)
       J4=IABS(K1-K4)
C
       IF(J1.EQ.0) GO TO 28
       DO 20 I1=1,J1
       U0=FLOAT(I1)+FLOAT(MIN(K1,K2))
       IF(ABS(U2-U1).LT.EP1) GO TO 20
       X1=DX*((U0-U1)/(U2-U1)+FLOAT(I-1))
       Y1=DY*FLOAT(J-1)
       X1=X0+X1
       Y1=Y0+Y1
C
       IF((U3-U0)*(U2-U0).GT.0.0) GO TO 22
       IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 22
       IF(ABS(U3-U2).LT.EP1) GO TO 22
       X2=DX*FLOAT(I)
       Y2=DY*((U0-U2)/(U3-U2)+FLOAT(J-1))
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   22  CONTINUE
C
       IF((U4-U0)*(U3-U0).GT.0.0) GO TO 24
       IF((U1-U0)*(U3-U0).GT.0.0) GO TO 24
       IF((U2-U0)*(U4-U0).GT.0.0) GO TO 24
       IF(ABS(U4-U3).LT.EP1) GO TO 24
       X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1))
       Y2=DY*FLOAT(J)
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   24  CONTINUE
C
       IF((U1-U0)*(U4-U0).GT.0.0) GO TO 26
       IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 26
       IF(ABS(U1-U4).LT.EP1) GO TO 26
       X2=DX*FLOAT(I-1)
       Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1))
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   26  CONTINUE
   20  CONTINUE
   28  CONTINUE
C
       IF(J2.EQ.0) GO TO 38
       DO 30 I1=1,J2
       U0=FLOAT(I1)+FLOAT(MIN(K2,K3))
       IF(ABS(U3-U2).LT.EP1) GO TO 30
       X1=DX*FLOAT(I)
       Y1=DY*((U0-U2)/(U3-U2)+FLOAT(J-1))
       X1=X0+X1
       Y1=Y0+Y1
C
       IF((U4-U0)*(U3-U0).GT.0.0) GO TO 32
       IF(U0-U1.GT.0.0.AND.U0-U3.GT.0.0) GO TO 32
       IF(ABS(U4-U3).LT.EP1) GO TO 32
       X2=DX*((U0-U4)/(U3-U4)+FLOAT(I-1))
       Y2=DY*FLOAT(J)
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   32  CONTINUE
C
       IF((U1-U0)*(U4-U0).GT.0.0) GO TO 34
       IF((U1-U0)*(U3-U0).GT.0.0) GO TO 34
       IF((U2-U0)*(U4-U0).GT.0.0) GO TO 34
       IF(ABS(U1-U4).LT.EP1) GO TO 34
       X2=DX*FLOAT(I-1)
       Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1))
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   34  CONTINUE
   30  CONTINUE
   38  CONTINUE
C
       IF(J3.EQ.0) GO TO 48
       DO 40 I1=1,J3
       U0=FLOAT(I1)+FLOAT(MIN(K3,K4))
       IF(ABS(U4-U3).LT.EP1) GO TO 40
       X1=DX*((U0-U4)/(U3-U4)+FLOAT(I-1))
       Y1=DY*FLOAT(J)
       X1=X0+X1
       Y1=Y0+Y1
C
       IF((U1-U0)*(U4-U0).GT.0.0) GO TO 42
       IF(U0-U2.GT.0.0.AND.U0-U4.GT.0.0) GO TO 42
       IF(ABS(U1-U4).LT.EP1) GO TO 42
       X2=DX*FLOAT(I-1)
       Y2=DY*((U0-U1)/(U4-U1)+FLOAT(J-1))
       X2=X0+X2
       Y2=Y0+Y2
       CALL PLOT(X1,Y1,3)
       CALL PLOT(X2,Y2,2)
   42  CONTINUE
   40  CONTINUE
   48  CONTINUE
C
  100  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP5M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,
     &                   KP,NB1,NB2,V0)
C
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
C
       NXY=NX*NY
       NX1=NX+1
       NY1=NY+1
       NX2=NX+2
       NY2=NY+2
       N1=NX2
       N2=N1*NY2
C
       DX=XL/FLOAT(NX1)
       DY=YL/FLOAT(NY1)
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 10 J=1,NY2
       DO 10 I=1,NX2
       Y=0.5*HY*FLOAT(2*J-NY-3)
       X=0.5*HX*FLOAT(2*I-NX-3+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 10
       I1=I+NX2*(J-1)
       X1=U(I1)
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
   10  CONTINUE
       IF(KP.EQ.0) GO TO 16
       IF(KP.EQ.1) GO TO 12
       IF(KP.EQ.-1) GO TO 14
       IF(KP.EQ.5) GO TO 15
       UMIN=VMIN
       UMAX=VMAX
       GO TO 16
   12  UMAX=VMAX
       GO TO 16
   14  UMIN=VMIN
       GO TO 16
   15  CONTINUE
       X=ABS(UMAX-V0)
       Y=ABS(V0-UMIN)
       X=AMAX1(X,Y)
       UMAX=V0+X
       UMIN=V0-X
   16  CONTINUE
C
       DO 18 I=1,N2
       IF(U(I).GT.UMAX) U(I)=UMAX
       IF(U(I).LT.UMIN) U(I)=UMIN
   18  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
  300  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
       X=X0
       Y=Y0+0.5*YL
       CALL PLOT(X,Y,3)
       X=X0+XL
       CALL PLOT(X,Y,2)
       X=X0+0.5*XL*(FLOAT(NX1-2*NXP))/FLOAT(NX1)
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       ULD=YL/(UMAX-UMIN)
       DDXL=2.0*DX
C
       DO 100 M=NB1,NB2
       CALL NEWPEN(3*M-3*NB1+3)
       I1=1+N1*(M-1)
       U1=(U(I1)-UMIN)*ULD
       X1=X0
       Y1=Y0+U1
       CALL PLOT(X1,Y1,3)
       IF(M.GE.2) DDL=DDXL
C
       DO 100 I=2,N1
       I1=I+N1*(M-1)
       U1=(U(I1)-UMIN)*ULD
       X1=X0+DX*FLOAT(I-1)
       Y1=Y0+U1
       CALL PLOT(X1,Y1,2)
C      IF(M.EQ.2) CALL DASHP(X1,Y1,DDL)
C      IF(M.EQ.3) CALL PLOT(X1,Y1,2)
C      IF(M.GE.4) CALL DASHP(X1,Y1,DDL)
C
  100  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP6M(NX,NY,NXP,X0,Y0,XL,YL,HX,HY,AR2,U,VMIN,VMAX,
     &                   KP,NB1,NB2,V0)
C
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
C
       NXY=NX*NY
       NX1=NX+1
       NY1=NY+1
       NX2=NX+2
       NY2=NY+2
       N1=NX2
       N2=N1*NY2
C
       DX=XL/FLOAT(NX1)
       DY=YL/FLOAT(NY1)
       UMAX=U(1)+1.0E-10
       UMIN=U(1)-1.0E-10
       DO 10 J=1,NY2
       DO 10 I=1,NX2
       Y=0.5*HY*FLOAT(2*J-NY-3)
       X=0.5*HX*FLOAT(2*I-NX-3+2*NXP)
       R=X*X+Y*Y
       IF(R.LT.AR2) GO TO 10
       I1=I+NX2*(J-1)
       X1=U(I1)
       UMAX=AMAX1(UMAX,X1)
       UMIN=AMIN1(UMIN,X1)
   10  CONTINUE
       IF(KP.EQ.0) GO TO 16
       IF(KP.EQ.1) GO TO 12
       IF(KP.EQ.-1) GO TO 14
       IF(KP.EQ.5) GO TO 15
       UMIN=VMIN
       UMAX=VMAX
       GO TO 16
   12  UMAX=VMAX
       GO TO 16
   14  UMIN=VMIN
       GO TO 16
   15  CONTINUE
       X=ABS(UMAX-V0)
       Y=ABS(V0-UMIN)
       X=AMAX1(X,Y)
       UMAX=V0+X
       UMIN=V0-X
   16  CONTINUE
C
       DO 18 I=1,N2
       IF(U(I).GT.UMAX) U(I)=UMAX
       IF(U(I).LT.UMIN) U(I)=UMIN
   18  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
  300  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
       X=X0
       Y=Y0+0.5*YL
C      CALL PLOT(X,Y,3)
       X=X0+XL
C      CALL PLOT(X,Y,2)
       X=X0+0.5*XL*(FLOAT(NX1-2*NXP))/FLOAT(NX1)
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       ULD=YL/(UMAX-UMIN)
       DDXL=2.0*DX
C
       DO 100 M=NB1,NB2
       CALL NEWPEN(3*M-3*NB1+3)
       I1=1+N1*(M-1)
       U1=(U(I1)-UMIN)*ULD
       X1=X0
       Y1=Y0+U1
       CALL PLOT(X1,Y1,3)
       IF(M.GE.2) DDL=DDXL
C
       DO 100 I=2,N1
       I1=I+N1*(M-1)
       U1=(U(I1)-UMIN)*ULD
       X1=X0+DX*FLOAT(I-1)
       Y1=Y0+U1
       CALL PLOT(X1,Y1,2)
C      IF(M.EQ.2) CALL DASHP(X1,Y1,DDL)
C      IF(M.EQ.3) CALL PLOT(X1,Y1,2)
C      IF(M.GE.4) CALL DASHP(X1,Y1,DDL)
C
  100  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP6G(IA,LB,LL,NX,NXX,X0,Y0,XL,YL,U,VMIN,VMAX,KP)
C
       PARAMETER XM=5.0, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*4/'MIN='/, M2*4/'MAX='/, M3*8, M4*8
       CHARACTER M5*3/'II='/,M6*6
C
       NX1=NX+1
       NX2=NX+2
       N1=NXX+2
C
       INX=NXX/NX
       L1=LB+LL-1
       DO 8 J=LB,L1
       DO 8 I=1,NX
       I1=I+N1*(J-1)
       I2=I*INX+N1*(J-1)
       U(I1)=U(I2)
    8  CONTINUE
       DX=XL/FLOAT(NX1)
       I1=1+N1*(LB-1)
       UMAX=U(I1)+1.0E-6
       UMIN=U(I1)-1.0E-6
       DO 10 J=LB,L1
       DO 10 I=1,NX
       I1=I+N1*(J-1)
       UMAX=AMAX1(UMAX,U(I1))
       UMIN=AMIN1(UMIN,U(I1))
   10  CONTINUE
       IF(KP.EQ.0) GO TO 16
       IF(KP.EQ.1) GO TO 12
       IF(KP.EQ.-1) GO TO 14
       IF(KP.EQ.5) GO TO 15
       UMIN=VMIN
       UMAX=VMAX
       GO TO 16
   12  UMAX=VMAX
       GO TO 16
   14  UMIN=VMIN
       GO TO 16
   15  UMIN=-UMIN
       UMAX=AMAX1(UMAX,UMIN)
       UMIN=-UMAX
   16  CONTINUE
C
       DO 18 J=LB,L1
       DO 18 I=1,NX
       I1=I+N1*(J-1)
       IF(U(I1).GT.UMAX) U(I1)=UMAX
       IF(U(I1).LT.UMIN) U(I1)=UMIN
   18  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
       WRITE(M6,302) IA
  302  FORMAT(I6)
  300  FORMAT(1PE8.1)
       X1=X0+XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M5,0.0,3)
       CALL SYMBOL(999.,Y1,0.25,M6,0.0,6)
       CALL SYMBOL(X1,Y1,0.25,M1,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       X1=X1+XM
       CALL SYMBOL(X1,Y1,0.25,M2,0.0,4)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
C
       DO 20 JJ=LB,L1
       X=X0
       I1=NX+N1*(JJ-1)
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,3)
       DO 30 I=1,NX
       I1=I+N1*(JJ-1)
       X=X+DX
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,2)
   30  CONTINUE
       I1=1+N1*(JJ-1)
       X=X+DX
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,2)
   20  CONTINUE
       RETURN
       END
       SUBROUTINE GRAP7G(IA,LB,LL,NX,NXX,X0,Y0,XL,YL,U,VMIN,VMAX,KP)
C
       PARAMETER XM=3.4, YM=-0.6
       DIMENSION U(1)
       CHARACTER M1*6/'  MIN='/, M2*6/'  MAX='/, M3*8, M4*8
       CHARACTER M5*3/'II='/,M6*6
C
       NX1=NX+1
       NX2=NX+2
       N1=NXX+2
C
       INX=NXX/NX
       L1=LB+LL-1
       DO 8 J=LB,L1
       DO 8 I=1,NX
       I1=I+N1*(J-1)
       I2=I*INX+N1*(J-1)
       U(I1)=U(I2)
    8  CONTINUE
       DX=XL/FLOAT(NX1)
       I1=1+N1*(LB-1)
       UMAX=U(I1)+1.0E-6
       UMIN=U(I1)-1.0E-6
       DO 10 J=LB,L1
       DO 10 I=1,NX
       I1=I+N1*(J-1)
       UMAX=AMAX1(UMAX,U(I1))
       UMIN=AMIN1(UMIN,U(I1))
   10  CONTINUE
       IF(KP.EQ.0) GO TO 16
       IF(KP.EQ.1) GO TO 12
       IF(KP.EQ.-1) GO TO 14
       IF(KP.EQ.5) GO TO 15
       UMIN=VMIN
       UMAX=VMAX
       GO TO 16
   12  UMAX=VMAX
       GO TO 16
   14  UMIN=VMIN
       GO TO 16
   15  UMIN=-UMIN
       UMAX=AMAX1(UMAX,UMIN)
       UMIN=-UMAX
   16  CONTINUE
C
       DO 18 J=LB,L1
       DO 18 I=1,NX
       I1=I+N1*(J-1)
       IF(U(I1).GT.UMAX) U(I1)=UMAX
       IF(U(I1).LT.UMIN) U(I1)=UMIN
   18  CONTINUE
C
       WRITE(M3,300) UMIN
       WRITE(M4,300) UMAX
       WRITE(M6,302) IA
  302  FORMAT(I6)
  300  FORMAT(1PE8.1)
       X1=X0+0.8*XM
       Y1=Y0+YM
       CALL SYMBOL(X0,Y1,0.25,M5,0.0,3)
       CALL SYMBOL(999.,Y1,0.25,M6,0.0,6)
       X1=X1
       CALL SYMBOL(x1,Y1,0.25,M1,0.0,6)
       CALL SYMBOL(999.,Y1,0.25,M3,0.0,8)
       X1=X1+XM
       CALL SYMBOL(x1,Y1,0.25,M2,0.0,6)
       CALL SYMBOL(999.,Y1,0.25,M4,0.0,8)
C
       CALL PLOT(X0,Y0,3)
       X=X0+XL
       Y=Y0
       CALL PLOT(X,Y,2)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
       X=X0
       CALL PLOT(X,Y,2)
       Y=Y0
       CALL PLOT(X,Y,2)
C
       X=X0
       Y=Y0+0.5*YL
       CALL PLOT(X,Y,3)
       X=X0+XL
       CALL PLOT(X,Y,2)
       X=X0+0.5*XL
       Y=Y0
       CALL PLOT(X,Y,3)
       Y=Y0+YL
       CALL PLOT(X,Y,2)
C
       DO 20 JJ=LB,L1
       X=X0
       I1=1+N1*(JJ-1)
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,3)
       DO 30 I=2,NX1
       I1=I+N1*(JJ-1)
       X=X+DX
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,2)
   30  CONTINUE
       I1=N1+N1*(JJ-1)
       X=X+DX
       Y=(U(I1)-UMIN)/(UMAX-UMIN)
       Y=Y0+YL*Y
       CALL PLOT(X,Y,2)
   20  CONTINUE
       RETURN
       END
       SUBROUTINE PHASP(III,NX,NY,NXX,VE,F,U)
       DIMENSION F(1),U(1)
       PI=3.1415926
       PI1=1.0/SQRT(2.0*PI)
       N1=NXX+2
       NY2=NY/2
       INX=NXX/NX
       N2=NX*NY
       DO 10 I=1,N2
       U(I)=0.0
   10  CONTINUE
C
       DO 20 I=1,NX
       DO 20 J=1,NY
       V=FLOAT(J-NY2)*VE/FLOAT(NY2)
       DO 20 II=1,2
       I1=I*INX+3*N1*(II-1)+6*N1*(III-1)
       XN=F(I1)
       XV=F(I1+N1)
       XVT2=F(I1+2*N1)
       XVT=SQRT(XVT2)
       X=0.5*(V-XV)*(V-XV)/XVT2
       A=0.0
       IF(X.LT.5.0) A=PI1*XN*EXP(-X)/XVT
       I2=I+NX*(J-1)
       U(I2)=U(I2)+A
   20  CONTINUE
       RETURN
       END
       SUBROUTINE SCALL1(JXG,JYG,NX2,YMIN,YMAX,P)
       DIMENSION P(1)
       YMIN=P(1)
       YMAX=P(1)
       DO 532 J=1,JYG
       DO 532 I=1,JXG
       I1=I+NX2*(J-1)
       Y=P(I1)
       YMIN=AMIN1(YMIN,Y)
       YMAX=AMAX1(YMAX,Y)
  532  CONTINUE
       RETURN
       END
===================== end of mhd2graph.f =========================