C FILE:  THIS FILE CONTAINS THE PORTABLE SUPPORT ROUTINES FOR
C MICROSCOPE.  THEY HAVE PASSED THE PFORT VERIFIER WITHOUT ANY
C ERROR MESSAGES.
C
      SUBROUTINE BNDRY(NH,NW,SHIFT,IOPSS,IOPWN,WIDTH,LEFT,RIGHT)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  BNDRY
C  PURPOSE:  GIVEN THE CURRENT SAMPLING PARAMETERS AND THE USER
C             SPECIFIED MODIFICATIONS (I.E. STEP-SIZE, WINDOW WIDTH, AND
C             SHIFT) THIS ROUTINE COMPUTES THE RESULTING LEFT AND RIGHT
C             BOUNDARIES.
C
      INTEGER             IFN
      INTEGER             ICENTR,      IOPSS,       IOPWN,       IXSIZE
      INTEGER             LEFT,        NH,          NHP,         NW
      INTEGER             NWP,         RIGHT,       SHIFT,       W2
      INTEGER             WIDTH
      DATA ICENTR,IXSIZE / 2689, 5377 /
      NHP   = IFN(NH,IOPSS)
      NWP   = IFN(NW,IOPWN)
      W2    = WIDTH/2
      LEFT  = ICENTR+(SHIFT-W2)*NHP-NWP
      RIGHT = ICENTR+(SHIFT+W2)*NHP+NWP
      RETURN
      END
      SUBROUTINE CHKCMP(LCOMP,LMAG,LSHIFT,LSAMPL,LZERO,NUM,DEN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  CHKCMP
C  PURPOSE:  TO CHECK THE NEW LIST OF OPTIONS TO SEE IF THE
C             INTERPOLATION FUNCTION (F) MUST BE SAMPLED.  IF SO,
C             THE LOGICAL VARIABLE LCOMP IS SET TO .TRUE..
C             CHKCMP ALSO CHECKS TO SEE IF THE SAMPLING ARRAY MUST
C             BE ALTERED DUE TO OPTION CHANGES.
C
      INTEGER             I,           W2,          NUM,         DEN
      INTEGER             IXSIZE,      ICENTR
      LOGICAL             LMAG,        LSHIFT,      LSAMPL,      LCOMP
      LOGICAL             OLDNML,      LZERO
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             M,           V,           SHIFT,       NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
      LOGICAL             LFO
      COMMON / FOOWN  /   LFO
      EQUIVALENCE (SHIFT,IOPSH(2)),(M,IOPSS(2)),(V,IOPWN(2)),
     X    (NDIM,IOPDM(2))
      DATA ICENTR,IXSIZE / 2689, 5377 /
C
C  SET THE LOGICAL VARIABLES TO .FALSE.
C
      LCOMP  = .FALSE.
      LMAG   = .FALSE.
      LSHIFT = .FALSE.
      LSAMPL = .FALSE.
      LZERO  = .FALSE.
C
C  CHECK TO SEE IF ANY COMPUTATION HAS BEEN DONE SO FAR
C
      IF ( OPCMP(1) ) GO TO 100
      LCOMP = .TRUE.
      LSAMPL= .TRUE.
C
C  CHECK TO SEE IF ANY HIGHER DERIVATIVES NEED TO BE CALCULATED THUS
C   REQUIRING MORE FUNCTION EVALUATIONS
C
  100      DO 200 I = 1,7
           IF ( .NOT.((OPDF1(I,2).OR.OPDF2(I,2)) .AND. .NOT.LDF(I)) )
     X          GO TO 200
           LCOMP = .TRUE.
  200      CONTINUE
C
C  CHECK TO SEE IF SHIFTING OR SCALING REQUIRES RECOMPUTATION
C
      IF ( IOPSH(1).NE.IOPSH(2) ) LCOMP = .TRUE.
      CALL BNDRY(NH(1),NW(1),SHIFT,M,V,WIDTH,ILEFT(2),IRIGHT(2))
      IF ( .NOT.(ILEFT(2).LT.1 .OR. IRIGHT(2).GT.IXSIZE) ) GO TO 300
      IF ( SHIFT.NE.0 ) LSHIFT = .TRUE.
      IF ( M.NE.1 .OR. V.NE.1 ) LMAG = .TRUE.
  300 IF ( .NOT.(ILEFT(2).LT.ILEFT(1) .OR. IRIGHT(2).GT.IRIGHT(1)) ) GO
     X     TO 400
      LCOMP = .TRUE.
  400 IF ( M.NE.1 .AND. SHIFT.NE.0 ) LSHIFT = .TRUE.
C
C  CHECK TO SEE IF THE STEP-SIZE HAS CHANGED DUE TO EITHER NUM OR DEN
C   NOT BEING 1.  IF SO, SET LMAG = .TRUE. AND LCOMP = .TRUE.
C
      IF ( NUM.EQ.1 .AND. DEN.EQ.1 ) GO TO 500
      LMAG  = .TRUE.
      LCOMP = .TRUE.
C
C  CHECK TO SEE IF THE STENCIL WIDTH HAS BEEN CHANGED
C
  500 IF ( NW(1).EQ.NW(2) ) GO TO 600
      LCOMP = .TRUE.
C
C  CHECK TO SEE IF SUBZO HAS BEEN CALLED
C
  600 IF (IOPSS(2).EQ.IOPSS(1)) GO TO 700
      LCOMP = .TRUE.
C
C  CHECK TO SEE IF THE STEP-SIZE ROPSTS(2) HAS BEEN CHANGED
C
  700 IF ( ROPSTS(2).EQ.ROPSTS(1) ) GO TO 800
      LCOMP  = .TRUE.
      LSAMPL = .TRUE.
      IF ( M.NE.1 .OR. V.NE.1 ) GO TO 800
      M      = 1
      LZERO  = .TRUE.
C
C  CHECK TO SEE IF THE DIRECTION VECTOR HAS BEEN CHANGED
C
  800 IF (ROPUDI(1,1).EQ.(-ROPUDI(1,2)).AND.ROPUDI(2,1).EQ.(-ROPUDI(2,2)
     X     ).AND. ROPUDI(3,1).EQ.(-ROPUDI(3,2)) ) GO TO 900
      IF ( ROPUDI(1,1).EQ.ROPUDI(1,2) .AND. ROPUDI(2,1).EQ.ROPUDI(2,2)
     X     .AND. ROPUDI(3,1).EQ.ROPUDI(3,2) ) GO TO 900
      LCOMP = .TRUE.
      LMAG  = .FALSE.
      LSHIFT= .FALSE.
      LZERO = .TRUE.
      IOPSH(2)  = 0
C
C  CHECK TO SEE IF THE CENTER POINT HAS BEEN CHANGED
C
  900 IF ( ROPPNT(1,1).EQ.ROPPNT(1,2) .AND. ROPPNT(2,1).EQ.ROPPNT(2,2)
     X     .AND. ROPPNT(3,1).EQ.ROPPNT(3,2) ) GO TO 1000
      LCOMP = .TRUE.
      LMAG  = .FALSE.
      LSHIFT= .FALSE.
      LZERO = .TRUE.
      IOPSH(2)  = 0
 1000 CONTINUE
C
C CHECK IF THE CROSS DERIVATIVE HAS CHANGED
C
      IF (.NOT.LCBD) GO TO 1100
      LCBD = .FALSE.
      LCOMP = .TRUE.
      LMAG = .FALSE.
      LSHIFT = .FALSE.
      IOPSH(2) = 0
      LSAMPL = .TRUE.
      LZERO = .TRUE.
 1100 CONTINUE
C
C CHECK IF  FO  HAS BEEN CALLED
C
      IF (.NOT.LFO) GO TO 1200
      LFO = .FALSE.
      CALL MAKSHF
      LCOMP = .TRUE.
      LMAG = .FALSE.
      LSHIFT = .FALSE.
      IOPSH(2) = 0
      LSAMPL = .TRUE.
      LZERO = .TRUE.
 1200 CONTINUE
C
C CHECK IF NORMALIZATION HAS CHANGED
C
      IF ((OLDNML.AND.NORMAL).OR.(.NOT.OLDNML.AND..NOT.NORMAL))
     X  GO TO 1300
      LCOMP = .TRUE.
      LSAMPL = .TRUE.
      LMAG  = .FALSE.
      LSHIFT= .FALSE.
      LZERO = .TRUE.
      OLDNML = NORMAL
 1300 CONTINUE
      RETURN
      END
      INTEGER FUNCTION COPY(M,N)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  COPY
C  PURPOSE:  FOR A SAMPLING SCHEME WHICH IS A POWER OF 2 (I.E.POWER = M)
C             FINER THAN THE EXISTING SCHEME, COPY RETURNS THE INDEX
C             OF THE ELEMENT IN THE FS() ARRAY WHICH ELEMENT N IN THE
C             FS() ARRAY MAPS TO.
C
      INTEGER             MOD
      INTEGER             ICENTR,      IXSIZE,      K,           L(10)
      INTEGER             M,           N
      DATA L(1),L(2),L(3),L(4),L(5)  /  1,  1,  3,  5, 11 /
      DATA L(6),L(7),L(8),L(9),L(10) / 21, 43, 85,171,341 /
      DATA ICENTR,IXSIZE / 2689, 5377 /
      K     = N-ICENTR
      COPY  = ICENTR+2**M*K
      IF ( MOD(K,2).EQ.0 ) RETURN
      COPY  = COPY+L(M)*(-1)**(K/2)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DERIV(F,NDIM,P,IORDER,DIR,HH)
C EVALUATE THE IORDER-TH DERIVATIVE OF THE FUNCTION F OF NDIM
C VARIABLES IN THE DIRECTION DIR AT THE POINT P USING STEP-SIZE
C H
C
C THIS FUNCTION MAY BE EMPLOYED BY A USER OF MICROSCOPE
C
      DOUBLE PRECISION    F,           HH
      DOUBLE PRECISION    H,           P(1),        W1(3),       H2
      DOUBLE PRECISION    W2(3),       H3,          W3(3),       H4
      DOUBLE PRECISION    W4(3),       W5(3),       H6,          W6(3)
      DOUBLE PRECISION    W7(3),       W8(3),       W9(3),       DIR(1)
      DOUBLE PRECISION    ZERO
      INTEGER             I,           IORDER,      NDIM
      IF (IORDER.GE.0.AND.IORDER.LE.6.AND.NDIM.GE.1.AND.NDIM.LE.3) GO
     X     TO 100
      DERIV = 0.0D0
      GO TO 1200
  100 CONTINUE
      H = 2.0D0*HH
      IF (IORDER.GT.0) GO TO 200
      DERIV = F(P)
      GO TO 1200
  200 CONTINUE
      H2 = HH
           DO 300 I = 1,NDIM
           W1(I) = P(I)-H2*DIR(I)
           W9(I) = P(I)+H2*DIR(I)
  300      CONTINUE
      IF (IORDER.GT.1) GO TO 400
      DERIV = (F(W9)-F(W1))/H
      GO TO 1200
  400 CONTINUE
           DO 500 I = 1,NDIM
           W5(I) = P(I)
  500      CONTINUE
      IF (IORDER.GT.2) GO TO 600
      DERIV = (F(W9)-2.0D0*F(W5)+F(W1))/(H2**2)
      GO TO 1200
  600 CONTINUE
      H4 = H/4.0D0
           DO 700 I = 1,NDIM
           W3(I) = P(I)-H4*DIR(I)
           W7(I) = P(I)+H4*DIR(I)
  700      CONTINUE
      IF (IORDER.GT.3) GO TO 800
      DERIV = (-2.0D0*F(W1)+4.0D0*F(W3)-4.0D0*F(W7)+2.0D0*F
     X     (W9))/(4.0D0*H4**3)
      GO TO 1200
  800 CONTINUE
      IF (IORDER.GT.4) GO TO 900
      DERIV = (F(W1)-4.0D0*F(W3)+6.0D0*F(W5)-4.0D0*F(W7)
     X     +F(W9))/H4**4
      GO TO 1200
  900 CONTINUE
      H3 = H/3.0D0
      H6 = H/6.0D0
           DO 1000 I = 1,NDIM
           W2(I) = P(I)-H3*DIR(I)
           W8(I) = P(I)+H3*DIR(I)
           W4(I) = P(I)-H6*DIR(I)
           W6(I) = P(I)+H6*DIR(I)
 1000      CONTINUE
      IF (IORDER.GT.5) GO TO 1100
      DERIV = (-F(W1)+4.0D0*F(W2)-5.0D0*F(W4)+5.0D0*F(W6
     X     )-4.0D0*F(W8)+F(W9))/(2.0D0*H6**5)
      GO TO 1200
 1100 CONTINUE
      DERIV = (F(W1)-6.0D0*F(W2)+15.0D0*F(W4)-20.0D0*F(W5)
     X     +15.0D0*F(W6)-6.0D0*F(W8)+F(W9))/H6**6
 1200 CONTINUE
      RETURN
      END
      SUBROUTINE DFAULT
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  DFAULT
C  PURPOSE:  TO SET CERTAIN VARIABLES IN THE COMMON BLOCKS  OPTION,
C             USER, AND CB TO THEIR DEFAULT VALUES.
C
      LOGICAL             FIRST
      INTEGER             I,           ICENTR,      IXSIZE,      J
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
C     COMMON BLOCK / CB     /
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      DOUBLE PRECISION    ETA
      INTEGER             IROUND,      N
      LOGICAL             ADD
      COMMON / USER   /   ETA,         IROUND,      N,           ADD
C
      DATA ICENTR,IXSIZE / 2689, 5377 /
C
C  FILL IN THE VALUES OF THE VARIABLES IN THE ORDER IN WHICH THEY ARE
C   LISTED IN THE COMMON BLOCK.
C
           DO 200 I = 1,2
           ROPSTS(I)   = 0.0001D0
           NH(I)       = 8
           NW(I)       = 48
           ILEFT(I)    = ICENTR
           IRIGHT(I)   = ICENTR
           IOPSH(I)    = 0
           IOPSS(I)    = 1
           IOPWN(I)    = 1
           OPDC(I)     = .FALSE.
           OPDS(I)     = .FALSE.
           OPDX(I)     = .FALSE.
           OPCMP(I)    = .FALSE.
           OPSPL(I)    = .FALSE.
           OPDF1(1,I)  = .TRUE.
           OPDF2(1,I)  = .FALSE.
                DO 100 J = 2,7
                OPDF1(J,I) = .FALSE.
                OPDF2(J,I) = .FALSE.
  100           CONTINUE
  200      CONTINUE
      LCBD = .FALSE.
      ICBD = 0
      CBDSW = 0.0D0
           DO 300 I = 1,3
           CBDU(I) = 1.0D0
  300      CONTINUE
      IROUND = 10
      N = 1
      ADD = .FALSE.
      ETA = 1.0D0
      FIRST = .FALSE.
      RETURN
      END
      DOUBLE PRECISION FUNCTION DFN(N,M)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  DFN
C  PURPOSE:  TO CALCULATE THE DOUBLE PRECISION VALUE FOR THE NEW N GIVEN
C             THE OLD VALUE FOR N AND THE USER SPECIFIED MAGNIFICATION M
C
      DOUBLE PRECISION    DM,          DN
      INTEGER             M,           N
      DN  = N
      DM  = M
      DFN = DN*DM
      IF ( M.LT.0 ) DFN = -DN/DM
      RETURN
      END
      INTEGER FUNCTION DIGIT(N)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  DIGIT
C  PURPOSE:  GIVEN A SINGLE INTEGER ARGUMENT, REPRESENTING A HOLLERITH
C             CHARACTER (I.E. CHARACTERS READ IN BY USING AN A1 FORMAT),
C             THE VALUE OF DIGIT IS THE INTEGER VALUE CORRESPONDING TO
C             THE CHARACTER (I.E. 0 THRU 9), AND IF THE CHARACTER IS NOT
C             A NUMBER, DIGIT HAS THE VALUE -1.
C
      INTEGER             M(10)
      INTEGER             I,           N
      DATA  M(1),M(2),M(3),M(4),M(5)  / 1H0, 1H1, 1H2, 1H3, 1H4 /
      DATA  M(6),M(7),M(8),M(9),M(10) / 1H5, 1H6, 1H7, 1H8, 1H9 /
      DIGIT = -1
           DO 100 I = 1,10
           IF ( N.NE.M(I) ) GO TO 100
           DIGIT = I-1
           GO TO 200
  100      CONTINUE
  200 RETURN
      END
      SUBROUTINE DIREAD(DEVICE,IARG,IARGP,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  DIREAD
C  PURPOSE:  TO READ IN TWO INTEGER NUMBERS, FREE FORMATTED,
C             FROM THE SPECIFIED DEVICE.
C
      INTEGER             LENGTH
      INTEGER             BEGIN,       CHAR(72),    DEVICE,      ENDE
      INTEGER             I,           I72,         IARG,        IARGP
      LOGICAL             ERR
      DATA  I72 / 72 /
C
C  SET ERR = .TRUE. (I.E. NO ERRORS)
C
      ERR  = .TRUE.
C
C  READ CHARACTER STRING FROM INPUT DEVICE AND RECOGNIZE NUMBER
C
      IARG = 0
      READ (DEVICE,10000) (CHAR(I),I=1,72)
      ENDE   = LENGTH(I72,CHAR)
      IF ( ENDE.EQ.0 ) GO TO 200
      BEGIN= 1
      CALL DRINT(I72,CHAR,BEGIN,ENDE,IARG,IARGP,ERR)
  200 RETURN
C
C  FORMAT STATEMENT
C
10000 FORMAT(72A1)
      END
      SUBROUTINE DRINT(N,CHAR,BEGIN,ENDE,NUMBER,NUMBRP,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  DRINT
C  PURPOSE:  SUBROUTINE TO RECOGNIZE CHARACTER STRINGS REPRESENTING
C             TWO INTEGERS.  IF THERE ARE ANY ILLEGAL CHARACTERS (NOT
C             INCLUDING BLANKS) ERR = .FALSE. .  OTHERWISE, IF THE
C             STRING DOES REPRESENT TWO INTEGER, ERR = .TRUE., AND
C             NUMBER AND NUMBRP ARE THE VALUE OF THE INTEGERS.
C
      INTEGER             DIGIT
      INTEGER             N
      INTEGER             BLANK,       MINUS,       PLUS,  COMMA
      INTEGER             BEGIN,       CHAR(N),     ENDE,   I
      INTEGER             K,           NUMBER,      SIGN,  NUMBRP
      INTEGER             ICOMMA
      LOGICAL             ERR
      DATA  BLANK,PLUS,MINUS,COMMA / 1H , 1H+, 1H-, 1H, /
      NUMBER = 0
      ERR    = .TRUE.
      SIGN   = 1
           DO 200 I = BEGIN,ENDE
           ICOMMA = I
           IF ( CHAR(I).EQ.ICOMMA) GO TO 200
           IF ( CHAR(I).EQ.BLANK.OR.CHAR(I).EQ.PLUS ) GO TO 200
           IF ( CHAR(I).NE.MINUS ) GO TO 100
           SIGN = -1
           GO TO 200
  100      K = DIGIT(CHAR(I))
           IF ( K.EQ.(-1) ) GO TO 300
           NUMBER = NUMBER*10+K
  200      CONTINUE
      NUMBER = SIGN*NUMBER
      RETURN
  300 ERR    = .FALSE.
      BEGIN = ICOMMA + 1
      CALL SRINT(N,CHAR,BEGIN,ENDE,NUMBRP,ERR)
      RETURN
      END
      DOUBLE PRECISION FUNCTION EVAL(F,P)
C EVALUATE THE FUNCTION F OR A DERIVATIVE AT THE POINT P
C     COMMON BLOCK / CB     /
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      DOUBLE PRECISION    P(1),        DERIV
      DOUBLE PRECISION F
      EXTERNAL F
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      NDIM = IOPDM(2)
      EVAL = DERIV(F,NDIM,P,ICBD,CBD,CBDSW)
      RETURN
      END
      INTEGER FUNCTION FINDCM(M,V)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  FINDCM
C  PURPOSE:  TO FIND A COMMON MULTIPLE FOR THE STEP-SIZE AND WINDOW
C             "MAGNIFICATIONS" M AND V.
C
      INTEGER             MOD
      INTEGER             CM,          IM,          IV,          M
      INTEGER             V
      IM = 1
      IV = 1
      IF ( M.LT.0 ) IM = -M
      IF ( V.LT.0 ) IV = -V
      CM = IM*IV
      IF ( IM.EQ.IV ) CM = IM
      IF ( MOD(IM,IV).EQ.0 ) CM = IM
      IF ( MOD(IV,IM).EQ.0 ) CM = IV
      FINDCM = CM
      RETURN
      END
      INTEGER FUNCTION FINDCO(N,CHAR,BEGIN,END)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  FINDCO
C  PURPOSE:  TO LOCATE THE FIRST OCCURRANCE OF A COMMA IN THE FIELD
C             (BEGIN,END) OF THE ARRAY CHAR().  IF NOT FOUND, FINDCO = 0
C
      INTEGER             N
      INTEGER             COMMA
      INTEGER             BEGIN,       CHAR(N),     END,         I
      DATA  COMMA / 1H, /
      FINDCO = 0
           DO 100 I = BEGIN,END
           IF ( CHAR(I).NE.COMMA ) GO TO 100
           FINDCO = I
           GO TO 200
  100      CONTINUE
  200 RETURN
      END
      INTEGER FUNCTION IFN(N,M)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  IFN
C  PURPOSE:  TO CALCULATE THE RESULTING VALUE N AFTER APPLYING
C             THE USER SPECIFIED "MAGNIFICATION" TERM M
C
      INTEGER             M,           N
      IFN = N*M
      IF ( M.LT.0 ) IFN = -N/M
      RETURN
      END
      SUBROUTINE INHELP
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  INHELP
C  PURPOSE:  TO INPUT PARTS OF THE HELP DATA FILE AND SET UP LOOKUP
C             TABLE (JHELP(99,2)) FOR LATER USE.
C
C  VARIABLES:
C
C     HELP         =  HELP DATA FILE
C     JHELP1       =  THE NUMBER OF LINES IN THE HELP SUMMARY LIST
C     JHELP2       =  THE NUMBER OF LINES IN THE DETAILED HELP LIST
C     JHELP3       =  THE NUMBER OF COMMANDS (= THE NUMBER OF LINES
C                      OF THE INTERACTIVE PROMPT COMMANDS)
C     JHELP(I,J)   =  LOOKUP TABLE
C                       JHELP(I,1)  SPECIFIES THE LINE WHICH THE
C                                    DETAILED HELP DOCUMENTATION ON
C                                    THE I'TH COMMAND BEGINS
C                       JHELP(I,2)  SPECIFIES THE NUMBER OF LINES
C                                    OF DOCUMENTATION FOR THE COMMAND
C     IHELP(I,J)   =  THE ARRAY CONTAINING THE PROMPT COMMANDS.
C                      I SPECIFIES THE I'TH COMMAND.
C
      INTEGER             I,           J,           IBEGIN,      ISKIP
      INTEGER             IDUM,        NLNS
      INTEGER             HELP,        JHELP1,      JHELP2,      JHELP3
      INTEGER             JHELP,       IHELP,       KHELP
      COMMON / HELPER /   HELP,        JHELP1,      JHELP2,      JHELP3
      COMMON / HELPER /   JHELP(99,2), IHELP(72,99)
C
      READ (HELP,10000) JHELP1,JHELP2,JHELP3
      READ (HELP,20000) (JHELP(I,2),I=1,JHELP3)
      NLNS = JHELP1+JHELP2
      NL = 6+NLNS
      CALL POS(HELP,NL)
           DO 200 I = 1,JHELP3
           READ (HELP,30000) (IHELP(J,I),J=1,72)
  200      CONTINUE
      IBEGIN = JHELP1+1
           DO 300 I = 1,JHELP3
           JHELP(I,1) = IBEGIN
           IBEGIN     = IBEGIN + JHELP(I,2)
           JHELP(I,2) = IBEGIN - 1
  300      CONTINUE
      RETURN
C
C  FORMAT STATEMENTS
C
10000 FORMAT(3I3)
20000 FORMAT(20I2)
30000 FORMAT(72A1)
      END
      INTEGER FUNCTION LENGTH(N,CHAR)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  LENGTH
C  PURPOSE:  TO DETERMINE THE LENGTH OF A CHARACTER STRING EMBEDDED IN
C             THE INTEGER ARRAY CHAR(N).
C
      INTEGER             N
      INTEGER             BLANK
      INTEGER             CHAR(N),     I,           M,           NP
      DATA  BLANK / 1H  /
      LENGTH = 0
      NP     = N+1
           DO 100 I = 1,N
           M = NP-I
           IF ( CHAR(M).EQ.BLANK ) GO TO 100
           LENGTH = M
           GO TO 200
  100      CONTINUE
  200 RETURN
      END
      SUBROUTINE MAKROM
C IF CROSS DERIVATIVES ARE BEING PLOTTED PREY ON THE GRAPHICS PART
C OF THE DISPLAY TO OBTAIN NUMERICAL OUTPUT
      INTEGER             IA
      INTEGER             I,           J
      LOGICAL             TOGGLE
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             ILPUSR,      IDSUSR
      COMMON / ROOM   /   ILPUSR,      IDSUSR
      DATA TOGGLE /.FALSE./
      DATA IA /1HA/
      IF (IDSPLA.GE.IPRMPT+7) GO TO 600
      IF (ICBD.EQ.0) GO TO 300
           DO 100 I = 1,7
           IF (.NOT.OPDF1(I,2)) GO TO 300
  100      CONTINUE
      TOGGLE = .TRUE.
           DO 200 I = 1,7
           LPLT(I) = .FALSE.
  200      CONTINUE
      ILP = ILPUSR - 1
      IDSPLA = IDSUSR + 1
      GO TO 600
  300 CONTINUE
      IF (.NOT.TOGGLE) GO TO 600
      TOGGLE = .FALSE.
           DO 400 I = 1,7
           LPLT(I) = .FALSE.
  400      CONTINUE
      ILP = ILPUSR
      IDSPLA = IDSUSR
           DO 500 J = 1,WIDTH
           ISCRN1(J,ILP) = IA
  500      CONTINUE
  600 CONTINUE
      RETURN
      END
      SUBROUTINE MAKSHF
C MAKE THE SHIFTED POINT THE NEW POINT
      INTEGER             I,           NDIM
      DOUBLE PRECISION    DNH,         SH,          FAC
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             M,           SHIFT,       V
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      IF (IOPSH(2).EQ.0) GO TO 200
      NDIM = IOPDM(2)
      DNH = NH(2)
      SH  = IOPSH(2)
      FAC = DNH*SH*ROPSTS(2)/8.0D0
           DO 100 I = 1,NDIM
           ROPPNT(I,2) = ROPPNT(I,2) + FAC*ROPUDI(I,2)
  100      CONTINUE
  200 CONTINUE
      RETURN
      END
      SUBROUTINE MINMAX
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  MINMAX
C  PURPOSE:  TO DETERMINE THE MINIMUM/MAXIMUM OF EACH DERIVATIVE TO BE
C             PLOTTED ON THE INDEX INTERVAL (ILEFT,IRIGHT).  THE DFMNMX
C             ARRAY IS THEN UPDATED AND A NEW SCALE FACTOR IS CALCULATED
C
      DOUBLE PRECISION    DMAX1,       DMIN1
      REAL                SNGL
      DOUBLE PRECISION    DILP,        FMAX,        FMIN
      INTEGER             I,           IDF
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
C
C  EVALUATE THE MINIMUM/MAXIMUM'S ONLY IF THE CORRESPONDING DERIVATIVES
C   HAVE BEEN CALCULATED AND THE DERIVATIVE IS TO BE PLOTTED.
C
      DILP = ILP
      DILP = DILP-0.1D0
           DO 200 IDF = 1,7
           IF ( LPLT(IDF) ) GO TO 200
           IF ( .NOT.(OPDF1(IDF,2) .OR. OPDF2(IDF,2)) ) GO TO 200
           IF ( .NOT.LDF(IDF) ) GO TO 200
C
C  DETERMINE THE MINIMUM/MAXIMUM OVER THE SPECIFIED INTERVAL
C
           FMIN = 1.0D30
           FMAX = -1.0D30
                DO 100 I = 1,WIDTH
                FMIN = DMIN1(FMIN,DF(I,IDF))
                FMAX = DMAX1(FMAX,DF(I,IDF))
  100           CONTINUE
           DFMNMX(1,IDF) = FMIN
           DFMNMX(2,IDF) = FMAX
           SCALE(IDF) = 0.0
           IF ( FMAX-FMIN.GT.0.0D0 ) SCALE(IDF) = SNGL(DILP/(FMAX-FMIN))
  200      CONTINUE
      RETURN
      END
      SUBROUTINE NRML
C TO TOGGLE THE FLAG INDICATING WHETHER THE DIRECTION OF DIFFERENTIATION
C TO BE NORMALIZED
      DOUBLE PRECISION    DSQRT
      DOUBLE PRECISION    SUM
      INTEGER             I
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
      EQUIVALENCE (NDIM,IOPDM(2))
      SUM = 0.0D0
           DO 100 I =1,NDIM
           SUM = SUM + ROPDI(I,2)**2
  100      CONTINUE
      IF (SUM.NE.0.0D0) GO TO 200
      ERR = .TRUE.
      ERRCOD = 25
      GO TO 600
  200 CONTINUE
      IF (NORMAL) GO TO 400
           DO 300 I = 1,NDIM
           ROPUDI(I,2) = ROPDI(I,2)
  300      CONTINUE
      GO TO 600
  400 CONTINUE
      SUM = DSQRT(SUM)
           DO 500 I = 1,NDIM
           ROPUDI(I,2) = ROPDI(I,2)/SUM
  500      CONTINUE
  600 CONTINUE
      RETURN
      END
      SUBROUTINE NRMLC
C PURPOSE: NORMALIZE THE CROSS DIRECTION IF REQUIRED
      DOUBLE PRECISION    DSQRT
      DOUBLE PRECISION    SUM
      INTEGER             I
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
      EQUIVALENCE (NDIM,IOPDM(2))
      SUM = 0.0D0
           DO 100 I = 1,NDIM
           SUM = SUM + CBDU(I)**2
  100      CONTINUE
      IF (SUM.EQ.0.0D0) GO TO 600
      IF (.NOT.NORMAL) GO TO 300
      SUM = DSQRT(SUM)
           DO 200 I = 1,NDIM
           CBD(I) = CBDU(I)/SUM
  200      CONTINUE
      GO TO 500
  300 CONTINUE
      DO 350 I = 1,NDIM
      CBD(I) = CBDU(I)
350   CONTINUE
      GO TO 500
  600 CONTINUE
      ERRCOD = 18
      ERR = .TRUE.
      ICBD = 0
500   CONTINUE
      RETURN
      END
      SUBROUTINE NUMDIG(ND)
C PURPOSE: DETERMINE THE MAXIMUM NUMBER OF SIGNIFICANT DIGITS
      DOUBLE PRECISION    DLOG10
      DOUBLE PRECISION    X,           OX
      INTEGER             ND
      X = 1.0D0
  100 CONTINUE
      X = X/2.0D0
      OX = 1.0D0 + X
      IF (OX.GT.1.0D0) GO TO 100
      X = X + X
      OX = -DLOG10(X)
      ND = OX
      ND = ND + 1
      RETURN
      END
      SUBROUTINE OPCOPY
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  OPCOPY
C  PURPOSE:  TO COPY THE NEW SET OF OPTIONS OVER THE OLD SET AND RESET
C             FLAGS.
C
      INTEGER             I
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
C
C  COPY BACK ALL VECTORS
C
           DO 100 I = 1,3
           ROPDI(I,1)   = ROPDI(I,2)
           ROPPNT(I,1)  = ROPPNT(I,2)
           ROPDR1(I,1)  = ROPDR1(I,2)
           ROPDR2(I,1)  = ROPDR2(I,2)
           ROPUDI(I,1)  = ROPUDI(I,2)
  100      CONTINUE
      ROPSTS(1)    = ROPSTS(2)
      ROPSTW(1)    = ROPSTW(2)
      NH(1)        = NH(2)
      NW(1)        = NW(2)
      ILEFT(1)     = ILEFT(2)
      IRIGHT(1)    = IRIGHT(2)
      IOPDM(1)     = IOPDM(2)
      IOPSH(1)     = IOPSH(2)
      IOPSS(1)     = 1
      IOPWN(1)     = 1
      IOPSS(2)     = 1
      IOPWN(2)     = 1
      OPDC(1)      = OPDC(2)
      OPDS(1)      = OPDS(2)
      OPDX(1)      = OPDX(2)
      OPCMP(1)     = OPCMP(2)
      OPSPL(1)     = OPSPL(2)
           DO 200 I = 1,7
           OPDF1(I,1)  = OPDF1(I,2)
           OPDF2(I,1)  = OPDF2(I,2)
  200      CONTINUE
      OPCMP(2)     = .FALSE.
      RETURN
      END
      INTEGER FUNCTION PWROF2(M)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  PWROF2
C  PURPOSE:  GIVEN AN INTEGER M, PWROF2 TAKES THE ABSOLUTE VALUE OF M
C             AND DETERMINES WHETHER OR NOT IT IS AN INTEGER POWER OF
C             2.  IF NOT, PWROF2 = 0.  IF IT IS A POWER OF 2, THEN
C             PWROF2 RETURNS THE EXPONENT.
C
      INTEGER             ISIGN,       MOD
      INTEGER             K,           L(10),       M,           MP
      DATA L(1),L(2),L(3),L(4),L(5) /  2,   4,   8,  16,  32 /
      DATA L(6),L(7),L(8),L(9),L(10)/ 64, 128, 256, 512,1024 /
      MP = M
      IF ( MOD(MP,2).NE.0 ) GO TO 300
      MP = ISIGN(MP,MP)
           DO 200 K = 1,10
           IF ( MP.NE.L(K) ) GO TO 200
           PWROF2 = K
           RETURN
  200      CONTINUE
  300 PWROF2 = 0
      RETURN
      END
      LOGICAL FUNCTION S0(NH,NW,IOPSS,IOPWN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  S0
C  PURPOSE:  TO CHECK IF THE NEW VALUES OF NH AND NW ARE NON-ZERO
C
      INTEGER             IFN
      INTEGER             IOPSS(2),    IOPWN(2),    NH(2),       NHP
      INTEGER             NW(2),       NWP
      S0    = .TRUE.
      NHP   = IFN(NH(1),IOPSS(2))
      NWP   = IFN(NW(1),IOPWN(2))
      IF ( NHP.EQ.0 .OR. NWP.EQ.0 ) S0 = .FALSE.
      RETURN
      END
      LOGICAL FUNCTION S1(NH,NW,IOPSS,IOPWN,WIDTH)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  S1
C  PURPOSE:  TO CHECK IF THE GIVEN MULTIPLICATION/DIVISION FACTORS
C             (IOPSS(2) AND IOPWN(2)) FOR THE STEP-SIZE AND STENCIL HALF
C             WIDTH (NH(2) AND NW(2)) RESULT IN NEW VALUES WHICH CAN BE
C             ACHIEVED WITHOUT CHANGING THE BASIC STEP SIZE H = ROPSTS
C
      INTEGER             IFN
      INTEGER             MOD
      INTEGER             ICENTR,      IOPSS(2),    IOPWN(2),    ISUM
      INTEGER             IXSIZE,      NH(2),       NHP,         NW(2)
      INTEGER             NWP,         WIDTH
      DATA ICENTR,IXSIZE / 2689, 5377 /
      S1   = .TRUE.
      NHP  = IFN(NH(1),IOPSS(2))
      NWP  = IFN(NW(1),IOPWN(2))
      ISUM = NHP*WIDTH+2*NWP
      IF ( MOD(NHP,4).NE.0 .OR. MOD(NWP,4).NE.0 .OR. ISUM.GE.IXSIZE )
     X     S1 = .FALSE.
      RETURN
      END
      SUBROUTINE S2(A,B)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  S2
C  PURPOSE:  TO DETERMINE THE BOUNDS, A AND B, ON THE POSSIBLE VALUES
C             OF GAMMA.
C
C             NOTE:            OLD VALUE OF H
C                    GAMMA  =  --------------
C                              NEW VALUE OF H
C
C                  IT IS USED WHEN THE ONLY WAY TO ACCOMODATE THE
C                  USER SPECIFIED CHANGES IN STEP-SIZE AND WINDOW SIZE
C                  IS TO CHANGE THE BASIC STEP-SIZE H = ROPSTS(2).
C
      DOUBLE PRECISION    DFN
      DOUBLE PRECISION    DMAX1,       DMIN1
      DOUBLE PRECISION    A,           B,           DNHP,        DNWP
      DOUBLE PRECISION    S,           W
      INTEGER             ICENTR,      IXSIZE
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DATA ICENTR,IXSIZE / 2689, 5377 /
      W    = WIDTH
      S    = IXSIZE
      DNHP = DFN(NH(1),IOPSS(2))
      DNWP = DFN(NW(1),IOPWN(2))
      A    = DMAX1(4.0D0/DNHP,4.0D0/DNWP)
      B    = DMIN1(S/(DNHP*W),S/(2.0D0*DNWP))
      RETURN
      END
      SUBROUTINE SAMPLE
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SAMPLE
C  PURPOSE:  TO SET UP THE SAMPLING ARRAY XS() AND RESET THE FS() AND
C             LDEF() ARRAYS.
C
      INTEGER             MOD
      DOUBLE PRECISION    DK,          SPACNG(8)
      INTEGER             I,           ICENTR,      INDX,        IXSIZE
      INTEGER             J
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      DOUBLE PRECISION    H
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE  (H,ROPSTS(2))
      DATA  ICENTR,IXSIZE/ 2689, 5377 /
C
C  SET UP THE SPACING ARRAY (SPACNG) WHICH INCORPORATES THE IRREGULAR
C   SPACING USED
C
      SPACNG(1) = 0.0D0
      SPACNG(2) = H/6.0D0
      SPACNG(3) = H/4.0D0
      SPACNG(4) = H/3.0D0
      SPACNG(5) = H/2.0D0
      SPACNG(6) = 2.0D0*H/3.0D0
      SPACNG(7) = 3.0D0*H/4.0D0
      SPACNG(8) = 5.0D0*H/6.0D0
C
C  SET UP THE SAMPLING ARRAY (XS) AND RESET FS() AND LDEF().
C
           DO 100 I = ICENTR,IXSIZE
           DK        = (I-ICENTR)/8
           INDX      = 2*ICENTR-I
           J         = MOD(I-ICENTR,8)+1
           XS(I)     = DK*H+SPACNG(J)
           XS(INDX)  = -XS(I)
  100      CONTINUE
      RETURN
      END
      SUBROUTINE SCMUPD(F,NCALLS)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SCMUPD
C  PURPOSE:  TO CALCULATE THE DIRECTIONAL DERIVATIVES OF THE USER
C             SUPPLIED FUNCTION F() UP TO THE 6TH DERIVATIVE.
C
      DOUBLE PRECISION    EVAL
      DOUBLE PRECISION F
      EXTERNAL F
      INTEGER             MOD
      DOUBLE PRECISION    D,           FP(9),       G0,          G1
      DOUBLE PRECISION    G2,          G3,          G4,          G5
      DOUBLE PRECISION    G6,          H1,          H2,          H3
      DOUBLE PRECISION    H4,          H5,          H6,          HP
      DOUBLE PRECISION    PS(3)
      INTEGER             I,           ICENTR,      IP,          ITEMP
      INTEGER             IXSIZE,      J,           JINDEX(9),   JP
      INTEGER             JS,          JTEMP,       K,           MODE
      INTEGER             NCALLS
      DOUBLE PRECISION    H
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             DMNSN,       M,           SHIFT
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
C     COMMON BLOCK / CB     /
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      EQUIVALENCE (H,ROPSTS(2)),(M,IOPSS(2)),(SHIFT,IOPSH(2))
      EQUIVALENCE (DMNSN,IOPDM(2))
      DATA ICENTR,IXSIZE / 2689, 5377 /
C
C  SET THE INTERNAL MODE DEPENDING ON WHAT DERIVATIVES ARE TO BE
C   CALCULATED.
C
      MODE = 0
      IF ( OPDF1(1,2) ) MODE = 1
      IF ( OPDF1(2,2).OR.OPDF1(3,2) ) MODE = 2
      IF ( OPDF1(4,2).OR.OPDF1(5,2) ) MODE = 3
      IF ( OPDF1(6,2).OR.OPDF1(7,2) ) MODE = 4
C
C  IF NO FUNCTION AT ALL IS TO BE PLOTTED THEN RETURN
C
      IF ( MODE.EQ.0 ) GO TO 1100
C
C  SET THE LOGICAL VALUES OF THE LDF() AND LPLT() ARRAYS
C
      LDF(1)   = .TRUE.
      LPLT(1)  = .FALSE.
      IF ( MODE.EQ.1 ) GO TO 100
      LDF(2)   = .TRUE.
      LDF(3)   = .TRUE.
      LPLT(2)  = .FALSE.
      LPLT(3)  = .FALSE.
      IF ( MODE.EQ.2 ) GO TO 100
      LDF(4)   = .TRUE.
      LDF(5)   = .TRUE.
      LPLT(4)  = .FALSE.
      LPLT(5)  = .FALSE.
      IF ( MODE.EQ.3 ) GO TO 100
      LDF(6)   = .TRUE.
      LDF(7)   = .TRUE.
      LPLT(6)  = .FALSE.
      LPLT(7)  = .FALSE.
C
C  CALCULATE THE VARIOUS POWERS OF THE SAMPLING STEP (HP) USED IN
C   CALCULATING THE DERIVATIVES.
C
C  SKIP THE CALCULATION OF UNNEEDED POWERS AS THEY MAY LEAD TO
C   FLOATING POINT OVERFLOWS OTHERWISE
  100 CONTINUE
      D  = NW(2)
      HP = D*H/4.0D0
      IF (MODE.EQ.1) GO TO 200
      H1 = 1.0D0/HP
      H2 = (2.0D0/HP)**2
      IF (MODE.EQ.2) GO TO 200
      H3 = (4.0D0/HP)**3/4.0D0
      H4 = (4.0D0/HP)**4
      IF (MODE.EQ.3) GO TO 200
      H5 = (6.0D0/HP)**5/2.0D0
      H6 = (6.0D0/HP)**6
  200 CONTINUE
C
C  SET UP THE INDEXING ARRAY  JINDEX()  WHICH IS USED TO INDEX THE
C   XS() AND FS() ARRAYS FOR PROPER EVALUATION OF THE DERIVATIVES.
C
      JTEMP = (NW(2)+1)/3
      JINDEX(5) = 0
      JINDEX(6) = JTEMP
      JINDEX(7) = NW(2)/2
      JINDEX(8) = NW(2)-JTEMP
      JINDEX(9) = NW(2)
      JINDEX(1) = -JINDEX(9)
      JINDEX(2) = -JINDEX(8)
      JINDEX(3) = -JINDEX(7)
      JINDEX(4) = -JINDEX(6)
      IP        = ICENTR+SHIFT*NH(1)
           DO 300 I = 1,9
           JINDEX(I) = JINDEX(I)+IP
  300      CONTINUE
C
C  COMPUTE INDEXING OF THE XS() AND FS() ARRAYS AND CHECK TO SEE IF F()
C   HAS ALREADY BEEN EVALUATED AT THE GIVEN ELEMENTS.  THEN CALCULATE
C    THE DERIVATIVES.
C
      ITEMP = WIDTH/2+1
           DO 1000 I = 1,WIDTH
           JP = (I-ITEMP)*NH(2)
C
C  EVALUATE THE INTERPOLATION FUNCTION F()
C
                DO 800 J = 1,9
                IF ( MODE.EQ.1.AND.J.NE.5 ) GO TO 800
                IF ( J.NE.1.AND.J.NE.5.AND.J.NE.9.AND.MODE.EQ.2 ) GO TO
     X               800
                IF ( MOD(J,2).EQ.0.AND.MODE.EQ.3 ) GO TO 800
                JS   = JINDEX(J)+JP
                IF ( .NOT. LDEF(JS) ) GO TO 500
                FP(J)= FS(JS)
                GO TO 800
  500           D    = XS(JS)
                     DO 600 K = 1,DMNSN
                     PS(K) = ROPPNT(K,2)+D*ROPUDI(K,2)
  600                CONTINUE
                FP(J)   = EVAL(F,PS)
                NCALLS  = NCALLS+ICBD+1
                FS(JS)  = FP(J)
                LDEF(JS)= .TRUE.
  800           CONTINUE
C
C  CALCULATE THE DERIVATIVES CORRESPONDING TO THE MODE VARIABLE
C
           DF(I,1) = FP(5)
           IF ( MODE.EQ.1 ) GO TO 1000
           G0      = FP(5)
           G1      = FP(1)+FP(9)
           G2      = -(FP(1)-FP(9))
           DF(I,2) = G2*H1
           DF(I,3) = (G1-2.0D0*G0)*H2
           IF ( MODE.EQ.2 ) GO TO 1000
           G3      = 4.0D0*(FP(3)+FP(7))
           G4      = 4.0D0*(FP(3)-FP(7))
           DF(I,4) = (2.0D0*G2+G4)*H3
           DF(I,5) = (G1-G3+6.0D0*G0)*H4
           IF ( MODE.EQ.3 ) GO TO 1000
           G5      =  4.0D0*(FP(2)-FP(8))+ 5.0D0*(FP(6)-FP(4))
           G6      = -6.0D0*(FP(2)+FP(8))+15.0D0*(FP(4)+FP(6))
           DF(I,6) = (G2+G5)*H5
           DF(I,7) = (G1+G6-20.0D0*G0)*H6
 1000      CONTINUE
      OPCMP(2) = .TRUE.
      CALL SFBNDS(ILEFT(2),IRIGHT(2))
 1100 CONTINUE
      RETURN
      END
      SUBROUTINE SCROLL(DEVICE,NCALLS,LGO)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SCROLL
C  PURPOSE:  THE SIMPLEST DISPLAY METHOD (I.E. SCROLLING).  THE ROUTINE
C             WRITES THE PLOT ARRAY ISCRN2() AND THE DISPLAYED DATA
C             ONTO UNIT = DEVICE.
C
      LOGICAL LGO
      INTEGER             DEVICE,      I,           J,           MODE
      INTEGER             NCALLS
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
      INTEGER             OUTPTD,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPTD,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
C
C  WRITE THE ISCRN2() ARRAY ONTO UNIT DEVICE
C
           DO 100 I = 1,ILP
           WRITE (DEVICE,10000) (ISCRN2(J,I),J=1,WIDTH)
  100      CONTINUE
      MODE = 1
      CALL SDDATA(DEVICE,MODE,NCALLS,LGO)
      RETURN
10000 FORMAT(1H ,135A1)
      END
      SUBROUTINE  SFBNDS(LEFT,RIGHT)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SFBNDS
C  PURPOSE:  TO FIND THE LEFT AND RIGHT EXTENT OF THE CALCULATED VALUES
C             IN THE FS() ARRAY.
C
      INTEGER             I,           ICENTR,      IP,          IXSIZE
      INTEGER             K,           LEFT,        RIGHT
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DATA ICENTR,IXSIZE / 2689, 5377 /
      LEFT   = IXSIZE+1
      RIGHT  = 0
           DO 100 I = 1,IXSIZE
           IF ( .NOT.LDEF(I) ) GO TO 100
           LEFT = I
           GO TO 200
  100      CONTINUE
  200 IP     = IXSIZE+1
           DO 300 I = 1,IXSIZE
           K     = IP-I
           IF ( .NOT.LDEF(K) ) GO TO 300
           RIGHT = K
           GO TO 400
  300      CONTINUE
  400 RETURN
      END
      SUBROUTINE SGAMMA(LGO,NUM,DEN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SGAMMA
C  PURPOSE:  THIS ROUTINE CHECKS THE NEW SET OF OPTIONS BEFORE LETTING
C             COMPUTATION BEGIN.  ITS PRIMARY PURPOSE IS TO CALCULATE
C
C                                        NUM         OLD STEP-SIZE
C                           GAMMA   =   -----   =   ---------------
C                                        DEN         NEW STEP-SIZE
C
      INTEGER             FINDCM,      PWROF2
      LOGICAL             S0,          S1
      INTEGER             MOD
      DOUBLE PRECISION    A,           B,           DI,          DJ
      DOUBLE PRECISION    DNW,         GAMMA
      INTEGER             BOTTOM(100), DEN,         I,           I1
      INTEGER             I2,          IBEST,       ICENTR,      IHIGH
      INTEGER             INDX,        ISUM,        ITYPE,       IXSIZE
      INTEGER             J,           JP,          MP,          NHB
      INTEGER             NHP,         NUM,         NWB,         NWP
      INTEGER             P,           RANK(100),   TOP(100),    VP
      INTEGER             W2,          W3
      LOGICAL             LGO
      REAL                D,           DDEN,        DNUM
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      DOUBLE PRECISION    H
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             M,           V
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPTD,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPTD,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      EQUIVALENCE (H,ROPSTS(2)),(M,IOPSS(2)),(V,IOPWN(2))
      DATA ICENTR,IXSIZE / 2689, 5377 /
      LGO   = .TRUE.
      W3    = 8*WIDTH
C
C  USE THE LOGICAL FUNCTIONS S0 AND S1 TO SEE IF THE EXISTING STEP-SIZE
C   H = ROPSTS(2) CAN BE USED FOR THE SPECIFIED VALUES OF M=IOPSS(2) AND
C   V = IOPWN(2).
C
      IF ( .NOT.S0(NH,NW,IOPSS,IOPWN) ) GO TO 500
      IF ( .NOT.S1(NH,NW,IOPSS,IOPWN,WIDTH) ) GO TO 500
      NH(2) = NH(1)*M
      IF ( M.LT.0 ) NH(2) = -NH(1)/M
      NW(2) = NW(1)*V
      IF ( V.LT.0 ) NW(2) = -NW(1)/V
      NUM = 1
      DEN = 1
      DNW = NW(2)
      ROPSTW(2) = ROPSTS(2)*DNW/4.0D0
      GO TO 1700
C
C  BY GETTING TO STATEMENT 1000 IT IS NOW ESTABLISHED THAT THE STEP-SIZE
C   H MUST BE CHANGED TO ACCOMODATE THE USER REQUESTS.  THE VARIABLE
C   "GAMMA" DENOTES THE RATIO OF THE OLD-SPACING (H) TO THE NEW-SPACING
C   (WHICH IS NOT YET KNOWN).  IT CAN BE SHOWN THAT GAMMA MUST LIE
C   BETWEEN THE BOUNDS A AND B CALCULATED IN THE ROUTINE S2
C
  500 CALL S2(A,B)
      IF ( A.LT.B ) GO TO 600
      ERRCOD = 8
      ERR    = .TRUE.
      ITYPE  = 1
      GO TO 1700
C
C  IF BOTH THE STEP-SIZE AND THE STENCIL-WIDTH ARE BEING SCALED BY
C   THE SAME AMOUNT, THEN UPDATE ROPSTS(2), NUM, DEN, AND RETURN.
C
  600 IF ( M.NE.V ) GO TO 800
      IF ( M.LT.0 ) GO TO 700
      NUM       = 1
      DEN       = M
      D         = M
      ROPSTS(2) = D*ROPSTS(2)
      DNW       = NW(2)
      ROPSTW(2) = ROPSTS(2)*DNW/4.0D0
      GO TO 1700
  700 NUM       = -M
      DEN       = 1
      D         = -M
      ROPSTS(2) = ROPSTS(2)/D
      DNW       = NW(2)
      ROPSTW(2) = ROPSTS(2)*DNW/4.0D0
      GO TO 1700
  800 INDX = 0
      NHB  = NH(1)/4
      NWB  = NW(1)/4
      P    = FINDCM(M,V)
      MP   = M*P
      IF ( M.LT.0 ) MP = -P/M
      VP   = V*P
      IF ( V.LT.0 ) VP = -P/V
C
C  FIND POSSIBLE CHOICES FOR THE NUMERATOR AND DENOMINATOR OF GAMMA
C
           DO 1100 I = 1,32
           I1 = I*MP*NHB
           I2 = I*VP*NWB
           DI = I
           IF ( MOD(I1,P).NE.0 .OR. MOD(I2,P).NE.0 ) GO TO 1100
                DO 1000 J = 1,32
                JP = J*P
                DJ = J
                IF ( MOD(I1,JP).NE.0 .OR. MOD(I2,JP).NE.0 ) GO TO 1000
                GAMMA = DI/DJ
                IF ( .NOT.(A.LE.GAMMA .AND. GAMMA.LT.B) ) GO TO 1000
                INDX         = INDX+1
                TOP(INDX)    = I
                BOTTOM(INDX) = J
                IF ( INDX.EQ.100 ) GO TO 1200
 1000           CONTINUE
 1100      CONTINUE
 1200 IF ( INDX.GT.0 ) GO TO 1300
      ERRCOD = 9
      ERR    = .TRUE.
      ITYPE  = 1
      GO TO 1700
C
C  NEXT, LOOK AT ALL THE POSSIBLE COMBINATIONS OF NUMERATORS AND
C   DENOMINATORS AND RANK THEM IF THEY DO NOT BLOW THE ARRAY BOUNDS.
C
 1300      DO 1400 I = 1,INDX
           NHP     = TOP(I)*MP*NH(1)/(P*BOTTOM(I))
           NWP     = TOP(I)*VP*NW(1)/(P*BOTTOM(I))
           ISUM    = NHP*WIDTH+2*NWP
           RANK(I) = 0
           IF ( ISUM.GT.IXSIZE ) GO TO 1400
           RANK(I) = RANK(I)+1
           IF ( W3.LE.ISUM .AND. ISUM.LT.ICENTR ) RANK(I) = RANK(I)+1
           IF ( TOP(I).EQ.1 .AND. PWROF2(BOTTOM(I)).NE.0 ) RANK(I) =
     X          RANK(I)+1
           IF ( PWROF2(TOP(I)).NE.0 .AND. BOTTOM(I).EQ.1 ) RANK(I) =
     X          RANK(I)+1
           IF ( NHP.GT.4 .AND. NWP.GT.4 ) RANK(I) = RANK(I)+1
 1400      CONTINUE
C
C  FIND THE PAIR WITH THE HIGHEST RANK
C
      IBEST = 1
      IHIGH = RANK(1)
           DO 1500 I = 2,INDX
           IF ( RANK(I).LE.IHIGH ) GO TO 1500
           IBEST = I
           IHIGH = RANK(I)
 1500      CONTINUE
      NUM       = TOP(IBEST)
      DEN       = BOTTOM(IBEST)
      DNUM      = NUM
      DDEN      = DEN
      ROPSTS(2) = ROPSTS(1)*DDEN/DNUM
      NH(2)     = NUM*MP*NH(1)/(DEN*P)
      NW(2)     = NUM*VP*NW(1)/(DEN*P)
      DNW       = NW(2)
      ROPSTW(2) = ROPSTS(2)*DNW/4.0D0
      IF ( IHIGH.GT.0 ) GO TO 1700
      ERRCOD = 10
      ERR    = .TRUE.
      ITYPE  = 1
 1700 RETURN
      END
      SUBROUTINE SGRUPD(ICOM)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SGRUPD
C  PURPOSE:  TO SET UP THE GRAPHICS DATA, I.E. ISCRN2(), WHICH THE
C             SGRAPH ROUTINE ACCESSES FOR PLOTTING.
C
      INTEGER             IABS,        IDINT,       MOD
      INTEGER             AST,         BLANK,       COL,         DASH
      INTEGER             IY,          NUM(10),     PLUS
      DOUBLE PRECISION    RMAX,        SCAL
      INTEGER             CHAR,        I,           ICOM,        IDF
      INTEGER             ILPD2,       J,           JCENTR,      JLEFT
      INTEGER             JRIGHT,      WINDOW
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DATA BLANK,AST,COL,DASH,PLUS,IY /1H , 1H*, 1H:, 1H-, 1H+, 1HI/
      DATA NUM(1),NUM(2),NUM(3),NUM(4),NUM(5)/ 1H., 1H1, 1H2, 1H3, 1H4 /
      DATA NUM(6),NUM(7),NUM(8),NUM(9),NUM(10)/1H5, 1H6, 1H7, 1H8, 1H9 /
C
C  IF ICOM = 52   (I.E. REFRESH SCREEN WITHOUT UPDATING) THEN COPY THE
C   ISCRN1() ARRAY ONTO THE ISCRN2() ARRAY AND SKIP THE UPDATING.
C
      IF ( ICOM.NE.52 ) GO TO 200
           DO 100 I = 1,WIDTH
                DO 100 J = 1,ILP
                ISCRN2(I,J) = ISCRN1(I,J)
  100           CONTINUE
      RETURN
C
C  SET UP THE WINDOW AND X-AXIS DEPENDING ON THE LOGICAL FLAGS, ETC.
C
  200 WINDOW = NW(2)/NH(2)
      JCENTR = WIDTH/2+1
      JRIGHT = JCENTR+WINDOW
      JLEFT  = JCENTR-WINDOW
C
C  FIRST, CLEAR THE ISCRN2() ARRAY
C
           DO 300 I = 1,WIDTH
                DO 300 J = 1,ILP
                ISCRN2(I,J) = BLANK
  300           CONTINUE
           DO 400 I = 1,ILP
           ISCRN2(JCENTR,I) = COL
  400      CONTINUE
      IF (JRIGHT.GT.135.OR.JLEFT.LT.1) GO TO 600
           DO 500 I = 1,ILP
           ISCRN2(JRIGHT,I) = IY
           ISCRN2(JLEFT,I)  = IY
  500      CONTINUE
  600 CONTINUE
      ILPD2  = ILP/2+1
      IF ( .NOT.OPDX(2) ) GO TO 800
           DO 700 I = 1,WIDTH
           ISCRN2(I,ILPD2)  = DASH
  700      CONTINUE
  800 IF ( .NOT.OPDS(2) ) GO TO 1000
           DO 900 I = 1,WIDTH
           J               = IABS(I-JCENTR)
           J               = MOD(J,10)+1
           ISCRN2(I,ILPD2) = NUM(J)
  900      CONTINUE
C
C  DETERMINE THE SCALE FACTORS WHICH ARE USED TO COMPRESS/EXPAND
C   THE GRAPHS OF THE FUNCTION AND ITS DERIVATIVES SO THAT THEY TAKE
C   THE ENTIRE VERTICAL EXTENT OF THE GRAPH.
C
 1000 CALL MINMAX
C
C  APPLY THE SCALE FACTORS TO THE DATA WHICH IS TO BE DISPLAYED
C
           DO 1400 IDF = 1,7
           IF ( LPLT(IDF).OR..NOT.LDF(IDF) ) GO TO 1400
           IF ( .NOT.OPDF1(IDF,2).AND..NOT.OPDF2(IDF,2) ) GO TO 1400
           SCAL = SCALE(IDF)
           RMAX = DFMNMX(2,IDF)
           IF ( SCAL.NE.0.0D0 ) GO TO 1200
                DO 1100 J = 1,WIDTH
                IPLOT(J,IDF) = ILPD2
 1100           CONTINUE
           GO TO 1400
 1200           DO 1300 J = 1,WIDTH
                IPLOT(J,IDF) = 1+IDINT(SCAL*(RMAX-DF(J,IDF)))
 1300           CONTINUE
           LPLT(IDF) = .TRUE.
 1400      CONTINUE
C
C  COPY THE IPLOT() ARRAY ONTO THE ISCRN2() ARRAY FOR THE DERIVATIVES
C   WHICH THE USER WANTS TO PLOT.
C
           DO 1600 IDF = 1,7
           IF ( OPDF2(IDF,2) ) GO TO 1600
           IF ( .NOT.( LDF(IDF).AND.OPDF1(IDF,2) ) ) GO TO 1600
           CHAR = NUM(IDF)
                DO 1500 I = 1,WIDTH
                J           = IPLOT(I,IDF)
                ISCRN2(I,J) = CHAR
 1500           CONTINUE
 1600      CONTINUE
C
C  ACCENTUATE THE DERIVATIVES FOR WHICH OPDF2(N,2) = .TRUE.
C
           DO 1800 IDF = 1,7
           IF ( .NOT.( LDF(IDF).AND.OPDF2(IDF,2) ) ) GO TO 1800
                DO 1700 I = 1,WIDTH
                J           = IPLOT(I,IDF)
                ISCRN2(I,J) = AST
 1700           CONTINUE
 1800      CONTINUE
C
C  PLACE A MARK AT THE CENTER IF OPDC(2) = .TRUE.
C
      IF ( OPDC(2) ) ISCRN2(JCENTR,ILPD2) = PLUS
      RETURN
      END
      SUBROUTINE SIREAD(DEVICE,IARG,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SIREAD
C  PURPOSE:  TO READ IN A SINGLE INTEGER NUMBER, FREE FORMATTED,
C             FROM THE SPECIFIED DEVICE.
C
      INTEGER             LENGTH
      INTEGER             BEGIN,       CHAR(72),    DEVICE,      ENDE
      INTEGER             I,           I72,         IARG
      LOGICAL             ERR
      DATA  I72 / 72 /
C
C  SET ERR = .TRUE. (I.E. NO ERRORS)
C
      ERR  = .TRUE.
C
C  READ CHARACTER STRING FROM INPUT DEVICE AND RECOGNIZE NUMBER
C
      IARG = 0
      READ (DEVICE,10000) (CHAR(I),I=1,72)
      ENDE   = LENGTH(I72,CHAR)
      IF ( ENDE.EQ.0 ) GO TO 200
      BEGIN= 1
      CALL SRINT(I72,CHAR,BEGIN,ENDE,IARG,ERR)
  200 RETURN
C
C  FORMAT STATEMENT
C
10000 FORMAT(72A1)
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SRDP
C  PURPOSE:  SUBROUTINE TO RECOGNIZE A CHARACTER STRING REPRESENTING
C             A DOUBLE PRECISION VARIABLE (THE VALUE IS RETURNED IN THE
C             DOUBLE PRECISION VARIABLE NUMBER).  THIS ROUTINE IS
C             SIMILAR TO SRINT.
C
      SUBROUTINE SRDP(N,CHAR,BEGIN,END,NUMBER,ERR)
      INTEGER             DIGIT
      INTEGER             N
      INTEGER             BLANK,       D,           E,           MINUS
      INTEGER             PERIOD,      PLUS
      DOUBLE PRECISION    NUMBER,      Q,           SIGN
      INTEGER             BEGIN,       CHAR(N),     COUNT,       END
      INTEGER             EXPONT,      FIRST,       I,           LEFT
      LOGICAL             ERR,         LCOUNT
      DATA  BLANK,PERIOD,PLUS,MINUS,D,E / 1H , 1H., 1H+, 1H-, 1HD, 1HE /
      NUMBER = 0.0D0
      SIGN   = 1.0D0
      EXPONT = 0
      LCOUNT = .FALSE.
      COUNT  = 0
      ERR    = .TRUE.
C
C  FIND THE FIRST NON-BLANK CHARACTER
C
           DO 100 I = BEGIN,END
           IF ( CHAR(I).EQ.BLANK ) GO TO 100
           FIRST = I
           GO TO 200
  100      CONTINUE
      GO TO 900
C
C  CHECK FOR SIGN
C
  200 IF ( CHAR(FIRST).NE.PLUS.AND.CHAR(FIRST).NE.MINUS ) GO TO 400
      IF ( CHAR(FIRST).NE.MINUS ) GO TO 300
      SIGN = -1.0D0
  300 FIRST = FIRST+1
      IF ( FIRST.GT.END ) GO TO 900
C
C  START MAIN CHARACTER RECOGNITION LOOP
C
  400      DO 700 I = FIRST,END
           IF ( CHAR(I).EQ.BLANK  ) GO TO 700
           IF ( CHAR(I).NE.D.AND.CHAR(I).NE.E ) GO TO 500
           LEFT = I+1
           GO TO 800
  500      IF ( LCOUNT ) COUNT = COUNT+1
           IF ( CHAR(I).NE.PERIOD ) GO TO 600
           LCOUNT = .TRUE.
           GO TO 700
  600      Q = DIGIT(CHAR(I))
           IF ( Q.LT.0.0D0 ) GO TO 900
           NUMBER = NUMBER*10.0D0+Q
  700      CONTINUE
      NUMBER = SIGN*NUMBER/10.0D0**COUNT
      RETURN
C
C  PROCESS THE EXPONENT
C
  800 CALL SRINT(N,CHAR,LEFT,END,EXPONT,ERR)
      NUMBER = SIGN*NUMBER*10.0D0**(EXPONT-COUNT)
      RETURN
C
C  SET ERROR FLAG
C
  900 ERR = .FALSE.
      RETURN
      END
      SUBROUTINE SRINT(N,CHAR,BEGIN,END,NUMBER,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SRINT
C  PURPOSE:  SUBROUTINE TO RECOGNIZE CHARACTER STRINGS REPRESENTING
C             INTEGERS.  IF THERE ARE ANY ILLEGAL CHARACTERS (NOT
C             INCLUDING BLANKS) ERR = .FALSE. .  OTHERWISE, IF THE
C             STRING DOES REPRESENT AN INTEGER, ERR = .TRUE., AND
C             NUMBER IS THE VALUE OF THE INTEGER.
C
      INTEGER             DIGIT
      INTEGER             N
      INTEGER             BLANK,       MINUS,       PLUS
      INTEGER             BEGIN,       CHAR(N),     END,         I
      INTEGER             K,           NUMBER,      SIGN
      LOGICAL             ERR
      DATA  BLANK,PLUS,MINUS / 1H , 1H+, 1H- /
      NUMBER = 0
      ERR    = .TRUE.
      SIGN   = 1
           DO 200 I = BEGIN,END
           IF ( CHAR(I).EQ.BLANK.OR.CHAR(I).EQ.PLUS ) GO TO 200
           IF ( CHAR(I).NE.MINUS ) GO TO 100
           SIGN = -1
           GO TO 200
  100      K = DIGIT(CHAR(I))
           IF ( K.EQ.(-1) ) GO TO 300
           NUMBER = NUMBER*10+K
  200      CONTINUE
      NUMBER = SIGN*NUMBER
      RETURN
  300 ERR    = .FALSE.
      RETURN
      END
      SUBROUTINE SRREAD(DEVICE,DARG,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SRREAD
C  PURPOSE:  TO READ IN A SINGLE DOUBLE PRECISION NUMBER,FREE FORMATTED,
C             FROM THE SPECIFIED DEVICE.
C
      INTEGER             LENGTH
      DOUBLE PRECISION    DARG
      INTEGER             BEGIN,       CHAR(72),    DEVICE,      ENDE
      INTEGER             I,           I72
      LOGICAL             ERR
      DATA  I72 / 72 /
C
C  SET ERR = .TRUE. (I.E. NO ERRORS)
C
      ERR  = .TRUE.
C
C  READ CHARACTER STRING FROM INPUT DEVICE AND RECOGNIZE NUMBER
C
      DARG = 0.0D0
      READ (DEVICE,10000) (CHAR(I),I=1,72)
           DO 200 I = 1,72
           CALL LCUC(CHAR(I))
  200      CONTINUE
      ENDE   = LENGTH(I72,CHAR)
      IF ( ENDE.EQ.0 ) GO TO 300
      BEGIN= 1
      CALL SRDP(I72,CHAR,BEGIN,ENDE,DARG,ERR)
  300 RETURN
C
C  FORMAT STATEMENT
C
10000 FORMAT(72A1)
      END
      SUBROUTINE SSHIFT(ISHIFT)
C CARRY OUR A SHIFT BY ISHIFT DISPLAY UNITS
      INTEGER             MAX0,        MIN0
      INTEGER             I,           IBEG,        ICENTR,      IEND
      INTEGER             IFROM,       ILEFTP,      IRGHTP,      ISHIFT
      INTEGER             ITO,         IXSIZE,      LEFT,        RIGHT
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      DATA ICENTR,IXSIZE / 2689, 5377 /
C
C
      CALL SFBNDS(ILEFT(2),IRIGHT(2))
      LEFT   = ILEFT(2)-ISHIFT
      RIGHT  = IRIGHT(2)-ISHIFT
      IF ( .NOT.( (LEFT.LT.1.AND.RIGHT.LT.1).OR.
     X     (LEFT.GT.IXSIZE.AND.RIGHT.GT.IXSIZE))) GO TO 100
      CALL ZERO
      RETURN
  100 ILEFTP = MAX0(LEFT,1)
      IRGHTP = MIN0(RIGHT,IXSIZE)
      IF ( ISHIFT.LE.0 ) GO TO 600
      IEND = IRGHTP-ILEFTP+1
           DO 300 I = 1,IEND
           IFROM     = IRGHTP+1-I
           ITO       = IFROM+ISHIFT
           LDEF(ITO) = LDEF(IFROM)
           FS(ITO)   = FS(IFROM)
  300      CONTINUE
      IEND = ILEFTP+ISHIFT
           DO 400 I = 1,IEND
           LDEF(I)   = .FALSE.
           FS(I)     = 0.0D0
  400      CONTINUE
      IBEG = IRGHTP+ISHIFT+1
      IF ( IBEG.GT.IXSIZE ) GO TO 1100
           DO 500 I = IBEG,IXSIZE
           LDEF(I)   = .FALSE.
           FS(I)     = 0.0D0
  500      CONTINUE
      GO TO 1100
  600      DO 700 I = ILEFTP,IRGHTP
           ITO       = I+ISHIFT
           LDEF(ITO) = LDEF(I)
           FS(ITO)   = FS(I)
  700      CONTINUE
      IEND = ILEFTP+ISHIFT-1
      IF ( IEND.LT.1 ) GO TO 900
           DO 800 I = 1,IEND
           LDEF(I)   = .FALSE.
           FS(I)     = 0.0D0
  800      CONTINUE
  900 IBEG = IRGHTP+ISHIFT+1
      IF ( IBEG.GT.IXSIZE ) GO TO 1100
           DO 1000 I = IBEG,IXSIZE
           LDEF(I)   = .FALSE.
           FS(I)     = 0.0D0
 1000      CONTINUE
C
C  RESET THE LEFT AND RIGHT BOUNDARIES
C
 1100 ILEFT(2)  = ILEFTP+ISHIFT
      IRIGHT(2) = IRGHTP+ISHIFT
      RETURN
      END
      SUBROUTINE SSMPUD(LMAG,LSHIFT,NUM,DEN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SSMPUD
C  PURPOSE:  TO COPY THE FS() AND LDEF() ARRAYS INTO THEMSELVES
C             CORRESPONDING TO SHIFTS AND STEP-SIZE CHANGES
C
      INTEGER             COPY,        PWROF2
      DOUBLE PRECISION    D,           DM,          DNH,         DSHIFT
      INTEGER             DEN,         I,           IBEGIN,      IC2
      INTEGER             ICENTR,      IEND,        IFROM,       ILAST
      INTEGER             INDX,        INDXP,       IP,          ISHIFT
      INTEGER             ITEMP,       ITO,         IXSIZE,      J
      INTEGER             J1,          J2,          LEFT,        MM
      INTEGER             NUM,         RIGHT,       W2
      INTEGER             M,           SHIFT
      LOGICAL             LMAG,        LSHIFT
      DOUBLE PRECISION    H
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      EQUIVALENCE (SHIFT,IOPSH(2)),(M,IOPSS(2)),(H,ROPSTS(2))
      DATA ICENTR,IXSIZE / 2689, 5377 /
C
C  FIRST SHIFT THE DATA IF LSHIFT = .TRUE.
C
      IF ( .NOT.LSHIFT ) GO TO 300
      ISHIFT  = -SHIFT*NH(1)
      CALL SSHIFT(ISHIFT)
      DSHIFT  = ISHIFT
      DSHIFT  = DSHIFT*ROPSTS(1)/8.0D0
      IEND    = IOPDM(2)
           DO 200 I = 1,IEND
           ROPPNT(I,2) = ROPPNT(I,2)-DSHIFT*ROPUDI(I,2)
  200      CONTINUE
      SHIFT   = 0
C
C  NEXT, IF LMAG = .TRUE. COPY THE FS() AND LDEF() ARRAYS INTO
C   THEMSELVES, COMPRESSING OR EXPANDING, DEPENDING ON THE VALUE
C   OF GAMMA = NUM/DEN.
C
  300 IF ( .NOT.LMAG ) GO TO 2100
      IF ( SHIFT.EQ.0 ) GO TO 500
      ISHIFT = -SHIFT*NH(1)
      CALL SSHIFT(ISHIFT)
      DSHIFT = ISHIFT
      DSHIFT = DSHIFT*ROPSTS(1)/8.0D0
      IEND   = IOPDM(2)
           DO 400 I = 1,IEND
           ROPPNT(I,2) = ROPPNT(I,2)-DSHIFT*ROPUDI(I,2)
  400      CONTINUE
      SHIFT  = 0
C
C  DETERMINE WHETHER OR NOT GAMMA IS A POWER OF 2.  IF SO, THE PROGRAM
C    IS SET UP TO SAVE PREVIOUS VALUES OF FS() AND NOT RECALCULATE THEM.
C    IF GAMMA IS NOT A POWER OF 2, PWROF2()=0, AND ALL PREVIOUS DATA IS
C    LOST.  FOR GAMMA A POWER OF 2, THE PROGRAM TREATS GAMMA>1 FIRST AND
C    THEN GAMMA<1 SECOND.
C
  500 IF ( NUM .EQ. DEN ) GO TO 2100
      IF ( .NOT.(PWROF2(NUM).NE.0 .AND. DEN.EQ.1) .AND..NOT.(NUM.EQ.1
     X     .AND. PWROF2(DEN).NE.0) ) GO TO 1800
      IF ( NUM.LT.DEN ) GO TO 1400
      MM     = PWROF2(NUM)
      LEFT   = ILEFT(2)
      RIGHT  = IRIGHT(2)
      IF ( LEFT.EQ.1 ) GO TO 700
      INDX = COPY(MM,LEFT)-1
      IF ( INDX.LT.1 ) GO TO 700
           DO 600 I = 1,INDX
           LDEF(I) = .FALSE.
  600      CONTINUE
  700 IF ( RIGHT.EQ.IXSIZE ) GO TO 900
      INDX = COPY(MM,RIGHT)+1
      IF (INDX.GT.IXSIZE ) GO TO 900
           DO 800 I = INDX,IXSIZE
           LDEF(I) = .FALSE.
  800      CONTINUE
  900 IC2    = 2*ICENTR
      IP     = IXSIZE-ICENTR-1
      INDXP  = COPY(MM,IXSIZE)
           DO 1300 I = 1,IP
           IFROM  = IXSIZE-I
           INDX   = COPY(MM,IFROM)
           IF ( INDX.GT.IXSIZE ) GO TO 1200
           IF ( INDXP.GT.IXSIZE ) INDXP = IXSIZE+1
           LDEF(INDX) = LDEF(IFROM)
           FS(INDX)   = FS(IFROM)
           LDEF(IFROM)= .FALSE.
           IEND       = INDXP-1
           IBEGIN     = INDX+1
                DO 1000 J = IBEGIN,IEND
                LDEF(J) = .FALSE.
 1000           CONTINUE
           J1         = IC2-INDX
           J2         = IC2-IFROM
           LDEF(J1)   = LDEF(J2)
           FS(J1)     = FS(J2)
           LDEF(J2)     = .FALSE.
           ITEMP      = IBEGIN
           IBEGIN     = IC2-IEND
           IEND       = IC2-ITEMP
                DO 1100 J = IBEGIN,IEND
                LDEF(J) = .FALSE.
 1100           CONTINUE
 1200      INDXP     = INDX
 1300      CONTINUE
      GO TO 2100
C
C  NEXT, HANDLE THE CASE IN WHICH GAMMA<1.  IN THIS CASE, THE MAPPING
C   IS A CONTRACTION.
C
 1400 MM    = PWROF2(DEN)
      IC2   = 2*ICENTR
      IBEGIN= ICENTR+1
           DO 1500 I = IBEGIN,IXSIZE
           ILAST = I
           INDX  = COPY(MM,ILAST)
           IF ( INDX.GT.IXSIZE ) GO TO 1600
           LDEF(I)   = LDEF(INDX)
           FS(I)     = FS(INDX)
           INDX      = IC2-INDX
           ITO       = IC2-I
           LDEF(ITO) = LDEF(INDX)
           FS(ITO)   = FS(INDX)
 1500      CONTINUE
 1600      DO 1700 I = ILAST,IXSIZE
           J       = IC2-I
           LDEF(I) = .FALSE.
           LDEF(J) = .FALSE.
 1700      CONTINUE
      GO TO 2100
C
C  GAMMA IS NOT A POWER OF 2.  THEREFORE, CALL SAMPLE ROUTINE.
C   NOTE THAT NONE OF THE PREVIOUS DATA IS SAVED.
C
 1800 IF ( SHIFT.EQ.0 ) GO TO 2000
      DSHIFT = ISHIFT
      DSHIFT = DSHIFT*ROPSTS(1)/8.0D0
      IEND   = IOPDM(2)
           DO 1900 I = 1,IEND
           ROPPNT(I,2) = ROPPNT(I,2)-DSHIFT*ROPUDI(I,2)
 1900      CONTINUE
      SHIFT  = 0
 2000 CALL ZERO
C
 2100 RETURN
      END
      SUBROUTINE SUBAC(ITYPE,N,OPDF1,OPDF2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBAC
C  PURPOSE:  TO ACCENTUATE THE N'TH DERIVATIVE.
C
      INTEGER             ITYPE,       N,           NPLUS1
      LOGICAL             OPDF1(7,2),  OPDF2(7,2)
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 2
      IF ( N.GE.0.AND.N.LE.6 ) GO TO 100
      ERRCOD = 1
      ERR    = .TRUE.
      GO TO 300
  100 NPLUS1= N+1
      IF ( OPDF2(NPLUS1,2) ) GO TO 200
      OPDF2(NPLUS1,2) = .TRUE.
      OPDF1(NPLUS1,2) = .TRUE.
      GO TO 300
  200 OPDF2(NPLUS1,2) = .FALSE.
  300 RETURN
      END
      SUBROUTINE SUBCC(ITYPE,IWHICH,ICH,INPUT,OUTD,GRPHC,RCRD,RSTRT)
C
C PURPOSE: CHANGE APPROPRIATE I/O CHANNEL NUMBER
C
C INPUT PARAMETERS:
C
C    IWHICH: DESCRIBES WHICH CHANNEL IS TO BE CHANGED
C    IWHICH = 1 INPUT CHANNEL NUMBER
C             2 OUTPUT CHANNEL NUMBER
C             3 GRAPHIC CHANNEL NUMBER
C             4 RECORDING CHANNEL NUMBER
C             5 RESTART DEVICE CHANNEL NUMBER
C
C    ICH:     THE NEW CHANNEL NUMBER
C
C ON OUTPUT, ICH WILL HAVE BEEN ASSIGNED TO THE APPROPRIATE CHANNEL
      INTEGER             ICH,         RCRD,        OUTD,        GRPHC
      INTEGER             IWHICH,      ITYPE,       RSTRT,       INPUT
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      ITYPE = 4
      ERR = .FALSE.
      IF (IWHICH.GT.0.AND.IWHICH.LT.6) GO TO 100
      ERR = .TRUE.
      ERRCOD = 16
      GO TO 200
  100 CONTINUE
      IF (IWHICH.EQ.1) INPUT = ICH
      IF (IWHICH.EQ.2) OUTPUT = ICH
      IF (IWHICH.EQ.2) OUTD = ICH
      IF (IWHICH.EQ.3) GRPHC = ICH
      IF (IWHICH.EQ.4) RCRD = ICH
      IF (IWHICH.EQ.5) RSTRT = ICH
  200 RETURN
      END
      SUBROUTINE SUBCCH(ITYPE,I,DARG)
C PURPOSE: REPLACE THE I-TH COMPONENT OF THE CURRENT CROSS DIRECTION
      DOUBLE PRECISION    DARG,        DNW
      INTEGER             I,           ITYPE
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      DOUBLE PRECISION    SUM
      EQUIVALENCE (NDIM,IOPDM(2))
      IF (I.LE.NDIM) GO TO 100
      ERRCOD = 21
      ERR = .TRUE.
      GO TO 200
  100 CONTINUE
C CHECK IF ASSIGNMENT WOULD LEAD TO ZERO DIRECTION
      SUM = 0.0D0
      DO 150 J = 1,NDIM
      IF (I.EQ.J) SUM = SUM + DARG**2
      IF (I.NE.J) SUM = SUM + CBDU(J)**2
150   CONTINUE
      IF (SUM.NE.0.0D0) GO TO 180
      ERR = .TRUE.
      ERRCOD = 18
      GO TO 200
180   CONTINUE
      LCBD = .TRUE.
      IF (ICBD.NE.0) GO TO 190
      ICBD = 1
190   CONTINUE
      DNW = NW(2)
      IF (CBDSW.EQ.0.0D0) CBDSW = (DNW*ROPSTS(2))/8.0D0
      ITYPE = 4
      IF (CBDU(I).EQ.DARG) GO TO 200
      CBDU(I) = DARG
      ITYPE = 1
      CALL NRMLC
      LCBD = .TRUE.
      LSCRN = .FALSE.
  200 RETURN
      END
      SUBROUTINE SUBCD(ITYPE,DIR)
C LOAD THE CROSS DIRECTION OF DIFFERENTIATION
C
C
      DOUBLE PRECISION    DIR(3),      DNW
      INTEGER             I,           NDIM,        ITYPE
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
C
      NDIM = IOPDM(2)
           DO 100 I = 1,NDIM
           IF (DIR(I).NE.0.0D0) GO TO 200
  100      CONTINUE
      ERR = .TRUE.
      ERRCOD = 18
      GO TO 700
  200 CONTINUE
      CALL MAKSHF
      LCBD = .TRUE.
      IF (ICBD.NE.0) GO TO 300
      ICBD = 1
      GO TO 500
  300 CONTINUE
      NDIM = IOPDM(2)
           DO 400 I = 1,NDIM
           IF (CBDU(I).NE.DIR(I)) GO TO 500
  400      CONTINUE
      GO TO 700
  500 CONTINUE
      DNW = NW(2)
      IF (CBDSW.EQ.0.0D0) CBDSW = (DNW*ROPSTS(2))/8.0D0
           DO 600 I = 1,NDIM
           CBDU(I) = DIR(I)
  600      CONTINUE
      CALL NRMLC
      ITYPE = 1
  700 CONTINUE
      RETURN
      END
      SUBROUTINE SUBCH(ITYPE,H)
C INCORPORATE THE DISCRETIZATION PARAMETER FOR THE CROSS DERIVATIVE
      DOUBLE PRECISION    H
      INTEGER             IORD,        ITYPE
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      IF (H.NE.0.0D0) GO TO 100
      ERR = .TRUE.
      ERRCOD = 19
      GO TO 200
  100 CONTINUE
      ITYPE = 4
      IF (H.EQ.CBDSW) GO TO 200
      CALL MAKSHF
      LCBD = .TRUE.
      LSCRN = .FALSE.
      ITYPE = 1
      CBDSW = H
  200 CONTINUE
      IF (ICBD.GT.0) GO TO 300
      ICBD = 1
  300 CONTINUE
      RETURN
      END
      SUBROUTINE SUBCO(ITYPE,IORD)
C INCORPORATE THE ORDER OF THE CROSS DERIVATIVE
      DOUBLE PRECISION    DNW
      INTEGER             IORD,        ITYPE
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      IF (IORD.GE.0.AND.IORD.LE.6) GO TO 100
      ERR = .TRUE.
      ERRCOD = 17
      GO TO 200
  100 CONTINUE
      CALL MAKSHF
      DNW = NW(2)
      IF (CBDSW.EQ.0.0D0) CBDSW = (DNW*ROPSTS(2))/8.0D0
      ITYPE = 4
      IF (IORD.EQ.ICBD) GO TO 200
      LCBD = .TRUE.
      LSCRN = .FALSE.
      ITYPE = 1
      ICBD = IORD
  200 CONTINUE
      RETURN
      END
      SUBROUTINE SUBCW(ITYPE,N,IOPWN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBCW
C  PURPOSE:  TO CHANGE THE STENCIL WIDTH FOR CALCULATION OF THE
C             DERIVATIVES.
C
      INTEGER             ICENTR,      IOPWN(2),    ITYPE,       IXSIZE
      INTEGER             N
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DATA ICENTR,IXSIZE / 2689, 5377 /
      ITYPE    = 2
      IF ( N.NE.0 ) GO TO 100
      ERRCOD = 2
      ERR    = .TRUE.
      RETURN
  100 IOPWN(2) = N
      RETURN
      END
      SUBROUTINE SUBDC(ITYPE,OPDC)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDC
C  PURPOSE:  TO TOGGLE THE LOGICAL FLAG OPDC WHICH CONTROLS THE DISPLAY
C             OF A 'CENTER MARK'
C
      INTEGER             ITYPE
      LOGICAL             OPDC(2)
      ITYPE = 3
      IF ( OPDC(2) ) GO TO 100
      OPDC(2) = .TRUE.
      GO TO 200
  100 OPDC(2) = .FALSE.
  200 RETURN
      END
      SUBROUTINE SUBDCH(ITYPE,I,DARG)
C PURPOSE: REPLACE THE I-TH COMPONENT OF THE CURRENT POINT
      INTEGER ITYPE,I,NDIM
      DOUBLE PRECISION DARG
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE (NDIM,IOPDM(2))
      IF (I.LE.NDIM) GO TO 100
      ERRCOD = 20
      ERR = .TRUE.
      GO TO 300
  100 CONTINUE
      ITYPE = 4
      IF (ROPDI(I,2).EQ.DARG) GO TO 300
      ROPDI(I,2) = DARG
      ITYPE = 1
           DO 200 J = 1,NDIM
           IF (ROPDI(J,2).NE.0.0D0) GO TO 300
  200      CONTINUE
      ERRCOD = 12
      ERR = .TRUE.
  300 RETURN
      END
      SUBROUTINE SUBDG(ITYPE,N,OPDF1)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDG
C  PURPOSE:  TO SET THE LOGICAL FLAG OPDF1 WHICH CAUSES THE N'TH
C             DERIVATIVE TO BE DISPLAYED.
C
      INTEGER             ITYPE,       N
      LOGICAL             OPDF1(7,2)
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 2
      IF ( N.GE.0.AND.N.LE.6 ) GO TO 100
      ERRCOD = 4
      ERR    = .TRUE.
      GO TO 200
  100 OPDF1(N+1,2) = .TRUE.
  200 RETURN
      END
      SUBROUTINE SUBDI(ITYPE,N,IOPSS,IOPWN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDI
C  PURPOSE:  TO DIVIDE THE CURRENT STEP-SIZE INDEX NH(2) BY THE INTEGER
C             ARGUMENT N.  N MUST BE POSITIVE.  IF IT SO LARGE AS TO
C             FORCE THE PROGRAM TO DEFINE A NEW XS() ARRAY, THE PROGRAM
C             WILL RECALCULATE THE XS() AND FS() ARRAYS ACCORDINGLY.
C             NOTE THAT THE PROGRAM IS SET UP TO TREAT THE CASE N IS A
C             POWER OF 2 EFFICIENTLY SO THAT DATA CAN BE PRESERVED.
C             THE ROUTINE ALSO CHANGES THE WINDOW IN A SIMILAR MANNER.
C
      INTEGER             IOPSS(2),    IOPWN(2),    ITYPE,       N
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 2
      IF ( N.NE.0 ) GO TO 100
      ERRCOD = 5
      ERR    = .TRUE.
      RETURN
  100 IOPSS(2) = -N
      IOPWN(2) = -N
      RETURN
      END
      SUBROUTINE SUBDM(ITYPE,N,IOPDM,ROPDI,ROPUDI,ROPPNT)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDM
C  PURPOSE:  TO SET THE PROGRAM TO WORK WITH 1,2 OR 3 DIMENSIONAL DOMAIN
C
      DOUBLE PRECISION    DSQRT
      DOUBLE PRECISION    DMNSN,       ROPDI(3,2),  ROPPNT(3,2)
      DOUBLE PRECISION    ROPUDI(3,2)
      INTEGER             I,           IOPDM(2),    ITYPE,       N
      INTEGER             NP
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      ITYPE    = 1
      IF ( N.GT.0 .AND. N.LT.4 ) GO TO 100
      ERRCOD = 6
      ERR    = .TRUE.
      GO TO 400
  100 DMNSN    = N
      IOPDM(2) = N
           DO 200 I = 1,N
           ROPDI(I,2)  = 1.0D0
           ROPPNT(I,2) = 0.0D0
           CBDU(I) = 1.0D0
  200      CONTINUE
      ICBD = 0
      CBDSW = 0.
      LCBD = .FALSE.
      NP       = N+1
      IF ( NP.GT.3 ) GO TO 400
           DO 300 I = NP,3
           ROPDI(I,2)  = 0.0D0
           ROPUDI(I,2) = 0.0D0
           ROPPNT(I,2) = 0.0D0
           CBDU(I) = 0.0D0
  300      CONTINUE
  400 CONTINUE
      CALL NRML
      CALL NRMLC
      RETURN
      END
      SUBROUTINE SUBDO(ITYPE,IOPSS,IOPWN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDO
C  PURPOSE:  DOUBLE THE STEP-SIZE.  THEN THE MAIN PROGRAM WILL SET NH(2)
C             TO 2*NH(1) UNLESS THIS EXCEEDS THE XS() AND FS() ARRAY
C             BOUNDARIES, IN WHICH CASE, THE PROGRAM SETS UP A NEW
C             SAMPLING.  THE ROUTINE ALSO CHANGES THE WINDOW IN THE SAME
C             WAY.
C
      INTEGER             IOPSS(2),    IOPWN(2),    ITYPE
      ITYPE    = 2
      IOPSS(2) = 2
      IOPWN(2) = 2
      RETURN
      END
      SUBROUTINE SUBDS(ITYPE,OPDS)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDS
C  PURPOSE:  SET LOGICAL FLAG OPDS WHICH CONTROLS THE DISPLAYING OF THE
C             X-AXIS WITH A SCALE.
C
      INTEGER             ITYPE
      LOGICAL             OPDS(2)
      ITYPE = 3
      IF ( OPDS(2) ) GO TO 100
      OPDS(2) = .TRUE.
      GO TO 200
  100 OPDS(2) = .FALSE.
  200 RETURN
      END
      SUBROUTINE SUBDX(ITYPE,OPDX)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBDX
C  PURPOSE:  SET THE LOGICAL FLAG OPDX WHICH DISPLAYS THE X-AXIS.
C
      INTEGER             ITYPE
      LOGICAL             OPDX(2)
      ITYPE = 3
      IF ( OPDX(2) ) GO TO 100
      OPDX(2) = .TRUE.
      GO TO 200
  100 OPDX(2) = .FALSE.
  200 RETURN
      END
      SUBROUTINE SUBEG(ITYPE,N,OPDF1,OPDF2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBEG          (ERASE GRAPH)
C  PURPOSE:  TO SET THE LOGICAL FLAGS OPDF1(N+1,2) AND OPDF2(N+1,2) TO
C             .FALSE. SO THAT THE N'TH DERIVATIVES (N RANGING FROM 0 TO
C             6) WILL NOT BE DISPLAYED.
C
      INTEGER             ITYPE,       N
      LOGICAL             OPDF1(7,2),  OPDF2(7,2)
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 3
      IF ( N.GE.0.AND.N.LE.6 ) GO TO 100
      ERRCOD = 7
      ERR    = .TRUE.
      GO TO 200
  100 OPDF1(N+1,2) = .FALSE.
      OPDF2(N+1,2) = .FALSE.
  200 RETURN
      END
      SUBROUTINE SUBFL(ITYPE)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBFL
C  PURPOSE:  TO FLIP THE FUNCTION. THAT IS: THE DIRECTION VECTORS ARE
C             FLIPPED 180 DEGREES AND THE FS(),DF(S) ARRAYS ARE FLIPPED.
C             FURTHERMORE, THE ODD DERIVATIVES CHANGE SIGN.
C
      DOUBLE PRECISION    D,           TEMP
      INTEGER             I,           ICENTR,      IEND,        INDX
      INTEGER             IP,          IPP,         IPPP,        ITEMP
      INTEGER             ITYPE,       IXSIZE,      J,           JP
      INTEGER             JPP,         JPPP,        K,           W2
      INTEGER             W2P1
      LOGICAL             LTEMP
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DATA ICENTR,IXSIZE / 2689, 5377 /
      ITYPE = 2
C
C  FLIP THE DIRECTION VECTORS
C
      IEND = IOPDM(2)
           DO 100 I = 1,IEND
           ROPDI(I,2)  = -ROPDI(I,2)
           ROPUDI(I,2) = -ROPUDI(I,2)
           TEMP        = ROPDR1(I,2)
           ROPDR1(I,2) = ROPDR2(I,2)
           ROPDR2(I,2) = TEMP
  100      CONTINUE
C
C  FLIP THE DF() ARRAY
C
      W2   = WIDTH/2
      W2P1 = W2+1
      JP   = WIDTH+1
           DO 300 I = 1,7
           D           = (-1)**(I+1)
           DF(W2P1,I)  = D*DF(W2P1,I)
                DO 200 J = 1,W2
                K        = JP-J
                TEMP     = DF(J,I)
                DF(J,I)  = D*DF(K,I)
                DF(K,I)  = D*TEMP
  200           CONTINUE
           LPLT(I)     = .FALSE.
  300      CONTINUE
C
C  FLIP THE FS() AND LDEF() ARRAY
C
      IP = ICENTR-1
      JP = IXSIZE+1
           DO 400 I = 1,IP
           INDX       = JP-I
           TEMP       = FS(I)
           LTEMP      = LDEF(I)
           FS(I)      = FS(INDX)
           LDEF(I)    = LDEF(INDX)
           FS(INDX)   = TEMP
           LDEF(INDX) = LTEMP
  400      CONTINUE
      RETURN
      END
      SUBROUTINE SUBFO(ITYPE,LGO)
C PURPOSE: MAKE MICROSCOPE RECALCULATE EVERYTHING
      INTEGER             ITYPE
      LOGICAL             LGO
      LOGICAL             LFO
      COMMON / FOOWN  /   LFO
      LGO = .TRUE.
      ITYPE = 2
      LFO = .TRUE.
      RETURN
      END
      SUBROUTINE SUBGO(ITYPE,LGO)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBGO
C  PURPOSE:  TO SET PROGRAM SO THAT IT BEGINS COMPUTATION
C
      INTEGER             ITYPE
      LOGICAL             LGO
      ITYPE = 2
      LGO   = .TRUE.
      RETURN
      END
      SUBROUTINE SUBHA(ITYPE,IOPSS,IOPWN)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBHA
C  PURPOSE:  TO SET THE PROGRAM TO HALVE THE STEP-SIZE.
C
      INTEGER             IOPSS(2),    IOPWN(2),    ITYPE
      ITYPE    = 2
      IOPSS(2) = -2
      IOPWN(2) = -2
      RETURN
      END
      SUBROUTINE SUBID(ITYPE,V,ROPDI,ROPUDI,IOPDM)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBID
C  PURPOSE:  TO SET THE DIRECTION OF THE LINE OF INVESTIGATION
C             (I.E. ROPDI()) AND FORM THE CORRESPONDING UNIT VECTOR
C             ROPDUDI().
C
      DOUBLE PRECISION    DSQRT
      DOUBLE PRECISION    ROPDI(3,2),  ROPUDI(3,2), SUM,         V(3)
      INTEGER             I,           IEND,        IOPDM(2),    ITYPE
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 1
      IF ( V(1).NE.0.0D0 .OR. V(2).NE.0.0D0 .OR. V(3).NE.0.0D0) GO TO
     X     100
      ERRCOD = 12
      ERR    = .TRUE.
      RETURN
  100 SUM   = 0.0D0
C MAKE SHIFTED POINT NEW POINT IF NECESSARY
      CALL MAKSHF
      IEND  = IOPDM(2)
           DO 200 I = 1,IEND
           ROPDI(I,2)  = V(I)
  200      CONTINUE
      CALL NRML
      RETURN
      END
      SUBROUTINE SUBIH(ITYPE,H,ROPSTS,NH)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBIH
C  PURPOSE:  TO SET THE DISCRETIZATION PARAMETER (I.E. H = ROPSTS(2) )
C
      DOUBLE PRECISION    DNH,         H,           ROPSTS(2)
      INTEGER             ITYPE,       NH(2)
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 1
      IF ( H.GT.0.0D0 ) GO TO 100
      ERRCOD = 11
      ERR    = .TRUE.
      RETURN
  100 DNH       = NH(2)
      ROPSTS(2) = H*8.0D0/DNH
      RETURN
      END
      SUBROUTINE SUBII(ITYPE,V1,V2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBII
C  PURPOSE:  TO SET THE ENDPOINTS OF THE LINE OF INVESTIGATION AND
C             FROM THEM CALCULATE AND SET THE DIRECTION AND CENTER POINT
C             OF THE LINE OF INVESTIGATION.
C
      DOUBLE PRECISION    DSQRT
      DOUBLE PRECISION    DNH,         SUM,         V1(3),       V2(3)
      DOUBLE PRECISION    W
      INTEGER             DMNSN,       I,           ITYPE
      LOGICAL             OPDD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      ITYPE = 1
      LSCRN = .FALSE.
      DMNSN = IOPDM(2)
      SUM = 0.0D0
           DO 100 I = 1,DMNSN
           ROPDI(I,2)  = V2(I)-V1(I)
           SUM         = SUM+ROPDI(I,2)**2
           ROPPNT(I,2) = V1(I)+ROPDI(I,2)/2.0D0
           ROPDR1(I,2) = V1(I)
           ROPDR2(I,2) = V2(I)
  100      CONTINUE
      SUM    = DSQRT(SUM)
      IF ( SUM.NE.0.0D0 ) GO TO 200
      ERRCOD = 14
      ERR    = .TRUE.
      RETURN
  200 DNH       = NH(2)
      W         = WIDTH
      ROPSTS(2) = (SUM/(W-1.0D0))*8.0D0/DNH
      CALL NRML
      RETURN
      END
      SUBROUTINE SUBIP(ITYPE,V,ROPPNT,IOPDM)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBIP
C  PURPOSE:  TO SET THE COORDINATES OF THE CENTER POINT ALONG THE
C             LINE OF INVESTIGATION.  (I.E.ROPPNT())
C
      DOUBLE PRECISION    ROPPNT(3,2), V(3)
      INTEGER             I,           IEND,        IOPDM(2),    ITYPE
      ITYPE = 1
      IEND  = IOPDM(2)
           DO 100 I = 1,IEND
           ROPPNT(I,2) = V(I)
  100      CONTINUE
      RETURN
      END
      SUBROUTINE SUBLC(ITYPE,IARG)
C PURPOSE   LOAD A CROSS DIRECTION FROM DEVICE CLOAD
      DOUBLE PRECISION    V(3)
      INTEGER             I,           IARG,        IARG1,       IDUM
      INTEGER             NLINES,      ITYPE
      LOGICAL             OK
      REAL                DNW
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE (NDIM,IOPDM(2))
      CALL POS(CLOAD,IARG)
  200 CONTINUE
      CALL SVREAD(CLOAD,NLINES,NDIM,V,OK)
      IF (.NOT.OK) GO TO 600
      ITYPE = 4
           DO 300 I = 1,NDIM
           IF (V(I).NE.CBDU(I)) GO TO 400
  300      CONTINUE
      GO TO 700
  400 CONTINUE
      DO 450 I = 1,NDIM
      IF (V(I).NE.0.0D0) GO TO 460
450   CONTINUE
      GO TO 650
460   CONTINUE
      ITYPE = 1
           DO 500 I =1,NDIM
           CBDU(I) = V(I)
  500      CONTINUE
      IF (ICBD.EQ.0) ICBD = 1
      LCBD = .TRUE.
      LSCRN = .FALSE.
      CALL NRMLC
      GO TO 700
  600 CONTINUE
      ERR = .TRUE.
      ERRCOD = 24
      GO TO 800
  650 CONTINUE
      ERR = .TRUE.
      ERRCOD = 18
      GO TO 800
  700 CONTINUE
      DNW = NW(2)
      IF (CBDSW.EQ.0.0D0) CBDSW = (DNW*ROPSTS(2))/8.0D0
  800 CONTINUE
      RETURN
10000 FORMAT(A1)
      END
      SUBROUTINE SUBLD(ITYPE,IARG)
C PURPOSE   LOAD A DIRECTION FROM DEVICE DLOAD
      DOUBLE PRECISION    V(3)
      INTEGER             I,           IARG,        IARG1,       IDUM
      INTEGER             NLINES,      ITYPE
      LOGICAL             OK
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE (NDIM,IOPDM(2))
      CALL POS(DLOAD,IARG)
  200 CONTINUE
      CALL SVREAD(DLOAD,NLINES,NDIM,V,OK)
      IF (.NOT.OK) GO TO 600
      DO 250 I = 1,NDIM
      IF (V(I).NE.0.0D0) GO TO 260
250   CONTINUE
      ERRCOD = 12
      ERR = .TRUE.
      GO TO 700
260   CONTINUE
      ITYPE = 4
           DO 300 I = 1,NDIM
           IF (V(I).NE.ROPDI(I,2)) GO TO 400
  300      CONTINUE
      GO TO 700
  400 CONTINUE
      ITYPE = 1
           DO 500 I =1,NDIM
           ROPDI(I,2) = V(I)
  500      CONTINUE
      CALL NRML
      GO TO 700
  600 CONTINUE
      ERR = .TRUE.
      ERRCOD = 23
  700 CONTINUE
      RETURN
10000 FORMAT(A1)
      END
      SUBROUTINE SUBLO(ITYPE,IARG)
C CHANGE LOG CHANNEL NUMBER - ZERO FOR NO LOG
      INTEGER             IARG,        ITYPE
      INTEGER             LCHN
      COMMON / LOG    /   LCHN
      ITYPE = 3
      LCHN = IARG
      RETURN
      END
      SUBROUTINE SUBLP(ITYPE,IARG)
C PURPOSE   LOAD A POINT FORM DEVICE PLOAD
      DOUBLE PRECISION    V(3)
      INTEGER             I,           IARG,        IARG1,       IDUM
      INTEGER             NLINES,      ITYPE
      LOGICAL             OK
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      INTEGER             NDIM
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE (NDIM,IOPDM(2))
      CALL POS(PLOAD,IARG)
  200 CONTINUE
      CALL SVREAD(PLOAD,NLINES,NDIM,V,OK)
      IF (.NOT.OK) GO TO 600
      ITYPE = 4
           DO 300 I = 1,NDIM
           IF (V(I).NE.ROPPNT(I,2)) GO TO 400
  300      CONTINUE
      GO TO 700
  400 CONTINUE
      ITYPE = 1
           DO 500 I =1,NDIM
           ROPPNT(I,2) = V(I)
  500      CONTINUE
      GO TO 700
  600 CONTINUE
      ERR = .TRUE.
      ERRCOD = 22
  700 CONTINUE
      RETURN
10000 FORMAT(A1)
      END
      SUBROUTINE SUBMU(ITYPE,N,IOPSS,IOPWN)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBMU
C  PURPOSE:  TO SET THE PROGRAM TO MULTIPLY THE APPARENT STEP-SIZE (HP)
C             BY THE USER SUPPLIED INTEGER NUMBER N.  THIS COMMAND
C             KEEPS THE APPARENT WINDOW FOR CALCULATING THE DERIVATIVE
C             THE SAME.
C
      INTEGER             IOPSS(2),    IOPWN(2),    ITYPE,       N
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 2
      IF ( N.NE.0 ) GO TO 100
      ERRCOD = 13
      ERR    = .TRUE.
      RETURN
  100 IOPSS(2) = N
      IOPWN(2) = N
      RETURN
      END
      SUBROUTINE SUBNO(ITYPE)
C PURPOSE: TO TURN ON OR OFF THE NORMALIZTION OPTION
      INTEGER             ITYPE
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
      ITYPE = 1
      LSCRN = .FALSE.
      NORMAL = .NOT.NORMAL
      CALL NRML
      CALL NRMLC
      RETURN
      END
      SUBROUTINE SUBOU(ITYPE,RECORD,NCALLS,LGO)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBOU
C  PURPOSE:  TO WRITE AN IMAGE OF THE TERMINAL SCREEN ONTO THE OUTPUT
C             FILE.
C
      INTEGER             ITYPE,       NCALLS,      RECORD
      LOGICAL LGO
      ITYPE = 4
      CALL SCROLL(RECORD,NCALLS,LGO)
      RETURN
      END
      SUBROUTINE SUBPCH(ITYPE,I,DARG)
C PURPOSE: REPLACE THE I-TH COMPONENT OF THE CURRENT POINT
      INTEGER ITYPE,I,NDIM
      DOUBLE PRECISION DARG
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      EQUIVALENCE (NDIM,IOPDM(2))
      IF (I.LE.NDIM) GO TO 100
      ERRCOD = 3
      ERR = .TRUE.
      GO TO 200
  100 CONTINUE
      ITYPE = 4
      IF (ROPPNT(I,2).EQ.DARG) GO TO 200
      ROPPNT(I,2) = DARG
      ITYPE = 1
  200 RETURN
      END
      SUBROUTINE SUBRC(ITYPE,IARG)
C PURPOSE:  TO READ THE DEVICE NUMBER FOR LOADING CROSS DIRECTIONS
      INTEGER             IARG,        ITYPE
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      ITYPE = 3
      CLOAD = IARG
      RETURN
      END
      SUBROUTINE SUBRD(ITYPE,IARG)
C PURPOSE:  TO READ THE DEVICE NUMBER FOR LOADING DIRECTIONS
      INTEGER             IARG,        ITYPE
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      ITYPE = 3
      DLOAD = IARG
      RETURN
      END
      SUBROUTINE SUBRE(ITYPE,NCALLS,LGO)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBRE
C  PURPOSE:  TO RESTORE ALL ARRAYS IN THE FILE ASSIGNED TO DEVICE:
C             RSTRTD
C             SO THAT WORK CAN RESUME AT A LATER DATE.  THE COMMAND FOR
C             STORING IS THE RESTORE (RS) COMMAND.
C
C
C  THIS ROUTINE RESTORES THE CONTENTS OF ALL COMMON BLOCKS EXCEPT FOR:
C      ERRCOM (WHICH IS NEEDED ONLY WITHIN EACH STEP OF THE MAIN LOOP)
C      HELPER (WHICH IS ESSENTIALLY PART OF THE PROGRAM AND SHOULD
C              REMAIN UNCHANGED THROUGHOUT THE LIFTIME OF THIS VERSION
C              OF MICROSCOPE)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     NON-COMMON VARIABLES
C
C     SETS                ITYPE
C
      LOGICAL LGO
      INTEGER             NCALLS,      ITYPE
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / CB     /
C
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / IO     /
C
      INTEGER             INPUTD,      GRAPHD,      HELPD
      INTEGER             RECORD,      RSTRTD
      COMMON / IO     /   INPUTD,      GRAPHD,      HELPD
      COMMON / IO     /   RECORD,      RSTRTD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOG    /
C
      INTEGER             LCHN
      COMMON / LOG    /   LCHN
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOGCOM /
C
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / FUNCOM /
C
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / PLTCOM /
C
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOADIO /
C
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / SCREEN /
C
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / OPTION /
C
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / NRMLZE /
C
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / ROOM   /
C
      INTEGER             ILPUSR,      IDSUSR
      COMMON / ROOM   /   ILPUSR,      IDSUSR
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / PLOWN  /
C
C     SETS                BL,          BR,          NUMBR,       TIME
C     SETS                DATE,        LMARK,       LBLS,        ITITLE
C     SETS                IBOTTM
C
      INTEGER             ITITLE,      IBOTTM,      NUMBR,       BL
      INTEGER             BR
      LOGICAL             FRAME,COLOR,NUMRCL
      LOGICAL             LMARK,       LBLS,        DATE,        TIME
      COMMON / PLOWN  /
     X ITITLE(72),IBOTTM(72),NUMBR,FRAME,COLOR,NUMRCL,LMARK
      COMMON / PLOWN  /   LBLS,        DATE,        TIME,        BL(72)
      COMMON / PLOWN  /   BR(72)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / FOOWN  /
C
      LOGICAL             LFO
      COMMON / FOOWN  /   LFO
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / USER  /
C
      DOUBLE PRECISION    ETA
      INTEGER             IROUND,      N
      LOGICAL             ADD
      COMMON / USER   /   ETA,         IROUND,      N,           ADD
C
      ITYPE = 4
C
C  READ IN THE VARIBALES IN THE PLOWN COMMON BLOCK
C
      READ (RSTRTD) ITITLE,IBOTTM,NUMBR,FRAME,COLOR,NUMRCL,LMARK,LBLS,
     X     DATE,TIME,BL,BR
C
C  READ IN THE VARIABLE IN FOOWN
C
      READ (RSTRTD)  LFO
C
C  READ IN THE VARIABLES IN THE ROOM COMMON BLOCK
C
      READ (RSTRTD) ILPUSR,IDSUSR
C
C  READ IN THE NORMALIZTION CONTROL
C
      READ (RSTRTD) NORMAL
C
C  READ IN ALL THE VARIABLES IN THE SCREEN COMMON BLOCK
C
      READ (RSTRTD) OUTPUT,LINES,WIDTH,ILP,IDSPLA,IPRMPT,LSCRN
C
C  READ IN ALL THE VARIABLES IN THE LOADIO COMMON BLOCK
C
      READ (RSTRTD) PLOAD,DLOAD,CLOAD
C
C  READ IN ALL THE VARIABLES IN THE PLTCOM COMMON BLOCK
C
      READ (RSTRTD) SCALE,IPLOT,ISCRN1,ISCRN2
C
C  READ IN THE VARIABLES IN THE LOGCOM COMMON BLOCK
C
      READ (RSTRTD) LDF,LPLT,LDEF
C
C  READ IN THE CHANNEL NUMBER IN THE LOG COMMON BLOCK
C
      READ (RSTRTD) LCHN
C
C  READ IN ALL THE VARIABLES IN THE IO COMMON BLOCK
C
      READ (RSTRTD) INPUTD,GRAPHD,HELPD,RECORD,RSTRTD
C
C  READ IN ALL THE VARIABLES IN THE CB COMMON BLOCK
C
      READ (RSTRTD) CBDSW,CBD,CBDU,ICBD,LCBD
C
C  READ IN ALL THE VARIABLES IN THE OPTION COMMON BLOCK
C
      READ (RSTRTD) ROPDI,ROPPNT,ROPDR1,ROPDR2,ROPSTS,ROPUDI
      READ (RSTRTD) NH,NW,ILEFT,IRIGHT,IOPDM,IOPSH,IOPSS,IOPWN
      READ (RSTRTD) OPDC,OPDS,OPDX,OPCMP,OPSPL,OPDF1,OPDF2,ROPSTW
C
C  READ IN ALL THE VARIABLES IN THE FUNCOM COMMON BLOCK
C
      READ (RSTRTD) XS,FS,DF,DFMNMX
C
C  READ IN THE NUMBER OF FUNCTION EVALUATIONS
C
      READ (RSTRTD) NCALLS,LGO
C
C  READ IN THE USER COMMON BLOCK (USED FOR THE TEST PACKAGE)
C
      READ (RSTRTD) ETA,IROUND,N,ADD
      REWIND RSTRTD
      RETURN
      END
      SUBROUTINE SUBRO(ITYPE,V,ROPDI,ROPUDI,IOPDM)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBRO
C  PURPOSE:  TO "ROTATE" THE DIRECTION VECTOR ROPDI() BY ADDING THE
C             USER INPUT VECTOR V() TO THE CURRENT VECTOR ROPDI().
C
      DOUBLE PRECISION    V(3),        SUM,         ROPDI(3,2),  DSQRT
      DOUBLE PRECISION    ROPUDI(3,2)
      INTEGER             I,           IEND,        IOPDM(2),    ITYPE
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 1
      IEND  = IOPDM(2)
C MAKE SURE WE WILL NOT OBTAIN A ZERO DIRECTION OF INVESTIGATION
           DO 100 I = 1,IEND
           IF (V(I).NE.(-ROPDI(I,2))) GO TO 200
  100      CONTINUE
      ERRCOD  = 12
      ERR = .TRUE.
      GO TO 400
  200 CONTINUE
C MAKE SHIFT IF NECESSARY
      CALL MAKSHF
           DO 300 I = 1,IEND
           ROPDI(I,2)  = ROPDI(I,2)+V(I)
  300      CONTINUE
      IEND  = IOPDM(2)
      CALL NRML
  400 CONTINUE
      RETURN
      END
      SUBROUTINE SUBRP(ITYPE,IARG)
C PURPOSE:  TO READ THE DEVICE NUMBER FOR LOADING POINTS
      INTEGER             IARG,        ITYPE
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
      ITYPE = 3
      PLOAD = IARG
      RETURN
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBRS
C  PURPOSE:  TO SET FLAGS SO THAT THE SCREEN IMAGE IS REFRESHED, WITHOUT
C             BEING UPDATED.  THE PURPOSE OF THIS ROUTINE IS TO REFRESH
C             THE SCREEN EASILY WHEN OPERATING SYSTEM MESSAGES, ETC.
C             HAVE CAUSED IT TO BE OVERWRITTEN.
C
      SUBROUTINE SUBRS(ITYPE,LSCRN)
      INTEGER             ITYPE
      LOGICAL             LSCRN
      ITYPE = 3
      LSCRN = .FALSE.
      RETURN
      END
      SUBROUTINE SUBRW(ITYPE,IDEVCE)
C REWIND DEVICE NUMBER IDEVCE
      INTEGER             ITYPE,       IDEVCE
      ITYPE = 4
      REWIND IDEVCE
      RETURN
      END
      SUBROUTINE SUBSE(ITYPE)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBSE
C  PURPOSE:  TO RESET ALL THE PARAMETERS TO THE DEFAULTS
C
      INTEGER             ITYPE
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      LSCRN = .FALSE.
      ITYPE = 1
      CALL DFAULT
      CALL ZERO
      RETURN
      END
      SUBROUTINE SUBSH(ITYPE,N,IOPSH)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBSH
C  PURPOSE:  TO SET THE SHIFT PARAMETER (I.E. SHIFT = IOPSH(2))
C
      INTEGER             IOPSH(2),    ITYPE,       N
      ITYPE    = 2
      IOPSH(2) = IOPSH(2)-N
      RETURN
      END
      SUBROUTINE SUBST(ITYPE,NCALLS,LGO)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBST
C  PURPOSE:  TO STORE ALL RELEVANT INFORMATION SO THAT COMPUTATION
C             MAY RESUME AT A LATER TIME.  THE CORRESPONDING COMMAND
C             FOR THE RESTORATION OF THE INFORMATION IS THE RESTART
C             COMMAND
C
C  THIS ROUTINE RESTORES THE CONTENTS OF ALL COMMON BLOCKS EXCEPT FOR:
C      ERRCOM (WHICH IS NEEDED ONLY WITHIN EACH STEP OF THE MAIN LOOP)
C      HELPER (WHICH IS ESSENTIALLY PART OF THE PROGRAM AND SHOULD
C              REMAIN UNCHANGED THROUGHOUT THE LIFTIME OF THIS VERSION
C              OF MICROSCOPE)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     NON-COMMON VARIABLES
C
C     SETS                ITYPE
C
      LOGICAL LGO
      INTEGER             NCALLS,      ITYPE
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     COMMON BLOCK / CB     /
      DOUBLE PRECISION    CBDSW,       CBD,         CBDU
      INTEGER             ICBD
      LOGICAL             LCBD
      COMMON / CB     /   CBDSW,       CBD(3),      CBDU(3),     ICBD
      COMMON / CB     /   LCBD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / IO     /
C
      INTEGER             INPUTD,      GRAPHD,      HELPD
      INTEGER             RECORD,      RSTRTD
      COMMON / IO     /   INPUTD,      GRAPHD,      HELPD
      COMMON / IO     /   RECORD,      RSTRTD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOG    /
C
      INTEGER             LCHN
      COMMON / LOG    /   LCHN
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOGCOM /
C
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / FUNCOM /
C
      DOUBLE PRECISION    XS,          FS,          DF,          DFMNMX
      COMMON / FUNCOM /   XS(5377),    FS(5377),    DF(135,7)
      COMMON / FUNCOM /   DFMNMX(2,7)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / PLTCOM /
C
      INTEGER             IPLOT,       ISCRN1,      ISCRN2
      REAL                SCALE
      COMMON / PLTCOM /   SCALE(7),    IPLOT(135,7),ISCRN1(135,57)
      COMMON / PLTCOM /   ISCRN2(135,57)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / LOADIO /
C
      INTEGER             PLOAD,       DLOAD,       CLOAD
      COMMON / LOADIO /   PLOAD,       DLOAD,       CLOAD
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / SCREEN /
C
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / OPTION /
C
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / NRMLZE /
C
      LOGICAL             NORMAL
      COMMON / NRMLZE /   NORMAL
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / ROOM   /
C
      INTEGER             ILPUSR,      IDSUSR
      COMMON / ROOM   /   ILPUSR,      IDSUSR
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / PLOWN  /
C
C     SETS                BL,          BR,          NUMBR,       TIME
C     SETS                DATE,        LMARK,       LBLS,        ITITLE
C     SETS                IBOTTM
C
      INTEGER             ITITLE,      IBOTTM,      NUMBR,       BL
      INTEGER             BR
      LOGICAL             FRAME,COLOR,NUMRCL
      LOGICAL             LMARK,       LBLS,        DATE,        TIME
      COMMON / PLOWN  /
     X ITITLE(72),IBOTTM(72),NUMBR,FRAME,COLOR,NUMRCL,LMARK
      COMMON / PLOWN  /   LBLS,        DATE,        TIME,        BL(72)
      COMMON / PLOWN  /   BR(72)
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / FOOWN  /
C
      LOGICAL             LFO
      COMMON / FOOWN  /   LFO
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMMON BLOCK / USER  /
C
      DOUBLE PRECISION    ETA
      INTEGER             IROUND,      N
      LOGICAL             ADD
      COMMON / USER   /   ETA,         IROUND,      N,           ADD
C
      ITYPE = 4
C
C  WRITE OUT THE VARIBALES IN THE PLOWN COMMON BLOCK
C
      WRITE (RSTRTD) ITITLE,IBOTTM,NUMBR,FRAME,COLOR,NUMRCL,LMARK,LBLS,
     X     DATE,TIME,BL,BR
C
C  WRITE OUT THE VARIABLE IN FOOWN
C
      WRITE (RSTRTD)  LFO
C
C  WRITE OUT THE VARIABLES IN THE ROOM COMMON BLOCK
C
      WRITE (RSTRTD) ILPUSR,IDSUSR
C
C  WRITE OUT THE NORMALIZTION CONTROL
C
      WRITE (RSTRTD) NORMAL
C
C  WRITE OUT ALL THE VARIABLES IN THE SCREEN COMMON BLOCK
C
      WRITE (RSTRTD) OUTPUT,LINES,WIDTH,ILP,IDSPLA,IPRMPT,LSCRN
C
C  WRITE OUT ALL THE VARIABLES IN THE LOADIO COMMON BLOCK
C
      WRITE (RSTRTD) PLOAD,DLOAD,CLOAD
C
C  WRITE OUT ALL THE VARIABLES IN THE PLTCOM COMMON BLOCK
C
      WRITE (RSTRTD) SCALE,IPLOT,ISCRN1,ISCRN2
C
C  WRITE OUT THE VARIABLES IN THE LOGCOM COMMON BLOCK
C
      WRITE (RSTRTD) LDF,LPLT,LDEF
C
C  WRITE OUT THE CHANNEL NUMBER IN THE LOG COMMON BLOCK
C
      WRITE (RSTRTD) LCHN
C
C  WRITE OUT ALL THE VARIABLES IN THE IO COMMON BLOCK
C
      WRITE (RSTRTD) INPUTD,GRAPHD,HELPD,RECORD,RSTRTD
C
C  WRITE OUT ALL THE VARIABLES IN THE CB COMMON BLOCK
C
      WRITE (RSTRTD) CBDSW,CBD,CBDU,ICBD,LCBD
C
C  WRITE OUT ALL THE VARIABLES IN THE OPTION COMMON BLOCK
C
      WRITE (RSTRTD) ROPDI,ROPPNT,ROPDR1,ROPDR2,ROPSTS,ROPUDI
      WRITE (RSTRTD) NH,NW,ILEFT,IRIGHT,IOPDM,IOPSH,IOPSS,IOPWN
      WRITE (RSTRTD) OPDC,OPDS,OPDX,OPCMP,OPSPL,OPDF1,OPDF2,ROPSTW
C
C  WRITE OUT ALL THE VARIABLES IN THE FUNCOM COMMON BLOCK
C
      WRITE (RSTRTD) XS,FS,DF,DFMNMX
C
C  WRITE OUT THE NUMBER OF FUNCTION EVALUATIONS
C
      WRITE (RSTRTD) NCALLS,LGO
C
C  WRITE OUT THE USER COMMON BLOCK (USED FOR THE TEST PACKAGE)
C
      WRITE (RSTRTD) ETA,IROUND,N,ADD
      REWIND RSTRTD
      RETURN
      END
      SUBROUTINE SUBUN(ITYPE)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBUN
C  PURPOSE:  TO "UNDO" THE CHANGES MADE TO THE OPTION LIST.  THIS IS
C             ACCOMPLISHED BY COPYING THE "OLD OPTIONS" , (I.E. THOSE
C             WITH SUBSCRIPT 1) OVER THE NEWLY INPUTTED OPTIONS.
C
      INTEGER             I,           ITYPE
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      ITYPE   = 4
C
C  COPY BACK ALL VECTORS
C
           DO 100 I = 1,3
           ROPDI(I,2)   = ROPDI(I,1)
           ROPPNT(I,2)  = ROPPNT(I,1)
           ROPDR1(I,2)  = ROPDR1(I,1)
           ROPDR2(I,2)  = ROPDR2(I,2)
           ROPUDI(I,2)  = ROPUDI(I,1)
  100      CONTINUE
      ROPSTS(2)    = ROPSTS(1)
      NH(2)        = NH(1)
      NW(2)        = NW(1)
      ILEFT(2)     = ILEFT(1)
      IRIGHT(2)    = IRIGHT(1)
      IOPDM(2)     = IOPDM(1)
      IOPSH(2)     = IOPSH(1)
      IOPSS(2)     = IOPSS(1)
      IOPWN(2)     = IOPWN(1)
      OPDC(2)      = OPDC(1)
      OPDS(2)      = OPDS(1)
      OPDX(2)      = OPDX(1)
      OPCMP(2)     = OPCMP(1)
      OPSPL(2)     = OPSPL(1)
           DO 200 I = 1,7
           OPDF1(I,2)  = OPDF1(I,1)
           OPDF2(I,2)  = OPDF2(I,1)
  200      CONTINUE
      RETURN
      END
      SUBROUTINE SUBWA(ITYPE)
C PURPOSE:  MAKE MICROSCOPE WAIT FOR A   GO   COMMAND BEFORE FURTHER
C  CALCULATION OR DISPLAY
      INTEGER             ITYPE
      INTEGER             OUTPUT,      LINES,       WIDTH,       ILP
      INTEGER             IDSPLA,      IPRMPT
      LOGICAL             LSCRN
      COMMON / SCREEN /   OUTPUT,      LINES,       WIDTH,       ILP
      COMMON / SCREEN /   IDSPLA,      IPRMPT,      LSCRN
      LSCRN = .FALSE.
      ITYPE = 1
      RETURN
      END
      SUBROUTINE SUBZO(ITYPE,N,IOPSS)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SUBZO
C  PURPOSE:  TO "ZOOM". THAT IS: TO MULTIPLY/DIVIDE THE APPARENT STENCIL
C             WINDOW BY A USER SUPPLIED INTEGER FACTOR N. (N>0:MULTIPLY,
C             N<0:DIVIDE).  NOTE THAT THIS COMMAND DOES NOT ALTER THE
C             ACTUAL STENCIL WIDTH.  RATHER, IT ALTERS THE STEP-SIZE
C             MULTIPLICATION FACTOR (I.E. M = IOPSS(2) ).
C
      INTEGER             IOPSS(2),    ITYPE,       N
      INTEGER             ERRCOD
      LOGICAL             ERR
      COMMON / ERRCOM /   ERRCOD,      ERR
      ITYPE = 2
      IF ( N.NE.0 ) GO TO 100
      ERRCOD = 15
      ERR    = .TRUE.
      RETURN
  100 IOPSS(2) = -N
      RETURN
      END
      SUBROUTINE SVREAD(DEVICE,NLINES,NARG,V,ERR)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  SVREAD
C  PURPOSE:  TO READ IN DOUBLE PRECISION VECTORS WITHOUT FORMATTING,
C             FROM THE SPECIFIED DEVICE.
C
      INTEGER             FINDCO,      LENGTH
      DOUBLE PRECISION    DPVAR,       V(3)
      INTEGER             BEGIN,       CHAR(72),    COUNT,       DEVICE
      INTEGER             ENDE,        I,           I72,         K
      INTEGER             NARG,        NLINES,      RIGHT
      LOGICAL             ERR
      DATA  I72 / 72 /
C
C  SET ERR = .TRUE. (I.E. NO ERRORS)
C
      ERR = .TRUE.
C
C  INITIALIZE COUNTER OF VALUES READ IN AND START RECOGNITION LOOP
C
      COUNT = 0
      NLINES= 0
 1000 READ(DEVICE,8000)(CHAR(I),I=1,72)
      DO 1500 I = 1,72
      CALL LCUC(CHAR(I))
1500  CONTINUE
      NLINES= NLINES+1
      ENDE  = LENGTH(I72,CHAR)
      IF( ENDE.EQ.0 )GO TO 1000
           BEGIN = 1
 2000 IF( BEGIN.GT.ENDE )GO TO 1000
      RIGHT  = ENDE
      K      = FINDCO(I72,CHAR,BEGIN,ENDE)-1
      IF( K.GT.0 )RIGHT = K
      CALL SRDP(I72,CHAR,BEGIN,RIGHT,DPVAR,ERR)
      IF( .NOT.ERR )GO TO 3000
      COUNT  = COUNT+1
      V(COUNT) = DPVAR
      IF( COUNT.EQ.NARG )GO TO 3000
      BEGIN  = RIGHT+2
      GO TO 2000
C
 3000 RETURN
C
C  FORMAT STATEMENT
C
 8000 FORMAT(72A1)
      END
      SUBROUTINE ZERO
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ROUTINE:  ZERO
C  PURPOSE:  TO RESET THE LOGICAL VARIABLES WHICH IMPLY A GIVEN QUANTITY
C             HAS BEEN CALCULATED.
C
      INTEGER             I,           ICENTR,      IXSIZE
      LOGICAL             LDF,         LPLT,        LDEF
      COMMON / LOGCOM /   LDF(7),      LPLT(7),     LDEF(5377)
      DOUBLE PRECISION    ROPDI,       ROPPNT,      ROPDR1,      ROPDR2
      DOUBLE PRECISION    ROPSTS,      ROPSTW,      ROPUDI
      INTEGER             NH,          NW,          ILEFT,       IRIGHT
      INTEGER             IOPDM,       IOPSH,       IOPSS,       IOPWN
      LOGICAL             OPDC,        OPDS,        OPDX,        OPCMP
      LOGICAL             OPSPL,       OPDF1,       OPDF2
      COMMON / OPTION /   ROPDI(3,2),  ROPPNT(3,2), ROPDR1(3,2)
      COMMON / OPTION /   ROPDR2(3,2), ROPSTS(2),   ROPSTW(2)
      COMMON / OPTION /   ROPUDI(3,2), NH(2),       NW(2)
      COMMON / OPTION /   ILEFT(2),    IRIGHT(2),   IOPDM(2)
      COMMON / OPTION /   IOPSH(2),    IOPSS(2),    IOPWN(2),    OPDC(2)
      COMMON / OPTION /   OPDS(2),     OPDX(2),     OPCMP(2)
      COMMON / OPTION /   OPSPL(2),    OPDF1(7,2),  OPDF2(7,2)
      DATA ICENTR,IXSIZE / 2689, 5377 /
           DO 100 I = 1,IXSIZE
           LDEF(I) = .FALSE.
  100      CONTINUE
           DO 200 I = 1,7
           LDF(I)  = .FALSE.
           LPLT(I) = .FALSE.
  200      CONTINUE
      ILEFT(2)  = ICENTR
      IRIGHT(2) = ICENTR
      OPCMP(1)  = .FALSE.
      RETURN
      END