SUBROUTINE FPOUT1(VPRNT1,DMISS,*)
C
C     FLOATING POINT OUTPUT SUBROUTINE
C
C     BOB STOUT, JULY 1972
C     MODIFIED JANUARY, 1973 BY BOB STOUT FOR DATAOUT AND XPRINT
C      REVISED SEPT. 1974
C
C     PERFORMS DYNAMIC OUTPUT FORMATTING FOR FLOATING POINT DATA.
C     ALL DATA IS PRINTED WITH THREE SIGNIFICANT FIGURES.
C     ALL SUMMARY STATISTICS ARE PRINTED WITH 4 SIGNIFICANT FIGURES,
C     EXCEPT THE SUM OF SQUARED OBSERVATIONS WHICH IS PRINTED WITH 6
C     FIGURES.
C
C     SUPERVISOR MUST INITIALIZE OVNAM IN /IO1/ WITH AN 8-CHARACTER
C     NAME FOR EACH VARIABLE WHICH MAY BE PRINTED.  NOV MUST CONTAIN
C     THE NUMBER OF POSSIBLE OUTPUT VARIABLES.	THE ORDER OF VARIABLE
C     NAMES IN OVNAM AND THE ORDER OF VALUES IN VPRNT AND VALUES
C     MUST CORRESPOND.	COSTPT MUST ALSO BE INITIALIZED--SEE BELOW
C     FOR DETAILS.
C
C------------------------------
C
C     ENTRY FPOUT1 MUST BE CALLED BY SUBROUTINE MODEL BEFORE DATA IS
C     SIMULATED FOR THE FIRST SUBJECT IN EACH EXPERIMENTAL GROUP.
C     FPOUT1 PRINTS ANY OUTPUT HEADING NECESSARY, AND PERFORMS INITIAL-
C     IZATION.	RETURN 1 FROM FPOUT1 INDICATES A FATAL ERROR SUCH AS
C     THERE BEING NO VARIABLES TO BE PRINTED.
C
C     ARGUMENTS:
C     VPRNT1 L*1  DIM(NOV).  VECTOR CONTAINING A .TRUE. VALUE FOR EACH
C		  VARIABLE WHOSE VALUE IS TO BE PRINTED FOR THE GROUP
C		  ABOUT TO BE SIMULATED, IN ANY CONDITION.  VPRNT1
C		  SHOULD BE THE INCLUSIVE OR OF ALL THE VPRNT2
C		  VECTORS ACROSS ALL CONDITIONS IN A WITHIN-
C		  SUBJECTS MULTIVARIATE DESIGN.  FOR ALL OTHER
C		  DESIGNS, VPRNT1 AND VPRNT2 MUST BE IDENTICAL.
C     DMISS  R*4  DIM(NOV).  VECTOR CONTAINING VALUES WHICH INDICATE
C		  'MISSING DATA' FOR EACH OUTPUT VARIABLE.  IF NO DATA
C		  IS EVER GOING TO BE MISSING, SOME UNLIKELY VALUE LIKE
C		  -10.E+10 SHOULD APPEAR IN THE APPROPRIATE SLOT IN
C		  DMISS.  WHEN DESCRIPTIVE STATISTICS ARE PRINTED FOR
C		  A GROUP THESE STATISTICS REFER ONLY TO THAT SUBSET
C		  OF SUBJECTS FOR WHOM COMPLETE DATA IS AVAILABLE.
C
C------------------------------
C
C     ENTRY FPOUT2 SHOULD BE CALLED BY SUBROUTINE MODEL TO PASS THE
C     VALUES OF VARIABLES TO BE PRINTED.  DATA MUST BE PASSED TO FPOUT2
C     IN A SUBJECT-BY-SUBJECT FASHION.	WHEN A WITHIN-SUBJECTS DESIGN
C     IS BEING SIMULATED, FPOUT2 MUST BE CALLED IN THE FOLLOWING
C     KIND OF SEQUENCE:
C	1ST CALL:  SUBJ 1  1ST CONDITION
C	2ND CALL:  SUBJ 1  2ND CONDITION
C	. . .
C	KTH CALL:  SUBJ 1  KTH CONDITION
C     K+1ST CALL:  SUBJ 2  1ST CONDITION
C     AND SO ON.  FPOUT2 CHECKS TO MAKE SURE THAT THE SUBJECT NUMBER
C     DOES NOT CHANGE BEFORE DATA FOR ALL CONDITIONS HAS COME IN.
C     AN ERROR MESSAGE IS PRINTED AND A RETURN 1 EXECUTED IF ANYTHING
C     GOES WRONG.
C
C     ARGUMENTS:
C     SUBJ    I*4  SUBJECT NUMBER
C     NCOND   I*4  NUMBER INDICATING WHETHER THE DATA PASSED IN VALUES
C		   IS FOR THE 1ST, 2ND, OR KTH CONDITION IN WHICH THE
C		   CURRENT SUBJECT HAS SERVED (NCOND IS IGNORED FOR
C		   BETWEEN-SUBJECTS DESIGNS).  NCOND SHOULD ***NOT***
C		   CONTAIN THE 'CONDITION NUMBER'--THE NUMBER WHICH
C		   INDICATES WHAT SET OF VARIABLE SETTINGS IS TO BE
C		   USED--NCOND SHOULD INDICATE THAT THIS IS, SAY, THE
C		   THIRD TREATMENT (CONDITION) TO WHICH THIS SUBJECT
C		   HAS BEEN SUBJECTED (THE SUBJECT MAY HAVE BEEN SUB-
C		   JECTED TO THE SAME TREATMENT THREE TIMES).
C     VPRNT2  L*1  DIM(NOV).  VECTOR CONTAINING A .TRUE. VALUE FOR
C		   EVERY VARIABLE WHOSE VALUE IS TO BE PRINTED FOR THE
C		   CONDITION CURRENTLY BEING SIMULATED.  VPRNT2 MAY BE
C		   DIFFERENT FROM ONE CONDITION TO THE NEXT IN WITHIN-
C		   SUBJECTS MULTIVARIATE EXPERIMENTS.  IN ALL OTHER
C		   CASES, VPRNT1 AND VPRNT2 MUST BE THE SAME.
C     VALUES  R*4  DIM(NOV).  VECTOR CONTAINING THE VALUES FOR EACH
C		   OUTPUT VARIABLE FOR THE GIVEN SUBJECT AND CONDITION.
C		   SEE DISCUSSION OF DMISS FOR MISSING DATA HANDLING.
C
C------------------------------
C
C     ENTRY FPOUT3 MUST BE CALLED AFTER ALL THE DATA FOR AN EXPERIMENT-
C     AL GROUP HAS BEEN PROCESSED SO THAT OUTPUT CAN BE FINISHED AND
C     STATISTICS PRINTED.
C
C     ARGUMENTS:
C     COST    R*4  SHOULD CONTAIN THE COST OF RUNNING THIS EXPER-
C		   IMENTAL GROUP.  COST IS PRINTED ONLY IF COSTPT IN
C		   /IO1/ IS .TRUE.
C
C------------------------------
C
C
      LOGICAL*4 VPRNT1(6),VPRNT2(6),OK,TRUBL,MEAN,VAR,SS,COV,CORR
      LOGICAL*4 COSTPT, DATOUT
      INTEGER*4 OVNAM,LINE,CONDS,OVN(8,6),CLABEL(2,16),OVPTRS(6)
      INTEGER*4 SUBJ,ODEV1,ODEV2,ODEV3,ODEV4,OPTR,N(16),SINK
      REAL*4 DMISS(6),VALUES(6),OUTVEC(16),SX(6,16),SXX(6,16)
      REAL*4 SXY(15,16)
C
      COMMON /IO/ IDEV(4),ODEV1,ODEV2,ODEV3,ODEV4,LINE(80)
      COMMON /IO1/ NOV,SINK,OVNAM(8,6),COSTPT,DATOUT
      COMMON /CNDTNS/ NCD,CONDS(16)
      COMMON /STAT1/ MEAN,VAR,SS,COV,CORR
C
C------------------------------
C
      OK=.TRUE.
C
C     CHECK WHICH VARIABLES HAVE TO BE PRINTED
      NVOUT=0
      DO 1 I=1,NOV
      IF(.NOT.VPRNT1(I))GO TO 1
      NVOUT=NVOUT+1
C     ADD ENTRY TO LIST OF VARIABLE LABELS
      DO 2 J=1,8
 2    OVN(J,NVOUT)=OVNAM(J,I)
      OVPTRS(NVOUT)=I
 1    CONTINUE
C
C     SET UP LIST OF CONDITION LABELS
      DO 3 I=1,NCD
      J=CONDS(I)
 3    CALL CCHARS(J,CLABEL(1,I))
C
C     SET UP TO DO STATISTICS
      DO 4 I=1,NCD
      N(I)=0
      DO 5 J=1,6
      SX(J,I)=0.
      SXX(J,I)=0.
      K1=J+1
      IF(J.EQ.6)GO TO 5
      DO 6 K=K1,6
      K2=(K-1)*(K-2)/2+J
 6    SXY(K2,I)=0.
 5    CONTINUE
 4    CONTINUE
C
C     SHOULD HAVE AT LEAST ONE DEPENDENT VARIABLE
      IF(NVOUT.GT.0)GO TO 7
      WRITE(ODEV1,1001)
      OK=.FALSE.
      RETURN 1
C
C     DETERMINE OUTPUT FORMAT
 7    IF(NCD.GT.1 .AND. NVOUT.LE.1)GO TO 10
      IF(NCD.GT.1 .AND. NVOUT.GT.1)GO TO 11
      IF(NVOUT.GT.1)GO TO 12
C
C     NO. OF CONDITIONS = 1, NO. OF VARIABLES OUTPUT = 1
      ASSIGN 110 TO IX1
      ASSIGN 502 TO IX2
      NCOND1=1
      OPTR=1
      WRITE(SINK,1002)(OVN(I,1),I=1,8)
      RETURN
C
C     NO. OF CONDITIONS > 1, NO. OF VARIABLES OUT = 1
 10   ASSIGN 210 TO IX1
      ASSIGN 503 TO IX2
      OPTR=0
      WRITE(SINK,1003)(OVN(I,1),I=1,8),((CLABEL(J,K),J=1,2),
     1 K=1,NCD)
      RETURN
C
C     NO. OF CONDITIONS = 1, NO. OF VARIABLES OUT > 1
 12   ASSIGN 310 TO IX1
      ASSIGN 504 TO IX2
      NCOND1=1
      OPTR=0
      WRITE(SINK,1004)((OVN(I,J),I=1,8),J=1,NVOUT)
      RETURN
C
C     NO. OF CONDITIONS > 1, NO. OF VARIABLES OUTPUT > 1
 11   ASSIGN 410 TO IX1
      ASSIGN 503 TO IX2
      OPTR=0
      WRITE(SINK,1005)((OVN(I,J),I=1,8),J=1,NVOUT)
      RETURN
C
C------------------------------
C
C     ENTRY TO HAVE DATA PRINTED
      ENTRY FPOUT2(SUBJ,NCOND,VPRNT2,VALUES,*)
      IF(.NOT.OK)RETURN
      TRUBL=.FALSE.
C
C     TEST FOR LEGALITY
      IF(NCOND.GT.0 .AND. NCOND.LE.NCD)GO TO 101
C     ERROR IN SUBROUTINE MODEL
      WRITE(ODEV1,1006)NCOND
      RETURN 1
C
C     BRANCH ON TYPE OF DESIGN
 101  GO TO IX1,(110,210,310,410)
C
C------------------------------
C
C     BETWEEN-SS DESIGN, ONE DEP. VAR.
 110  I=OVPTRS(1)
      IF(DATOUT)WRITE(ODEV2,9000)VALUES(I)
      OUTVEC(OPTR)=VALUES(I)
      OPTR=OPTR+1
      IF(OPTR.LE.5)GO TO 900
C
C     PRINT LINE
      WRITE(SINK,1007)(OUTVEC(J),J=1,5)
      OPTR=1
      GO TO 900
C
C------------------------------
C
C     WITHIN-SS DESIGN, ONE DEPENDENT VARIABLE
 210  IF(OPTR.EQ.0)ISUBJ=SUBJ
      NCOND1=NCOND
      IF(SUBJ.EQ.ISUBJ)GO TO 211
C
C     DID NOT GET ALL DATA FOR LAST SUBJECT
      WRITE(ODEV1,1008)ISUBJ,OPTR
      WRITE(ODEV1,1009)ISUBJ,(OUTVEC(J),J=1,OPTR)
      OPTR=0
      TRUBL=.TRUE.
      GO TO 210
C
 211  I=OVPTRS(1)
      OUTVEC(NCOND)=VALUES(I)
      OPTR=OPTR+1
      IF(OPTR.LT.NCD)GO TO 900
C
C     HAVE ALL VALUES FOR THIS SUBJECT--PRINT LINE
      WRITE(SINK,1009)SUBJ,(OUTVEC(J),J=1,NCD)
      IF(DATOUT)WRITE(ODEV2,9000)(OUTVEC(J),J=1,NCD)
      OPTR=0
      GO TO 900
C
C------------------------------
C
C     HAVE BETWEEN-SS DESIGN, SEVERAL VARIABLES
 310  DO 311 I=1,NVOUT
      J=OVPTRS(I)
 311  OUTVEC(I)=VALUES(J)
      WRITE(SINK,1009)SUBJ,(OUTVEC(J),J=1,NVOUT)
      IF(DATOUT)WRITE(ODEV2,9000)(OUTVEC(J),J=1,NVOUT)
      GO TO 900
C
C------------------------------
C
C     WITHIN-SS DESIGN, SEVERAL VARIABLES
 410  DO 411 I=1,NVOUT
      J=OVPTRS(I)
      OUTVEC(I)=VALUES(J)
      IF(.NOT.VPRNT2(J))OUTVEC(I)=-0.
 411  CONTINUE
      NCOND1=NCOND
      WRITE(SINK,1011)SUBJ,CLABEL(1,NCOND),CLABEL(2,NCOND),(OUTVEC(J),
     1 J=1,NVOUT)
      IF(DATOUT)WRITE(ODEV2,9001)NCOND,(OUTVEC(J),J=1,NVOUT)
      GO TO 900
C
C------------------------------
C
C     STATISTICS-GATHERING SECTION
 900  DO 901 I=1,NOV
      IF(.NOT.VPRNT2(I))GO TO 901
      IF(VALUES(I).EQ.DMISS(I))GO TO 910
 901  CONTINUE
      N(NCOND1)=N(NCOND1)+1
      DO 902 I=1,NOV
      IF(.NOT.VPRNT2(I))GO TO 902
      V=VALUES(I)
      SX(I,NCOND1)=SX(I,NCOND1)+V
      SXX(I,NCOND1)=SXX(I,NCOND1)+V*V
      IF(I.EQ.NOV)GO TO 902
      J=I+1
      DO 903 K=J,NOV
      IF(.NOT.VPRNT2(K))GO TO 903
      K1=(K-1)*(K-2)/2+I
      SXY(K1,NCOND1)=SXY(K1,NCOND1)+V*VALUES(K)
 903  CONTINUE
 902  CONTINUE
C
 910  IF(TRUBL)RETURN 1
      RETURN
C
C------------------------------
C
C     ENTRY TO FINISH OUTPUT AND PRINT STATISTICS
      ENTRY FPOUT3(COST)
      IF(.NOT.OK)RETURN
      GO TO IX2,(502,503,504)
C
C     FINISH PRINTING LINE OF VALUES IF NECESSARY
 502  OPTR=OPTR-1
      IF(OPTR.GT.0)WRITE(SINK,1007)(OUTVEC(J),J=1,OPTR)
      GO TO 503
C
C     MAKE SURE ALL DATA IN FOR LAST SUBJECT
 504  IF(OPTR.EQ.0)GO TO 503
      WRITE(ODEV1,1008)ISUBJ,OPTR
      WRITE(ODEV1,1009)ISUBJ,(OUTVEC(J),J=1,OPTR)
C
C     SEE IF ANY STATISTICS NEED TO BE PRINTED
 503  IF(.NOT.(MEAN.OR.SS.OR.VAR.OR.COV.OR.CORR))GO TO 550
C
C     PRINT STATISTICS CONDITION BY CONDITION
      DO 510 IC=1,NCD
      WRITE(SINK,1502)
      IF(NCD.GT.1)WRITE(SINK,1503)CLABEL(1,IC),CLABEL(2,IC),IC
      WRITE(SINK,1501)N(IC)
      IF(N(IC).LE.1)GO TO 510
      XN=N(IC)
      XN1=1./(XN-1.)
C
C     ARE ANY INDIVIDUAL VARIABLE STATISTICS TO BE PRINTED?
      IF(.NOT.(MEAN.OR.SS.OR.VAR))GO TO 511
      DO 512 I=1,NOV
      IF(.NOT.VPRNT1(I))GO TO 512
C     PRINT VARIABLE NAME
      WRITE(SINK,1504)(OVNAM(K,I),K=1,8)
      XBAR=SX(I,IC)/XN
      V=(SXX(I,IC)-XBAR*SX(I,IC))*XN1
      SD=SQRT(V)
      IF(MEAN)WRITE(SINK,1505)XBAR
      IF(VAR)WRITE(SINK,1506)V,SD
      IF(SS)WRITE(SINK,1507)SXX(I,IC)
 512  CONTINUE
C
C     PRINT CORRELATION STATISTICS AS APPROPRIATE
 511  IF(.NOT.(COV.OR.CORR) .OR. XN.LT.3. .OR. NVOUT.LE.1)GO TO 510
      I1=NOV-1
      DO 515 I=1,I1
      IF(.NOT.VPRNT1(I))GO TO 515
      XBAR=SX(I,IC)/XN
      SD=(SXX(I,IC)-XBAR*SX(I,IC))*XN1
      SD=SQRT(SD)
      J1=I+1
      DO 516 J=J1,NOV
      IF(.NOT.VPRNT1(J))GO TO 516
      WRITE(SINK,1502)
      K1=(J-1)*(J-2)/2+I
      CV=(SXY(K1,IC)-XBAR*SX(J,IC))*XN1
      IF(COV)WRITE(SINK,1509)(OVNAM(K,I),K=1,8),(OVNAM(K1,J),K1=
     1 1,8),CV
      SD2=(SXX(J,IC)-SX(J,IC)*SX(J,IC)/XN)*XN1
      SD2=SQRT(SD2)*SD
      R=0.
      IF(SD2.GT.10.0E-10)R=CV/SD2
      IF(CORR)WRITE(SINK,1510)(OVNAM(K,I),K=1,8),(OVNAM(K1,J),K1=
     1	1,8),R
 516  CONTINUE
 515  CONTINUE
 510  CONTINUE
C
C     QUIT
 550  IF(COSTPT)WRITE(SINK,1512)COST
      RETURN
C
C-----------------------------
C
 1001 FORMAT(1H0,'NO DEPENDENT VARIABLE?')
 1002 FORMAT(1H0,4X,8A1,' SCORES')
 1003 FORMAT(1H0,10X,8A1,' SCORES'/1X,'SUBJ  COND:',2X,6(2A1,8X)/
     1 15X,6(2A1,8X)/15X,6(2A1,8X))
 1004 FORMAT(1H0,'SUBJ',4X,6(2X,8A1))
 1005 FORMAT(1H0,'SUBJ COND',6(2X,8A1))
 1006 FORMAT(1H0,'BAD CALL TO FPOUT2; NCOND=',I4)
 1007 FORMAT(10X,5G10.3)
 1008 FORMAT(1H0,'DATA INCOMPLETE FOR SUBJ',I4,' # VALUES=',I3)
 1009 FORMAT(1X,I4,5X,6G10.3,(/11X,6G10.3))
 1011 FORMAT(1X,I4,3X,2A1,1X,6G10.3)
C
 1501 FORMAT(5X,'NO. OF SS WITH COMPLETE DATA:',I5)
 1502 FORMAT('	')
 1503 FORMAT(5X,'CONDITION ',2A1,' (',I2,')')
 1504 FORMAT(1H0,4X,'VARIABLE: ',8A1)
 1505 FORMAT(5X,'MEAN: ',G10.4)
 1506 FORMAT(5X,'VARIANCE: ',G10.4/5X,'STD. DEVIATION: ',G10.4)
 1507 FORMAT(5X,'SUM OF SQUARED OBSERVATIONS: ',G15.6)
 1509 FORMAT(5X,'COVARIANCE  OF ',8A1,' AND ',8A1,' IS: ',G12.4)
 1510 FORMAT(5X,'CORRELATION OF ',8A1,' AND ',8A1,' IS:',F8.4)
 1512 FORMAT(1H0,4X,'COST: ',G12.4)
 9000 FORMAT(8G10.3)
C     FMT 9001 IS NECESSITATED BY INCOMPETENT IMPLEMENTATION OF
C     G FORMATS IN IBM FORTRAN IV.
 9001 FORMAT(I10,7G10.3/(8G10.3))
C
      END
      SUBROUTINE URAND(X)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS UNIFORMLY DISTRI-
C     BUTED OVER THE INTERVAL (0,1)
C
      X=RAN(DUMMY)
      RETURN
      END
      SUBROUTINE NRAND(MEAN,SIGMA,X)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS WITH AN APPROXIMATE
C     NORMAL DISTRIBUTION WITH LOCATION PARAMETER MEAN AND STANDARD DEVIA-
C     TION SIGMA.  VALUE GENERATED IS RETURNED IN X.
C     REQUIRES SUBROUTINE URAND.
C
      REAL MEAN
C
C
      X=0.
      DO 1 I=1,12
      CALL URAND(Y)
 1    X=X+Y
      X=X-6.
      X=X*SIGMA
      X=X+MEAN
      RETURN
      END
      SUBROUTINE URAND1(XL1,XL2,X)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS DISTRIBUTED UNIFORMLY
C     OVER THE INTERVAL (XL1,XL2).  XL1 AND XL2 DO NOT HAVE TO BE
C     PROPERLY ORDERED.
C     REQUIRES SUBROUTINE URAND.
C
C
      X1=AMIN1(XL1,XL2)
      DIFF=XL2-XL1
      DIFF=ABS(DIFF)
      CALL URAND(X)
      X=DIFF*X
      X=X+X1
      RETURN
      END
      SUBROUTINE URAND2(L1,L2,IX)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM INTEGERS DISTRIBUTED
C     UNIFORMLY OVER THE INTERVAL (L1,L2).  L1 AND L2 DO NOT HAVE
C     TO BE PROPERLY ORDERED.
C     REQUIRES SUBROUTINES URAND1 AND URAND.
C
C
      L=MIN0(L1,L2)
      IDIFF=L2-L1
      IDIFF=IABS(IDIFF)
      DIFF=IDIFF+1
      CALL URAND1(0.,DIFF,X)
      IX=X
      IF(IX.GE.(IDIFF+1))IX=IDIFF
      IX=L+IX
      RETURN
      END
      SUBROUTINE BINOM(N,P,IX)
C
C     SUBROUTINE TO GENERATE BINOMIAL RANDOM VARIABLES
C     VALUES ARE RETURNED IN IX.
C     NORMAL APPROXIMATION USED IF N*P>15
C     POISSON APPROXIMATION USED IF P<.1 AND N>25
C     OTHERWISE, NUMBERS GENERATED BY N DRAWS
C     REQUIRES SUBROUTINE URAND
C
C
      IX=0
      IF(N.EQ.0 .OR. P.EQ.0.)RETURN
      IF(N.LT.0 .OR. P.LT.-0.00001 .OR. P.GT.1.00001)STOP '3901'
C
      P1=1.-P
      P2=AMIN1(P,P1)
      XN=N
      IF((P2*XN).GT.15.)GO TO 1
      IF(P2.LT.0.1 .AND. N.GT.25)GO TO 3
C
C     USE REAL LIVE BINOMIAL PROCESS TO GENERATE NUMBERS
      DO 2 I=1,N
      CALL URAND(X)
      IF(X.LE.P)IX=IX+1
 2    CONTINUE
      RETURN
C
C     USE NORMAL APPROXIMATION
 1    XM=P*XN
      S=P*P1*XN
      S=SQRT(S)
      CALL NRAND(XM,S,X)
      X=X+.5
      IX=X
      IF(IX.LT.0)IX=0
      IF(IX.GT.N)IX=N
      RETURN
C
C     USE POISSON APPROXIMATION
 3    XM=P2*XN
      CALL POISSN(XM,IX)
      IF(IX.GT.N)IX=N
      IF(P.EQ.P2)RETURN
      IX=N-IX
      RETURN
      END
      SUBROUTINE POISSN(LAMBDA,IX)
C
C     SUBROUTINE TO GENERATE RANDOM INTEGER NUMBERS WITH A POISSON
C     DISTRIBUTION WITH INTENSITY PARAMETER LAMBDA
C     REQUIRES SUBROUTINES URAND AND NRAND
C     NORMAL APPROXIMATION IS USED FOR LAMBDA>16
C
C
      REAL*4 LAMBDA
C
C     TAKE CARE OF PROGRAMMING ERRORS, EXTREME CASES
      IF(LAMBDA.LT.0.)STOP '3903'
      IX=0
      IF(LAMBDA.LT.1.0E-8)RETURN
C
      IF(LAMBDA.GT.16.)GO TO 1
C
C     DO IT THE HARD WAY
      CALL URAND(UNIFRV)
      P=EXP(-LAMBDA)
      CUMP=P
      X=0.
      DO 2 I=1,100
      X=X+1.
      IF(UNIFRV.LE.CUMP)GO TO 3
      P=P*LAMBDA/X
 2    CUMP=CUMP+P
C
C     IF WE EVER FINISH THE DO LOOP, THERE IS PROBABLY SOMETHING WRONG
C
 3    IX=I-1
      RETURN
C
C     USE NORMAL APPROXIMATION
 1    P=SQRT(LAMBDA)
      CALL NRAND(LAMBDA,P,X)
      IX=X+.5
      IF(IX.LT.0)IX=0
      RETURN
      END
      SUBROUTINE MULNOM(PTABLE,N1,N2,N3,N4,IV1,IV2,IV3,IV4)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS WITH A MULTIVARIATE
C     MULTINOMIAL DISTRIBUTION SPECIFIED BY A TABLE OF PROBABILITIES.
C     1-4 VARIABLES MAY BE GENERATED.  REQUIRES SUBROUTINE URAND.
C
C
C     ARGUMENTS:
C     PTABLE  REAL*4  DIMENSIONS (N1,N2,N3,N4)
C     TABLE OF PROBABILITIES.  PTABLE(I,J,K,L) SHOULD CONTAIN THE PROB-
C     ABILITY THAT VARIABLE 1 = I, VARIABLE 2 = J, ETC.
C
C     N1, N2, N3, N4  INTEGER*4
C     DIMENSIONS OF PTABLE.  FOR TABLES OF M DIMENSIONS (M<4), USE N1
C     THRU NM AND SET N(M+1) THRU N4 TO 1.
C
C     IV1, IV2, IV3, IV4  INTEGER*4
C     THE VALUES GENERATED ARE RETURNED IN THESE LOCATIONS.
C
C
C     ERROR HALTS:
C     STOP 3904  ONE OF N1-N4 IS LESS THAN 1
C     STOP 3905  A VALUE IN PTABLE IS LESS THAN -.000001 OR GREATER
C		 THAN 1.000001
C     STOP 3906  THE VALUES IN PTABLE SUM TO LESS THAN .999--THIS COND-
C		 ITION IS NOT NECESSARILY DISCOVERED ON THE FIRST CALL
C		 TO MULNOM
C
C------------------------------
C
      REAL*4 PTABLE(N1,N2,N3,N4)
C
C     ERROR TESTING
      IF(N1.LT.1 .OR. N2.LT.1 .OR. N3.LT.1 .OR. N4.LT.1)STOP '3904'
C
      CALL URAND(X)
      P=0.
      DO 1 IV4=1,N4
      DO 1 IV3=1,N3
      DO 1 IV2=1,N2
      DO 1 IV1=1,N1
      P1=PTABLE(IV1,IV2,IV3,IV4)
      IF(P1.LT.-0.000001 .OR. P1.GT.1.000001)STOP '3905'
      P=P+P1
      IF(X.LE.P) RETURN
 1    CONTINUE
      IF(P.LT.0.999)STOP '3906'
      RETURN
      END
      SUBROUTINE CMULNM(PTABLE,N1,N2,N3,N4,IV1,IV2,IV3,IV4)
C
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS WITH A COND-
C     ITIONAL MULTIVARIATE MULTINOMIAL DISTRIBUTION AS SPECIFIED BY A
C     TABLE OF PROBABILITIES.  CAN HANDLE UP TO 4 VARIABLES; VALUES
C     RETURNED MAY BE CONDITIONED ON ANY COMBINATION OF 0 TO 4 OF THE
C     VARIABLES.  REQUIRES SUBROUTINE URAND.
C
C
C     ARGUMENTS:
C     PTABLE  REAL*4  DIMENSIONS (N1,N2,N3,N4)
C     TABLE OF PROBABILITIES.  PTABLE(I,J,K,L) SHOULD CONTAIN THE
C     PROBABILITY THAT VARIABLE 1 = I, VARIABLE 2 = J, ETC.
C
C     N1, N2, N3, N4  INTEGER*4
C     DIMENSIONS OF PTABLE.  FOR TABLES OF M DIMENSIONS (M<4) USE N1
C     THRU NM AND SET N(M+1) THRU N4 TO 1.
C
C     IV1, IV2, IV3, IV4  INTEGER*4
C     VALUES GENERATED ARE RETURNED IN THESE LOCATIONS.  THE VALUES IN
C     IV1-IV4 AT THE TIME CMULNM IS CALLED INDICATE WHICH VARIABLES
C     ARE TO BE ASSIGNED VALUES BY THE PROGRAM AND WHICH ARE TO COND-
C     ITION THE SELECTION OF THE OTHERS.  IF IVJ IS GREATER THAN 0,
C     ITS VALUE IS USED TO CONDITION THE SELECTION OF THE VARIABLES
C     WHOSE VALUES ARE 0 OR LESS.
C
C
C     ERROR HALTS:
C     STOP 3907  ONE OF N1-N4 IS LESS THAN 1
C     STOP 3908  A VALUE IN PTABLE IS LESS THAN -0.000001 OR GREATER
C		 THAN 1.000001
C     STOP 3909  THE VALUES IN PTABLE SUM TO MORE THAN 1.001--THIS
C		 CONDITION IS NOT NECESSARILY DISCOVERED ON THE FIRST
C		 CALL TO CMULNM
C     STOP 3910  THE PROBABILITIES IN ONE HYPERPLANE OF PTABLE SUM TO
C		 LESS THAN 10**-20.  THIS CONDITION IS NOT NECESSARILY
C		 DISCOVERED ON THE FIRST CALL
C
C--------------------------
      REAL*4 PTABLE(N1,N2,N3,N4)
C
C     ERROR TESTING
      IF(N1.LT.1 .OR. N2.LT.1 .OR. N3.LT.1 .OR. N4.LT.1)STOP '3907'
C
      IF(IV1.GT.0 .AND. IV2.GT.0 .AND. IV3.GT.0 .AND. IV4.GT.0)RETURN
C
C     COMPUTE MULTIPLIER FACTOR
      I1A=1
      I1B=N1
      IF(IV1.LT.1)GO TO 1
      I1A=IV1
      I1B=IV1
 1    I2A=1
      I2B=N2
      IF(IV2.LT.1)GO TO 2
      I2A=IV2
      I2B=IV2
 2    I3A=1
      I3B=N3
      IF(IV3.LT.1)GO TO 3
      I3A=IV3
      I3B=IV3
 3    I4A=1
      I4B=N4
      IF(IV4.LT.1)GO TO 4
      I4A=IV4
      I4B=IV4
C
 4    PT=0.
      DO 5 L=I4A,I4B
      DO 5 K=I3A,I3B
      DO 5 J=I2A,I2B
      DO 5 I=I1A,I1B
      P1=PTABLE(I,J,K,L)
      IF(P1.LT.-0.000001 .OR. P1.GT.1.000001)STOP '3908'
      PT=PT+P1
 5    CONTINUE
      IF(PT.GT.1.001)STOP '3909'
      IF(PT.LT.10.0E-20)STOP '3910'
      PT=1./PT
C
C     SELECT VALUES FOR THE UNSPECIFIED VARIABLES
      CALL URAND(X)
      P=0.
      DO 6 IV4=I4A,I4B
      DO 6 IV3=I3A,I3B
      DO 6 IV2=I2A,I2B
      DO 6 IV1=I1A,I1B
      P=P+PT*PTABLE(IV1,IV2,IV3,IV4)
      IF(X.LE.P)RETURN
 6    CONTINUE
      RETURN
      END
      SUBROUTINE REXP(BETA,VALUE)
C     SUBROUTINE TO GENERATE RANDOM NUMBERS WITH AN EXPONENTIAL DIST-
C     RIBUTION WITH SCALE PARAMETER BETA.  FOR STANDARD EXPONENTIALS,
C     SET BETA=1.  IF BETA IS LESS THAN -.00001, AN ERROR HALT
C     (STOP 3620) OCCURS.  THE RANDOM NUMBER GENERATED IS RETURNED
C     IN VALUE.
C
C
C
C
      IF(BETA .LT. -.00001)STOP 3620
      IF(BETA .LT. 0.)BETA=0.
C
C     EXPONENTIALS ARE EASY
      CALL URAND(X)
      VALUE=-BETA*ALOG(X)
      RETURN
      END
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS HAVING A
C     DOUBLE EXPONENTIAL DISTRIBUTION WITH LOCATION PARAMETER
C     ALPHA AND SCALE PARAMETER BETA.  IF BETA < -0.00001 AN
C     ERROR HALT (STOP 3617) OCCURS.  THE RANDOM NUMBER GEN-
C     ERATED IS RETURNED IN VALUE.
C
C------------------------------
C
      SUBROUTINE RDEXP(ALPHA, BETA, VALUE)
C
C
      IF(BETA .LT. -0.00001)STOP 3617
C
      CALL URAND(X)
      VALUE=-BETA*ALOG(X)
      CALL URAND(X)
      IF(X .GT. 0.5)VALUE=-VALUE
      VALUE=VALUE+ALPHA
      RETURN
      END
      SUBROUTINE RGAMMA(IALPHA, BETA, VALUE)
C     SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS HAVING A
C     GAMMA DISTRIBUTION WITH SHAPE PARAMETER IALPHA AND SCALE
C     PARAMETER BETA.  THE RANDOM NUMBER GENERATED IS RETURNED
C     IN VALUE.  IF IALPHA IS LESS THAN OR EQUAL TO ZERO, AN
C     ERROR HALT (STOP 3621) OCCURS.  IF BETA IS LESS THAN -.00001
C     A STOP 3622 OCCURS.
C
C     NOTE: FOR IALPHA < 31 THE GENERATION PROCESS INVOLVES
C     GENERATING UP TO 30 RANDOM EXPONENTIALS, SO FOR 5<IALPHA<31
C     THIS ALGORITHM IS RATHER SLOW, ALTHOUGH QUITE ACCURATE.
C     FOR IALPHA GREATER THAN OR EQUAL TO 31, A NORMAL APPROXI-
C     MATION IS USED WHICH IS FAST BUT DOES NOT HAVE TAILS OF
C     THE PROPER SHAPE, SO IF THE SHAPE OF THE TAILS IS CRITICAL
C     THIS ALGORITHM SHOULD NOT BE USED FOR IALPHA GREATER THAN
C     OR EQUAL TO 31.
C
C
      INTEGER IALPHA
      REAL  BETA, GAMMA
      REAL*8 C1, E1, PI
	GAMMA=GAMMA!DEC-10 BUG REQUIRES THIS
	E1=E1!DITTO
	C1=C1!DEC-10 BUG
	PI=PI!DITTO
C
C     TRAP ERRORS
      IF(IALPHA .LE. 0)STOP 3621
      IF(BETA .LT. -.00001)STOP 3622
      IF(BETA .LT. 0.)BETA=0.
C
      IF(IALPHA .GE. 31)GO TO 20
C
C     GENERATE VALUE BY SUM OF EXPONENTIALS
      VALUE=0.
C
C     ADD APPROPRIATE NO. OF STANDARD RANDOM EXPONENTIALS
      DO 10 J=1,IALPHA
      CALL URAND(X)
      X=-ALOG(X)
 10   VALUE=VALUE+X
      VALUE=VALUE*BETA
      RETURN
C
C     USE APPROXIMATION FOR LARGE VALUES OF SHAPE PARAMETER
 20   X1=IALPHA
      X1=X1*BETA
      X2=X1*BETA
      X2=SQRT(X2)
 21   CALL NRAND(X1, X2, VALUE)
      IF(VALUE .LE. 0.)GO TO 21
      RETURN
      END
      SUBROUTINE PLFN(TABLE,N,XIN,YOUT)
C     SUBROUTINE TO CALCULATE A PIECEWISE LINEAR FUNCTION OF A
C     CONTINUOUS-VALUED VARIABLE.
C
C     ARGUMENTS:
C
C     TABLE  REAL*4
C     THE DIMENSIONS OF TABLE MUST BE (N,2), WHERE N IS THE NUMBER OF
C     GRID POINTS.  XIN IS COMPARED AGAINST THE VALUES IN TABLE(1-N,1)
C     AND THE VALUE OF YOUT IS CALCULATED BY LINEAR INTERPOLATION
C     USING THE APPROPRIATE VALUES IN TABLE(1-N,2).  NOTE: NUMBERS
C     IN THE UPPER HALF OF TABLE (TABLE(1-N,1)) MUST INCREASE FROM
C     LEFT TO RIGHT; THAT IS, TABLE(1,1) MUST BE LESS THAN TABLE(2,1)
C     AND SO ON.  THERE IS, OF COURSE, NO SUCH RESTRICTION ON THE
C     NUMBERS IN THE LOWER HALF OF TABLE (TABLE(1-N,2)).  THE VALUES
C     IN THE UPPER HALF OF TABLE DO NOT HAVE TO BE EQUALLY SPACED.
C
C     N  INTEGER*4
C     FIRST DIMENSION OF TABLE; SEE ABOVE.  IF N IS LESS THAN 2
C     AN ERROR HALT (STOP 3701) OCCURS.
C
C     XIN  REAL*4
C     X-VALUE FOR WHICH A CORRESPONDING Y-VALUE IS TO BE CALCULATED.
C     VALUES OF XIN WHICH ARE LESS THAN THE VALUE IN TABLE(1,1) ARE
C     TREATED AS THOUGH THEY WERE EQUAL TO TABLE(1,1); SIMILARLY,
C     VALUES OF XIN WHICH ARE GREATER THAN THE VALUE IN TABLE(N,1)
C     ARE TREATED AS THOUGH THEY WERE EQUAL TO TABLE(N,1).  THIS
C     MEANS THAT THE Y-VALUES IN THE ENDS OF THE TABLE ARE TREATED
C     BY THE PROGRAM AS ASYMPTOTIC VALUES OF THE FUNCTION.
C
C     YOUT  REAL*4
C     CONTAINS ON RETURN THE ESTIMATED VALUE OF THE FUNCTION, THE
C     RESULT OF THE CALCULATION.
C
C
      REAL TABLE(N,2)
C
C
      IF(N.LE.1)STOP 3701
      X=XIN
      IF(X.LT.TABLE(1,1))X=TABLE(1,1)
C
      DO 1 I=2,N
      IF(X.LE.TABLE(I,1))GO TO 2
 1    CONTINUE
      X=TABLE(N,1)
 2    D=(X-TABLE(I-1,1))/(TABLE(I,1)-TABLE(I-1,1))
      YOUT=D*(TABLE(I,2)-TABLE(I-1,2))+TABLE(I-1,2)
      RETURN
      END
      SUBROUTINE INTPOL(TABLE,N1,N2,N3,N4,X1,X2,X3,X4,VALUE)
C     SUBROUTINE TO PERFORM LINEAR INTERPOLATION IN TABLES OF UP
C     TO 4 DIMENSIONS.
C
C     ARGUMENTS:
C
C     TABLE  REAL*4
C     TABLE IN WHICH INTERPOLATING IS TO BE DONE.  TABLE MUST BE DIM-
C     ENSIONED (N1,N2,N3,N4).
C
C     N1,N2,N3,N4  INTEGER*4
C     DIMENSIONS OF TABLE.  IF ANY OF N1-4 IS LESS THAN 1 OR GREATER
C     THAN 10000 AN ERROR MESSAGE WILL BE PRINTED AND THE PROGRAM WILL
C     HALT (STOP 3700).
C
C     X1,X2,X3,X4  REAL*4
C     X1-4 ARE THE SUBSCRIPT VALUES TO BE USED IN CALCULATING A
C     VALUE FROM THE TABLE.  *****ALWAYS REMEMBER THAT X1-4 ARE
C     REAL NUMBERS, NOT INTEGERS*****  IF ANY OF X1-4 EXCEEDS
C     100,000 IN ABSOLUTE VALUE, AN ERROR MESSAGE WILL BE PRINTED
C     AND THE PROGRAM WILL HALT (STOP 3700).  OTHERWISE, HOWEVER,
C     AN X-VALUE LESS THAN 1 IS TREATED AS BEING EQUAL TO 1, AND
C     A VALUE FOR ANY XJ WHICH IS GREATER THAN THE CORRESPONDING
C     NJ IS TREATED AS EQUAL TO NJ.  THUS, THE VALUES IN THE
C     EXTREME PARTS OF THE TABLE ARE TREATED AS ASYMPTOTIC VALUES.
C
C     VALUE  REAL*4
C     THE NUMBER CALCULATED BY INTERPOLATION IS RETURNED IN VALUE.
C
C------------------------------
C
      REAL TABLE(N1,N2,N3,N4), Y(4), T(2,2,2)
      INTEGER M(4), L(4,2), ODEV1
      LOGICAL*4 GIVEUP
      COMMON /IO/IDEV(4),ODEV1
C
C------------------------------
C
C     ISN'T FORTRAN AWFUL?
      M(1)=N1
      M(2)=N2
      M(3)=N3
      M(4)=N4
      Y(1)=X1
      Y(2)=X2
      Y(3)=X3
      Y(4)=X4
C
C     ERROR CHECKING + OTHER TASKS
      GIVEUP=.FALSE.
      DO 1 I=1,4
      IF(M(I).GT.0 .AND. M(I).LE.10000)GO TO 2
C
C     INVALID DIMENSION
      WRITE(ODEV1,1001)I,M(I)
      GIVEUP=.TRUE.
C
 2    IF(ABS(Y(I)).LE.100000.)GO TO 3
C
C     INVALID SUBSCRIPT
      WRITE(ODEV1,1002)I,Y(I)
      GIVEUP=.TRUE.
      GO TO 1
C
C     TRUNCATION, ETC.
 3    IF(Y(I).LT.1.)Y(I)=1.
      Z=M(I)
      IF(Y(I).GT.Z)Y(I)=Z
      L(I,1)=Y(I)+0.000001
      L(I,2)=L(I,1)+1
      IF(L(I,2).GT.M(I))L(I,2)=M(I)
 1    CONTINUE
      IF(GIVEUP)STOP 3700
C
C------------------------------
C
C     COLLAPSE TABLE STARTING WITH 4RTH DIMENSION
      D=L(4,1)
      D=Y(4)-D
      L1=L(4,1)
      L2=L(4,2)
      DO 14 K=1,2
      K1=L(3,K)
      DO 14 J=1,2
      J1=L(2,J)
      DO 14 I=1,2
      I1=L(1,I)
 14   T(I,J,K)=D*(TABLE(I1,J1,K1,L2)-TABLE(I1,J1,K1,L1))
     1	+TABLE(I1,J1,K1,L1)
C     AREN'T YOU GLAD THAT'S OVER?
C
C     3RD DIMENSION
      D=L(3,1)
      D=Y(3)-D
      DO 13 I=1,2
      DO 13 J=1,2
 13   T(I,J,1)=D*(T(I,J,2)-T(I,J,1))+T(I,J,1)
C
C     2ND DIMENSION
      D=L(2,1)
      D=Y(2)-D
      DO 12 I=1,2
 12   T(I,1,1)=D*(T(I,2,1)-T(I,1,1))+T(I,1,1)
C
C     AND THE FIRST DIMENSION SHALL BE THE LAST
      D=L(1,1)
      D=Y(1)-D
      VALUE=D*(T(2,1,1)-T(1,1,1))+T(1,1,1)
      RETURN
C
C
 1001 FORMAT(1H0,'INTPOL ERROR--INVALID DIMENSION (',I1,'): ',I8)
 1002 FORMAT(1H0,'INTPOL ERROR--INVALID SUBSCRIPT FOR DIM. ',
     1	I1,': ',G15.6)
      END
      SUBROUTINE CHANCE(P, *)
C     SUBROUTINE TO SELECT ONE OF TWO ALTERNATIVE PATHS PROB-
C     ABILISTICALLY
C
C     P IS THE PROBABILITY THAT THE PROGRAM SHOULD BRANCH TO THE
C     STATEMENT SPECIFIED BY THE STATEMENT NUMBER IN THE CALLING
C     SEQUENCE; WITH PROBABILITY 1-P THE PROGRAM WILL EXECUTE THE
C     STATEMENT FOLLOWING THE CALL.  FOR EXAMPLE, THE CODE BELOW
C     WOULD CAUSE THE PROGRAM TO BRANCH TO THE STATEMENT LABELED
C     100 WITH PROBABILITY .8, AND WOULD PROCEED INSTEAD TO
C     STATEMENT 200 AFTER THE CALL WITH PROBABILITY .2.
C
C	   CALL CHANCE(.8, &100)
C      200 . . .
C	   . . .
C      100 . . .
C
C
C     IF P IS OUTSIDE THE RANGE (-.00001, 1.00001) AN ERROR HALT
C     (STOP 3603) WILL OCCUR.
C
C------------------------------
C
C
      IF(P .LT. -.00001 .OR. P .GT. 1.00001)STOP 3603
      CALL URAND(X)
      IF(X .LE. P)RETURN 1
      RETURN
      END
      SUBROUTINE ROUND(X, UNIT)
C     SUBROUTINE TO ROUND FLOATING POINT NUMBERS TO ANY SPECIFIED
C     DEGREE OF PRECISION.  THE NUMBER TO BE ROUNDED SHOULD BE IN X,
C     AND THE MAGNITUDE OF THE UNIT TO BE USED IN ROUNDING X SHOULD
C     BE IN UNIT.  THUS, IF YOU WANT X TO VARY IN STEPS OF .25, SET
C     UNIT=.25; FOR STEPS OF 25, SET UNIT=25., AND SO ON.  THE
C     ROUNDED NUMBER IS RETURNED IN X.	IF UNIT IS LESS THAN OR EQUAL
C     TO 10**-20, AN ERROR HALT (STOP 3704) WILL OCCUR.
C
C
      REAL X, UNIT
C
C     PREVENT FLOATING POINT OVERFLOWS, DIVIDING BY 0, ETC.
      IF(UNIT .LE. 1.E-20)STOP 3704
      Y=X/UNIT+.5
      Y=AINT(Y)
      X=Y*UNIT
      RETURN
      END
      SUBROUTINE STETTR(DFN,DX,DY,DXN,DH)
C     STETTR
C
C     NUMERICAL INTEGRATION SUBROUTINE FOR THE SOLUTION OF
C     Y'=DFN(DX,DY) USING A TWO STEP ALGORITHM DUE TO HANS J. STETTER.
C
C     THE ORDER FOUR RUNGE-KUTTA METHOD IS USED FOR OBTAINING INITIAL
C     VALUES.  THE STETTER ALGORITHM IS ALSO OF ORDER FOUR.  THE
C     PRINCIPAL ADVANTAGE OF THE STETTER ALGORITHM IS THAT IT IS
C     ROUGHLY TWICE AS FAST AS THE ORDER 4 RUNGE-KUTTA METHOD.
C
C     REFERENCE:
C	STETTER, H. J.	STABILIZING PREDICTORS FOR WEAKLY UNSTABLE
C	CORRECTORS.  MATHEMATICS OF COMPUTATION, 1965, VOL. 19,
C	84-89.
C
C
C     EXECUTION IS TERMINATED IF ERRORS IN THE CALLING PARAMETERS ARE
C     DETECTED.
C
C     ROBERT L. STOUT
C     DEPARTMENT OF PSYCHOLOGY
C     UNIVERSITY OF MICHIGAN
C     NOV. 25, 1968
C     MODIFIED MAY 23, 1973
C     MODIFIED SEPT. 23, 1973
C
C
C     ARGUMENTS (ALL IN DOUBLE PRECISION):
C	DFN  FUNCTION BEING INTEGRATED (MUST BE A DOUBLE PRECISION
C	     FUNCTION)
C	DX   X-VALUE FROM WHICH INTEGRATION IS TO START (CONTAINS
C	     X-VALUE AT WHICH INTEGRATION STOPPED ON RETURN)
C	DY   INITIAL VALUE OF THE INTEGRAL (CONTAINS THE FINAL
C	     VALUE OF THE INTEGRAL ON RETURN)
C	DXN  X-VALUE AT WHICH INTEGRATION IS TO STOP
C	DH   INTEGRATION STEPSIZE
C
C------------------------------
C
C
      IMPLICIT REAL*8(D)
      EXTERNAL DFN
C
C     CHECK LEGALITY
      IF(DH .LE. 0.D0 .OR. DXN .LE. DX)GO TO 1
C
C     INITIALIZATION
      DF0=DFN(DX,DY)
      DY0=DY
      DX1=DX+DH
      CALL RUNGK(DFN,DX,DY,DX1,DH)
      DF1=DFN(DX1,DY)
      DY1=DY
      DX=DX1
      D1=.5D0*DH
      D2=2.D0*DH
      D3=DH/3.0D0
C
C     MAIN INTEGRATION LOOP
 2    DX=DX+DH
C     ARE WE DONE?
      IF(DX.GE.(DXN+D1))GO TO 3
C     PROTECT THE CUSTOMER FROM POLES (NOT THE EAST EUROPEAN KIND)
      IF(DX.GT.DXN)DX=DXN
C
C     PREDICT
      DY2=-4.0D0*DY1+5.0D0*DY0+D2*(2.0D0*DF1+DF0)
C
C     CORRECT
      DF2=DFN(DX,DY2)
      DY2=DY0+D3*(DF2+4.0D0*DF1+DF0)
C
C     PUSH DOWN ARGUMENTS
      DY0=DY1
      DY1=DY2
      DF0=DF1
      DF1=DFN(DX,DY2)
      GO TO 2
C
C
C     DONE
 3    DY=DY1
      RETURN
C
C
C     ERROR CONDITION
 1    WRITE(6,1001)DH,DX,DY,DXN
      STOP 3703
C
C
 1001 FORMAT(1H0,9X,'ERROR CONDITION IN STETTR'/10X,'DH=',G15.8/10X,
     1'DX=',G15.8/10X,'DY=',G15.8/10X,'DXN=',G15.8)
      END
      SUBROUTINE RUNGK(DFN,DX,DY,DXN,DH)
C
C     FOURTH ORDER RUNGE-KUTTA INTEGRATION SUBROUTINE
C
C     DFN IS THE FUNCTION NAME, DX AND DY ITS ARGUMENTS.  ALL NAMES
C     BEGINNING WITH D ARE DOUBLE PRECISION NAMES.
C     DX AND DY MUST BE INITIALIZED BY THE CALLING PROGRAM.  DXN
C     MUST ALSO BE SET TO THE TERMINATING X VALUE, AND DH TO THE
C     STEPSIZE.
C
C     ROBERT L. STOUT
C     DEPARTMENT OF PSYCHOLOGY
C     UNIVERSITY OF MICHIGAN
C     NOVEMBER 25, 1968
C     REVISED MAY 23, 1973
C     REVISED SEPT. 23, 1973
C
C------------------------------
C
      IMPLICIT REAL*8(D)
C
C     CHECK LEGALITY
      IF(DH .LE. 0.0D0 .OR. DXN .LE. DX)GO TO 1
C
C     INITIALIZE
      DK1=DFN(DX,DY)
      DH2=.5D0*DH
      DZ=DH/6.0D0
C
C     INTEGRATE
 2    IF(DX.GE.(DXN+DH2))RETURN
      DA=DX+DH2
      DB=DY+DH2*DK1
      DK2=DFN(DA,DB)
      DB=DY+DH2*DK2
      DK3=DFN(DA,DB)
      DX=DX+DH
      DB=DY+DH*DK3
      DK4=DFN(DX,DB)
      DY=DY+DZ*(DK1+2.0D0*DK2+2.0D0*DK3+DK4)
      DK1=DFN(DX,DY)
      GO TO 2
C
C     BAD CALL
   1  WRITE(6,1001)DX,DY,DXN,DH
      STOP 3702
C
C
 1001 FORMAT(1H0,'BAD CALL TO RUNGK'/1X,'DX=',G15.8,2X,'DY=',G15.8,2X,
     1 'DXN=',G15.8,2X,'DH=',G15.8)
      END
      SUBROUTINE NEXT(LPTR,NAME,BC,NUM,FNUM,SFLAG,*,*,*)
C
C     SUBROUTINE TO GET NEXT INTERESTING THING FROM LINE STARTING
C     AT LPTR
C
C     A SEPARATE SUBROUTINE RETURN IS PROVIDED FOR EACH CLASS
C     OF INTERESTING THINGS:
C	RETURN 0  ALPHAMERIC NAME (1-8 CHARACTERS)
C	RETURN 1  NUMBER
C	RETURN 2  NON-BLANK BREAK CHARACTER
C	RETURN 3  END OF LINE
C
C     ARGUMENTS:
C	LPTR  INTEGER*4  POINTER TO THE NEXT CHARACTER TO BE
C			 PROCESSED BY NEXT.  POINTS TO CHARACTER
C			 AFTER THE LAST ONE PROCESSED BY NEXT
C			 ON RETURN.
C	NAME  INTEGER*4  8 CHARACTER VECTOR IN WHICH ALPHA-
C			 MERIC NAMES CAN BE RETURNED.
C	BC    INTEGER*4  CONTAINS BREAK CHARACTER AFTER RETURN 2.
C	NUM   INTEGER*4  CONTAINS ROUNDED INTEGER VALUE OF NUMBER
C			 AFTER RETURN 1.
C	FNUM  REAL*4	 CONTAINS FLOATING POINT VALUE FOR NUMBER
C			 AFTER RETURN 1.
C
C     ALPHAMERIC NAMES 9 OR MORE CHARACTERS LONG ARE TRUNCATED TO
C     8 CHARACTERS WITH A WARNING MESSAGE
C
C
      IMPLICIT INTEGER*4 (A-E)
      INTEGER*4 NAME(8), LINE, BLANK
      INTEGER*4 ODEV1, ODV
      LOGICAL*4 SFLAG
      COMMON /IO/IDEV(4), ODEV1, ODV(3), LINE(80)
      DATA BLANK/' '/
C
C     INITIALIZATION
      ISIGN=1
      NUM=0
      FNUM=0.
      LWP1=1
      LWP2=80
      BC=BLANK
C
      DO 1 I=1,8
 1    NAME(I)=BLANK
C
C     MAKE SURE LPTR IS REASONABLE
      IF(LPTR.LE.0)LPTR=1
C     CHECK FOR END-OF-LINE CONDITION
 2    IF(LPTR.GT.80)RETURN 3
C
C     SEARCH FOR START OF FIRST INTERESTING THING
      C=LINE(LPTR)
      LPTR=LPTR+1
      CALL CINT(C,ICT,IDV)
      IF(SFLAG .AND. ICT.EQ.7)ICT=1
      GO TO (3,4,2,6,6,6),ICT
C
C     HAVE A BREAK CHARACTER OF SOME SORT
 7    BC=C
      RETURN 2
C
C     HAVE '+', '-', OR '.' WHICH MAY BE BREAK CHARACTERS OR
C     START OF NUMBER
C     IF A NUMBER IS COMING IN, NEXT CHAR SHOULD BE DIGIT IF 1ST CHAR
C     WAS '.', DIGIT OR . IF 1ST CHAR WAS + OR -
 6    IF(LPTR.GT.80)GO TO 7
      CALL CINT(LINE(LPTR), ICT1,IDV)
      IF(SFLAG .AND. ICT1.EQ.7)ICT1=1
      GO TO (7,8,7,9,7,7,7,7),ICT1
C
C     2ND CHAR WAS '.', WHICH IS OK IF 1ST CHAR WAS + OR -
 9    IF(ICT.EQ.4)GO TO 7
C     OK, HAVE + OR - FOLLOWED BY . SO FAR
C     DIGIT MUST BE NEXT, OR ELSE
      I=LPTR+1
      IF(I.GT.80)GO TO 7
      CALL CINT(LINE(I),ICT1,IDV)
      IF(ICT1.NE.2)GO TO 7
C     FINALLY, IT IS CLEAR THAT A NUMBER IS COMING
C
C     SET UP TO START CONVERTING NUMBER
 8    IF(ICT.EQ.6)ISIGN=-1
      IF(ICT.EQ.4)GO TO 10
      GO TO 20
C
C
C     A NUMBER IS BEING CONVERTED.  NO DECIMAL POINT ENCOUNTERED YET
 4    FNUM=10.*FNUM
      X=IDV
      FNUM=FNUM+X
C
 20   IF(LPTR.GT.80)GO TO 11
      C=LINE(LPTR)
      CALL CINT(C,ICT,IDV)
C     CHECK FOR DECIMAL POINT
      IF(ICT.EQ.4)GO TO 12
      IF(ICT.NE.2)GO TO 11
C     HAVE ANOTHER DIGIT
      LPTR=LPTR+1
      GO TO 4
C
C     HAVE HIT DECIMAL POINT IN NUMBER
 12   LPTR=LPTR+1
 10   XMULT=.1
C
 13   IF(LPTR.GT.80)GO TO 11
      C=LINE(LPTR)
      CALL CINT(C,ICT,IDV)
      IF(ICT.NE.2)GO TO 11
C     ANOTHER DIGIT
      X=IDV
      FNUM=FNUM+XMULT*X
      XMULT=.1*XMULT
      LPTR=LPTR+1
      GO TO 13
C
C     END OF NUMBER
 11   IF(ISIGN.LT.0)FNUM=-FNUM
      X=FNUM
      IF(X.LE.0.)X=X-1.
      X=X+.5
      NUM=X
C     RETURN WITH CONVERTED NUMBER
      RETURN 1
C
C
C     LOOKS LIKE AN ALPHAMERIC NAME COMING UP
 3    NAME(1)=C
      LWP1=LPTR-1
      NPTR=2
C
 15   IF(LPTR.GT.80)RETURN
      C=LINE(LPTR)
      LPTR=LPTR+1
      CALL CINT(C,ICT,IDV)
      IF(SFLAG .AND. ICT.EQ.7)ICT=1
      GO TO (14,14),ICT
C
C     END OF ALPHAMERIC NAME
      LPTR=LPTR-1
      RETURN
C
C     ADD CHARACTER TO NAME UNLESS TOO LONG
 14   IF(NPTR.GT.8)GO TO 16
      NAME(NPTR)=C
      NPTR=NPTR+1
      GO TO 15
C
C     LONG WORD--MORE THAN 8 CHARACTERS
C     FIND END OF LONG WORD
 16   LWP2=LPTR-1
      IF(LPTR.GT.80)GO TO 17
      C=LINE(LPTR)
      CALL CINT(C,ICT,IDV)
      IF(SFLAG .AND. ICT.EQ.7)ICT=1
      IF(ICT.GT.2)GO TO 17
      LPTR=LPTR+1
      GO TO 16
C
C     PRINT MESSAGE THAT WORD HAS BEEN SHORTENED
 17   WRITE(ODEV1,1001)(LINE(I),I=LWP1,LWP2)
      WRITE(ODEV1,1002)NAME
      RETURN
C
 1001 FORMAT(' UNDULY LONG WORD: ',80A1)
 1002 FORMAT(' SHORTENED TO: ',8A1)
      END
      SUBROUTINE CINT(C,ICT,IDV)
C
C     CHARACTER DECODING SUBROUTINE
C
C     ARGUMENT C MUST BE A CHARACTER IN THE LEFTMOST BYTE OF AN
C     INTEGER*4 WORD, WITH A BLANK IN THE RIGHT BYTE
C
C     ARGUMENT ICT IS AN INTEGER*4 WORD WHICH CONTAINS A NUMBER
C     INDICATING THE TYPE OF CHARACTER WHEN CINT RETURNS (SEE
C     BELOW FOR RETURN CODES)
C
C     IDV IS AN INTEGER*4 WORD USED TO RETURN THE NUMERIC VALUE
C     WHEN C IS AN INTEGER
C
C     RETURN CODES:
C	1  C IS ALPHABETIC
C	2  C IS DIGIT; NUMERIC VALUE OF DIGIT IS IN IDV
C	3  C IS A BLANK
C	4  C IS A PERIOD
C	5  C IS A +
C	6  C IS A -
C	7  C IS A *
C	8  C IS NONE OF THE ABOVE
C
C
      IMPLICIT INTEGER*4 (C)
      INTEGER*4 SPEC(5),ALPHA(26),NUM(10)
      DATA ALPHA/'A','B','C','D','E','F','G','H','I','J','K','L',
     1     'M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
      DATA NUM/'0','1','2','3','4','5','6','7','8','9'/
      DATA SPEC/' ','.','+','-','*'/
C
C
      IDV=0
C     SEE IF C IS ALPHABETIC
      ICT=1
      DO 10 I=1,26
  10  IF (C.EQ.ALPHA(I)) RETURN
C     SEE IF C IS NUMERIC
      DO 20 I=1,10
   20 IF (C.EQ.NUM(I)) GO TO 21
      GO TO 22
   21 ICT = 2
      IDV = I-1
      RETURN
   22  CONTINUE
C
C     CHECK FOR SPECIAL CHARACTER
 2    DO 3 I=1,5
      IF(C.NE.SPEC(I))GO TO 3
      ICT=I+2
      RETURN
 3    CONTINUE
C     CHARACTER C IS NONE OF THE ABOVE
      ICT=8
      RETURN
      END
      SUBROUTINE SFIND(S1,S2,N1,N2,START,SPTR,EPTR,*)
C
C     SUBROUTINE TO FIND THE FIRST OCCURRENCE, IF ANY, OF A
C     GIVEN STRING OF CHARACTERS IN A SECOND STRING, STARTING
C     AT A GIVEN POINT IN THE SECOND STRING.  RETURNS POINTERS
C     TO THE BEGINNING AND END POSITIONS OF THE MATCHED STRING
C     IN THE SECOND STRING.  NORMAL RETURN INDICATES THAT A
C     MATCH HAS OCCURRED, RETURN 1 INDICATES NO MATCH.
C
C     ARGUMENTS:
C	S1     INTEGER*4  VECTOR OF LENGTH N1 CONTAINING STRING
C			  TO BE MATCHED
C	S2     INTEGER*4  VECTOR OF LENGTH N2 CONTAINING STRING
C			  TO BE SEARCHED
C	N1     INTEGER*4  LENGTH OF S1
C	N2     INTEGER*4  LENGTH OF S2
C	START  INTEGER*4  POINTER TO 1ST CHARACTER TO BE EXAMINED
C			  IN S2
C	SPTR   INTEGER*4  POINTS TO START OF MATCH IN S2 ON
C			  NORMAL RETURN
C	EPTR   INTEGER*4  POINTS TO END OF MATCH IN S2 ON NORMAL
C			  RETURN
C
C
      IMPLICIT INTEGER*4 (A-H,O-Z)
      INTEGER*4 START,SPTR,EPTR,ODEV1,ODV
      DIMENSION S1(N1), S2(N2)
      COMMON /IO/IDEV(4),ODEV1,ODV(3),ALINE(80)
C
      SPTR=0
      EPTR=0
C     CHECK STRING LENGTHS FOR LEGALITY
      IF(N1.GT.0 .AND. N2.GT.0)GO TO 1
C     BAD VALUE FOR N1 AND/OR N2
      WRITE(ODEV1,1001)N1,N2,START
 1001 FORMAT(' BAD CALL TO SFIND  N1=',I8,'  N2=',I8,
     1	'  START=',I8)
      STOP 767
C
 1    IF(START.LE.0)START=1
      I=START
C
C     CHECK FOR END OF 2ND STRING
 2    IF(I.GT.N2)RETURN 1
C
C     TRY TO MATCH 1ST CHAR OF S1 WITH A CHAR IN S2
      IF(S1(1).EQ.S2(I))GO TO 3
C
C     NO LUCK
 4    I=I+1
      GO TO 2
C
C     SEE IF CAN MATCH REST OF S1 WITH SUCCEEDING CHARS IN S2
 3    IF(N1.EQ.1)GO TO 5
      DO 6 J=2,N1
      K=I+J-1
      IF(K.GT.N2)RETURN 1
      IF(S2(K).NE.S1(J))GO TO 4
 6    CONTINUE
 5    SPTR=I
      EPTR=I+N1-1
      RETURN
      END
      SUBROUTINE INPUT(IDEVNO,ECHO,*,*)
C
C     SUBROUTINE TO READ A LINE OF TEXTUAL DATA FROM A GIVEN DEVICE
C     AND STORE IT IN LINE.
C
C     NORMAL RETURN INDICATES THAT THE LINE WAS SUCCESSFULLY
C     READ, AND NO SUPERVISOR COMMAND PREFIX CHARACTER PAIRS
C     ('>>', '<<', OR '&&') WERE FOUND.
C     RETURN 1 INDICATES AN END-OF-FILE WAS ENCOUNTERED
C     RETURN 2 INDICATES THE LINE CONTAINS COMMAND PREFIX CHARS
C
C     ARGUMENT:
C	IDEVNO	INTEGER*4  NUMBER OF DEVICE FROM WHICH INPUT
C			   IS TO BE READ.  LEGAL VALUES RANGE
C			   FROM 0 TO 99.
C	ECHO	LOGICAL*4  IF .TRUE., CAUSES INPUT LINES TO BE PRINTED
C			   ON ODEV1.
C
C
      IMPLICIT INTEGER*4 (L), INTEGER*4 (A-K,M-Z)
      COMMON /IO/IDEV(4),ODEV(4),LINE(80)
C
C     CPFX CONTAINS THE ALLOWABLE COMMAND PREFIX CHARACTER PAIRS.
C     NCPFX IS THE NUMBER OF SUCH PAIRS.  TO CHANGE PREFIX CHAR-
C     ACTERS, ALTER THE DATA BELOW AND ALSO IN THE SUPERVISOR MAIN
C     PROGRAM.
      INTEGER*4 CPFX(2,3)
!!      DATA NCPFX/3/, CPFX/'>','>','<','<','&','&'/
C
      LOGICAL*4 ECHO
      EQUIVALENCE (ODEV1,ODEV(1))
      DATA NCPFX/3/, CPFX/'>','>','<','<','&','&'/
C
C     TEST LEGALITY OF CALL
      IF(IDEVNO.GE.0 .AND. IDEVNO.LT.100)GO TO 1
C     BAD CALL
      WRITE(ODEV1,1001)IDEVNO
      STOP '6769'
C
C     TRY TO READ LINE
 1    READ(IDEVNO,1002,END=2,ERR=3)LINE
C
      IF(ECHO) JMAX=JMAXX(LINE,80)
      IF(ECHO) WRITE(ODEV1,1005) (LINE(III),III=1,JMAX)
C
C     SEARCH LINE FOR SUPERVISOR COMMAND
      DO 4 I=1,NCPFX
      CALL SFIND(CPFX(1,I),LINE,2,80,1,J,K,&4)
C     THERE IS A SUPERVISOR COMMAND IN THE LINE
      RETURN 2
C
 4    CONTINUE
C     NO COMMAND--NORMAL RETURN
      RETURN
C
C     END-OF-FILE ENCOUNTERED.	TELL THE WORLD.
 2    WRITE(ODEV1,1003)IDEVNO
      RETURN 1
C
C     DEVICE ERROR--TIME TO QUIT BEFORE THINGS GET WORSE
 3    WRITE(ODEV1,1004)IDEVNO
      STOP 6770
C
 1001 FORMAT(' BAD CALL TO INPUT:  DEVICE NO.=',I8)
 1002 FORMAT(80A1)
 1003 FORMAT(1H0,'END-OF-FILE ON DEVICE',I3)
 1004 FORMAT(1H0,'DEVICE',I3,' HAS FAILED.')
 1005 FORMAT(1X,80A1)
      END
      SUBROUTINE NMATCH(NAME1,NAME2,*)
C
C     SUBROUTINE TO DETERMINE WHETHER OR NOT TWO 8-CHARACTER
C     ALPHAMERIC NAMES MATCH
C
C     BOTH NAME1 AND NAME2 MUST CONTAIN ONLY ALPHAMERIC CHAR-
C     ACTERS AND TRAILING BLANKS, EXCEPT THAT NAME1 MAY CONTAIN
C     ASTERISKS TO DENOTE VARIABLE STRING ELEMENTS.  NAME1 MAY
C     CONTAIN UP TO 4 *'S; IT IS NOT LEGAL TO PUT TWO ASTERISKS
C     IN A ROW.
C
C     NMATCH RETURNS NORMALLY IF THE TWO NAMES MATCH, DOES A
C     RETURN 1 IF THEY DO NOT.
C
C     ARGUMENTS:
C	NAME1  INTEGER*4  8 CHARACTER VECTOR, MAY CONTAIN *'S
C	NAME2  INTEGER*4  8 CHARACTER VECTOR, MAY NOT CONTAIN *'S
C
C
      IMPLICIT INTEGER*4 (A-H,O-Z)
      INTEGER*4 NAME1(8),NAME2(8)
C
      DATA STAR/'*'/, BLANK/' '/
C
C     TRY TO MATCH 1ST CHUNK OF NAME1 UP TO 1ST * OR END
      N1PTR=0
      N2PTR=0
      DO 1 I=1,8
      IF(NAME1(I).NE.STAR)GO TO 1
      IF(I.EQ.1)GO TO 3
      I2=I-1
C
C     MUST HAVE EXACT MATCH FROM 1 TO I2
 2    DO 4 J=1,I2
      IF(NAME1(J).NE.NAME2(J))RETURN 1
 4    CONTINUE
C     SO FAR SO GOOD
      N1PTR=I2
      N2PTR=I2
      GO TO 3
 1    CONTINUE
C     NO STARS MEANS EXACT MATCH NECESSARY
      I2=8
      GO TO 2
C
C
C     SEE IF DONE
 3    IF(N1PTR.LT.8)GO TO 5
      IF(N2PTR.EQ.8 .OR. NAME2(N2PTR+1).EQ.BLANK)RETURN
      RETURN 1
C
 5    IF(NAME1(N1PTR+1).EQ.STAR)GO TO 11
      IF(N2PTR.EQ.8)RETURN
      IF(NAME2(N2PTR+1).EQ.BLANK)RETURN
      RETURN 1
C
C     FIND NEXT CLUMP OF NON-BLANK CHARS FOLLOWING *, IF ANY
 11   I1=N1PTR+2
C     IF LAST SIGNIFICANT CHARACTER OF NAME1 IS *, CAN QUIT
      IF(I1.GT.8)RETURN
      IF(NAME1(I1).EQ.BLANK)RETURN
C
C     MORE ALPHAMERIC CHARACTERS TO BE MATCHED
      DO 6 I=I1,8
      IF(NAME1(I).NE.STAR .AND. NAME1(I).NE.BLANK)GO TO 6
C     HAVE END OF CLUMP OF CHARACTERS
      I2=I-1
      GO TO 7
 6    CONTINUE
      I2=8
C
C     SCAN THRU NAME2 FOR SEQUENCE OF CHARS MATCHING CLUMP IN NAME1
 7    N2=N2PTR+1
      IF(N2.GT.8)RETURN 1
C
C     TRY TO MATCH 1ST CHAR OR CLUMP
      N1I1=NAME1(I1)
      DO 8 I=N2,8
      IF(NAME2(I).NE.N1I1)GO TO 8
C     HAVE 1ST CHAR MATCH--TRY FOR REST
      IF(I2.EQ.I1)GO TO 9
C
      I11=I1+1
      DO 10 J=I11,I2
      K=I+J-I1
      IF(K.GT.8)RETURN 1
      IF(NAME2(K).NE.NAME1(J))GO TO 8
 10   CONTINUE
C
C     IF CLUMP IN NAME1 IS AT END OF WORD, NAME2 CLUMP SHOULD BE THERE TOO
 9    IF(I2.LT.8 .AND. NAME1(I2+1).NE.BLANK)GO TO 12
C     NOW AT END OF WORD IN NAME1--BETTER BE AT END IN NAME2
      K=I+I2-I1
      IF(K.LT.8 .AND. NAME2(K+1).NE.BLANK)GO TO 8
C
C     SUCCESS
 12   N1PTR=I2
      N2PTR=I+I2-I1
      GO TO 3
 8    CONTINUE
      RETURN 1
      END
      SUBROUTINE NSRCH(NAME,NTABLE,START,STOP,NNUM,*)
C
C     SUBROUTINE TO DETERMINE WHETHER A GIVEN NAME IS IN A GIVEN
C     RANGE OF A TABLE OF NAMES, AND, IF SO, WHICH ONE IT IS
C
C     A NORMAL RETURN IS EXECUTED WHEN THE NAME IS FOUND IN THE
C     TABLE, RETURN 1 IF IT IS NOT
C
C     ARGUMENTS:
C	NAME   INTEGER*4  8 CHARACTER VECTOR CONTAINING NAME TO
C			  BE LOOKED FOR IN NTABLE.  MAY NOT CON-
C			  TAIN *'S.
C	NTABLE INTEGER*4  TABLE OF NAMES AND ABBREVIATIONS DIMEN-
C			  SIONED (8,N), WHERE N IS GREATER THAN OR
C			  OR EQUAL TO STOP.  NAMES IN NTABLE MAY
C			  CONTAIN *'S.
C	START  INTEGER*4  POINTER TO THE 1ST NAME IN NTABLE TO BE
C			  SCANNED FOR A MATCH WITH NAME.
C	STOP   INTEGER*4  POINTER TO THE LAST NAME TO BE SCANNED.
C	NNUM   INTEGER*4  CONTAINS POINTER TO THE FIRST NAME WHICH
C			  MATCHED NAME ON NORMAL RETURN.  EQUALS
C			  ZERO ON RETURN 1.
C
C
      IMPLICIT INTEGER*4 (A-Z)
      INTEGER*4 NAME(8),NTABLE(8,STOP),LINE,NT(8)
      COMMON /IO/IDEV(4),ODEV(4),LINE(80)
      EQUIVALENCE (ODEV1,ODEV(1))
C
C
C     TEST LEGALITY OF ARGS
      IF(STOP.GE.1 .AND. START.LE.STOP)GO TO 1
C     ILLEGAL SITUATION
      WRITE(ODEV1,1001)START,STOP
      STOP '6768'
C
 1    IF(START.LE.1)START=1
      DO 2 I=START,STOP
      DO 3 J=1,8
 3    NT(J)=NTABLE(J,I)
      CALL NMATCH(NT,NAME,&2)
C
C     HAVE MATCH
      NNUM=I
      RETURN
 2    CONTINUE
C     NO LUCK
      NNUM=0
      RETURN 1
C
 1001 FORMAT(' BAD CALL TO NSRCH  START=',I8,'	STOP=',I8)
      END
      SUBROUTINE ATQ(BC,LPTR,*)
C
C     SUBROUTINE TO TEST FOR PRESENCE OF '@END' IN LINE(80)
C     ENTRY ATQ TESTS BC TO SEE IF IT IS '@', AND, IF SO,
C     TESTS FOR E* AFTER THE '@'.
C     ENTRY ATEND ASSUMES THAT AN '@' HAS ALREADY BEEN DETECTED
C     AND JUST CHECKS FOR E*.
C
C     BOTH ENTRIES GIVE NORMAL RETURNS IF '@END' IS NOT FOUND,
C     RETURN 1 IF IT IS.
C     IN NO CASE IS LPTR CHANGED.
C
C     ARGUMENTS:
C	LPTR  INTEGER*4  POINTER INTO LINE(80)
C	BC    INTEGER*4  BREAK CHAR TO BE TESTED FOR '@'
C
C     THIS VERSION OF ATQ AND ATEND ACCEPTS '@' ALONE AS
C     END-OF-STRING WHETHER OR NOT FOLLOWED BY E*.
C
C
      IMPLICIT INTEGER*4 (O), INTEGER*4 (A-H)
      INTEGER*4 N(8),ENDX(8),LINE
      COMMON /IO/IDEV(4),ODEV1,ODEV(3),LINE(80)
!!      DATA AT/'@'/
!!      DATA ENDX/'E','*',6*' '/
      LOGICAL*4 F
      DATA AT/'@'/
      DATA ENDX/'E','*',6*' '/
      DATA F/.FALSE./
C
C
C     SEE IF BC IS AN '@'
      IF(BC.EQ.AT)GO TO 3
      RETURN
C
C
      ENTRY ATEND(LPTR,*)
C
 3    I=LPTR
      CALL NEXT(I,N,BC,J,X,F,&1,&1,&1)
C     SEE IF NEXT THING AFTER '@' IS E*
      CALL NMATCH(ENDX,N,&1)
C
C     OK--HAVE CLEAR EOS
      RETURN 1
C
C     '@' NOT FOLLOWED BY E*
 1    CONTINUE
      RETURN 1
C
      END
      SUBROUTINE CVALS(NC)
C
C     SUBROUTINE TO GET THE VALUES OF ALL SIMULATION VARIABLES
C     INTO THE 'CURRENT VALUES' VECTORS IN /VARS1/
C
C     ARGUMENT:
C     NC  I*4  NUMBER OF EXPTL CONDITION
C
C------------------------------
C
      INTEGER*4 SPECI(12,32),SPECF(12,32),SI(12),SF(12)
      INTEGER*4 IVARS(12,32),IV(12)
      REAL*4 FVARS(12,32),FV(12)
      LOGICAL*4 FLAGS(12,32),FLGS(12),ENTERD(36)
      COMMON /VARS/ IVARS,FVARS,SPECI,SPECF,FLAGS,ENTERD
      COMMON /VARS1/ IV,FV,SI,SF,FLGS
C
      IF(NC.GT.0 .AND. NC.LE.32)GO TO 1
      STOP '1973'
 1    DO 2 I=1,12
      FLGS(I)=FLAGS(I,NC)
      IV(I)=IVARS(I,NC)
      FV(I)=FVARS(I,NC)
      SI(I)=SPECI(I,NC)
 2    SF(I)=SPECF(I,NC)
      RETURN
      END
      SUBROUTINE EVAL(VNO,TYPE,VSPTR,NC,X,*)
C
C     SUBROUTINE USED IN MAKING LEGALITY TESTS ON INPUT
C     GETS FP VALUE FOR SPECIFIED VARIABLE; IF VARIABLE NUMBER IS
C     NEGATIVE, THE 'SPECIAL' VALUE OF THE VARIABLE IS WHAT IS RETURNED
C     RETURN 1 INDICATES THAT THE VARIABLE HAS A SPECIAL VALUE WHEN A
C     NORMAL VALUE IS CALLED FOR, OR VICE VERSA.  IN THAT CASE, A ZERO
C     VALUE IS RETURNED IN X.
C
C     ARGUMENTS:
C     VNO    I*2  VARIABLE NUMBER (SEE ABOVE)
C     TYPE   I*2  INTERNAL TYPE OF VARIABLE
C     VSPTR  I*2  VARIABLE STORAGE POINTER
C     NC     I*4  EXPTL CONDITION NUMBER
C     X      R*4  CONTAINS VALUE OF VARIABLE ON RETURN
C
C
      INTEGER*4 VNO,TYPE,VSPTR
C
C     GET FP VALUE FOR VARIABLE
      CALL VALUE(VSPTR,TYPE,NC,X,IX,&1)
      IF((VNO.GE.0).AND.(IX.EQ.0))RETURN
      IF((VNO.GE.0).AND.(IX.NE.0))GO TO 2
      IF((VNO.LT.0).AND.(IX.EQ.0))GO TO 2
      X=IX
      RETURN
 2    X=0.
      RETURN 1
C
C     ERROR
 1    STOP '1974'
      END
      SUBROUTINE ABBREV(LIST,NPTR,LAST,IDEV,LPTR,ECHO,*)
C
C     SUBROUTINE TO READ A LIST OF NAMES/ABBREVIATIONS ENCLOSED IN
C     PARENTHESES
C
C     NAMES ARE STORED SEQUENTIALLY IN LIST STARTING AT SLOT NPTR+1
C     LIST OF NAMES MUST BE TERMINATED BY ')'; END-OF-FILE OR @END
C     WILL CAUSE ERROR MESSAGE AND ERROR RETURN
C     IF TOO MANY ERRORS ARE ENCOUNTERED, ABBREV ASSUMES THAT THE
C     RIGHT PARENTHESIS IS MISSING AND QUITS PROCESSING
C     INITIAL LEFT PARENTHESIS IS NOT NECESSARY
C
C     RETURN 1 INDICATES FATAL ERROR(S)
C
C     ARGUMENTS:
C     LIST  I*2  LIST IN WHICH NAMES/ABBREVNS ARE TO BE STORED.
C		 DIM(8,LAST)
C     NPTR  I*4  POINTER TO THE LAST SLOT USED IN LIST.  ON RETURN,
C		 POINTS TO THE LAST SLOT USED BY ABBREV
C     LAST  I*4  MAXIMUM NUMBER OF NAMES THAT CAN BE STORED IN LIST
C     IDEV  I*4  INPUT DEVICE NUMBER
C     LPTR  I*4  POINTER INTO LINE TO THE 1ST CHARACTER TO BE SCANNED
C		 BY ABBREV.  ON RETURN, POINTS TO THE CHARACTER AFTER
C		 THE LAST ONE PROCESSED BY ABBREV
C     ECHO  L*1  WHEN ECHO=.TRUE. INPUT IS ECHOED ON ODEV1
C
C------------------------------
C
      IMPLICIT INTEGER*4 (A-H,P-Z), INTEGER*4 (O)
      INTEGER*4 LIST(8,LAST), ABV(8), LINE
      LOGICAL*4 T,ECHO
      REAL*4 FNUM
!!      DATA T/.TRUE./, RPAREN/') '/
      COMMON /IO/IDV(4),ODEV1,ODV(3),LINE(80)
      DATA T/.TRUE./, RPAREN/') '/
C
C------------------------------
C
C     TEST FOR PROGRAM ERROR
      IF(NPTR.LT.0)GO TO 112
      IF(LAST.LE.0)GO TO 112
      IF(IDEV.LT.0 .OR. IDEV.GT.99)GO TO 112
      NERR=0
C
C------------------------------
C
C     GET NEXT THING FROM LIST OF ABBREVIATIONS
 1    CALL NEXT(LPTR,ABV,BC,NUM,FNUM,T,&100,&10,&11)
C
C     HAVE ANOTHER ABBREV--ADD TO LIST IF HAVE ROOM
      IF(NPTR.LT.LAST)GO TO 2
C
C     LIST OVERFLOW
      WRITE(ODEV1,1001)ABV
      GO TO 110
C
 2    NPTR=NPTR+1
      DO 3 I=1,8
 3    LIST(I,NPTR)=ABV(I)
      GO TO 1
C
C     BREAK CHAR MAY BE ')' OR '@'; IGNORE IF NOT
 10   CALL ATQ(BC,LPTR,&111)
      IF(BC.NE.RPAREN)GO TO 1
C     END OF ABBREVN LIST
      IF(NERR.GT.0)RETURN 1
      RETURN
C
C     END OF LINE
 11   CALL INPUT(IDEV,ECHO,&102,&120)
 120  LPTR=1
      GO TO 1
C
C------------------------------
C
C     NUMBER IN INPUT
 100  WRITE(ODEV1,1002)FNUM
      GO TO 110
C
C     END-OF-FILE
 102  WRITE(ODEV1,1003)
      GO TO 111
C
C     IF TOO MANY ERRORS, PROBABLY HAVE UNCLOSED PARENTHESIS
 110  NERR=NERR+1
      IF(NERR.LT.5)GO TO 1
C
C     GIVE UP
 111  WRITE(ODEV1,1004)LINE
      IF(NPTR.LE.0)RETURN 1
      WRITE(ODEV1,1006)((LIST(I,IS1),I=1,8),IS1=1,NPTR)
      RETURN 1
C
C     BAD CALL
 112  WRITE(ODEV1,1007)NPTR,LAST,IDEV,LPTR
      STOP '1921'
C
C------------------------------
C
 1001 FORMAT(1H0,'NAME LIST OVERFLOW.  NAME:',8A1)
 1002 FORMAT(1H0,'NUMBER IN LIST OF NAMES:',F12.4)
 1003 FORMAT(1H0,'ABBREV ABORT')
 1004 FORMAT(1H0,'MISSING RT PARENTHESIS.  LAST LINE:'/1X,80A1)
 1006 FORMAT(' NAME LIST:'/(1X,7(8A1,', ')))
 1007 FORMAT(1H0,'BAD CALL TO ABBREV'/'   NPTR	LAST  IDEV'
     1 ,'  LPTR'/1X,4I6)
      END
      SUBROUTINE CCHARS(NC,CC)
C
C     SUBROUTINES DEALING WITH ALPHABETIC CONDITION LABELS
C     CCHARS PRODUCES A LABEL GIVEN A CONDITION NUMBER
C     CCNUM PRODUCES A CONDITION NUMBER GIVEN A LABEL
C
C     CCHARS
C     WORKS ONLY FOR CONDITION NUMBERS BETWEEN 1 AND 26**2
C
C     ARGUMENTS:
C     NC  I*4  NUMBER OF CONDITION
C     CC  I*2  DIM(2).	RETURN VECTOR FOR PAIR OF CONDITION CHARACTERS
C	       CHARACTERS ARE RIGHT JUSTIFIED IN CC.
C
C------------------------------
C
      INTEGER*4 CC(2), ABC(26), BLANK
      INTEGER*4 NAME(8)
      INTEGER*4 ODEV1
      COMMON /IO/IDEV(4),ODEV1
        DATA BLANK/'  '/
        DATA ABC/'A','B','C','D','E','F','G','H','I','J','K','L',
     1    'M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
C
      I=0
      NC1=NC
      CC(1)=BLANK
      IF(NC1.LE.26)GO TO 1
 2    I=I+1
      NC1=NC1-26
      IF(NC1.GT.26)GO TO 2
      CC(1)=ABC(I)
 1    CC(2)=ABC(NC1)
      RETURN
C
C------------------------------
C
C     CCNUM
C     PRINTS AN ERROR MESSAGE AND DOES RETURN 1 IF CONDITION LABEL IS
C     TOO LONG OR SPECIFIES A CONDITION NOT DEFINED BY THE STUDENT
C
C     ARGUMENTS:
C     NAME    I*2  DIM(8).  CHARACTER VECTOR CONTAINING CONDITION LABEL
C		   LEFT JUSTIFIED WITH AT LEAST ONE TRAILING BLANK
C     NUM     I*4  CONDITION NUMBER IS RETURNED IN NUM
C     NCONDS  I*4  MUST CONTAIN THE NUMBER OF CONDITIONS DEFINED BY THE
C		   STUDENT
C
C------------------------------
C
      ENTRY CCNUM(NAME,NUM,NCONDS,*)
!!      INTEGER*4 NAME(8)
!!      INTEGER*4 ODEV1
!!      COMMON /IO/IDEV(4),ODEV1
C
      NUM=0
      DO 10 I=1,3
      IF(NAME(I).EQ.BLANK)GO TO 12
      NUM=26*NUM
      DO 11 J=1,26
      IF(NAME(I).EQ.ABC(J))GO TO 13
 11   CONTINUE
C     UNIDENTIFIABLE CHARACTER
      GO TO 20
C
 13   NUM=NUM+J
 10   CONTINUE
C     TOO MANY CHARACTERS
      GO TO 20
C
C     SEE IF NUMBER IS TOO BIG
 12   IF(NUM.LE.0 .OR. NUM.GT.NCONDS)GO TO 20
      RETURN
C
C     ILLEGAL CONDITION LABEL
 20   WRITE(ODEV1,1001)NAME
      NUM=1
      RETURN 1
C
 1001 FORMAT(1H0,'ILLEGAL CONDITION LABEL: ',8A1)
      END
      SUBROUTINE NXT(LPTR,NAME,NUM,FNUM,ECHO,IDEV,*,*,*,*,*,*)
C
C     AUXILIARY SUBROUTINE FOR USE BY SIM DATA LOADER
C     ADAPTS SUBROUTINE NEXT FOR SPECIAL APPLICATION
C
C     RETURNS:
C     0  NAME
C     1  NUMBER
C     2  SLASH
C     3  @END
C     4  (
C     5  )
C     6  EOF
C
C------------------------------
C
      INTEGER*4 NAME(8),BC,SLASH,RPAREN,LPAREN
      LOGICAL*4 F,ECHO
      DATA F/.FALSE./, SLASH/'/'/, LPAREN/'('/, RPAREN/')'/
C
C------------------------------
C
 6    CALL NEXT(LPTR,NAME,BC,NUM,FNUM,F,&1,&2,&3)
      RETURN
 1    RETURN 1
 2    CALL ATQ(BC,LPTR,&5)
      IF(BC.EQ.SLASH)RETURN 2
      IF(BC.EQ.LPAREN)RETURN 4
      IF(BC.EQ.RPAREN)RETURN 5
      GO TO 6
C
 3    CALL INPUT(IDEV,ECHO,&7,&8)
 8    LPTR=1
      GO TO 6
C
 5    RETURN 3
 7    RETURN 6
      END
      SUBROUTINE VALUE(PTR,TYPE,NC,X,IX,*)
C
C     SUBROUTINE TO GET VALUE OF VARIABLE
C
      INTEGER*4 PTR,TYPE
      INTEGER*4 SPECI(12,32),SPECF(12,32)
      LOGICAL*4 FLAGS(12,32),ENTERD(36)
      COMMON /VARS/IVARS(12,32),FVARS(12,32),SPECI,SPECF,FLAGS,ENTERD
C
C
      I=TYPE
      GO TO (1,2,1,1,2,1,9,1,2,9),I
C     INTEGER VARIABLE
 1    X=IVARS(PTR,NC)
      IX=SPECI(PTR,NC)
      RETURN
C     FLOATING POINT VARIABLE
 2    X=FVARS(PTR,NC)
      IX=SPECF(PTR,NC)
      RETURN
C     FLAG VARIABLE
 9    RETURN 1
      END
      SUBROUTINE RJUST(NVEC,NVP,VLIST,VPTR)
C
C     SUBROUTINE TO TAKE A NAME FROM A LIST OF LEFT-JUSTIFIED NAMES,
C     RIGHT-JUSTIFY IT, AND PLACE IT IN A LIST OF RIGHT-JUSTIFIED NAMES
C
C     ARGUMENTS:
C     NVEC   I*2  DIM(8,NVP+).	VECTOR OF NAMES INTO WHICH RIGHT-
C		  JUSTIFIED NAME IS TO BE PLACED
C     NVP    I*4  POINTER INTO NVEC INDICATING WHERE RIGHT-JUSTIFIED
C		  NAME SHOULD BE PUT
C     VLIST  I*2  DIM(8,75 OR SO).  VECTOR OF LEFT-JUSTIFIED NAMES
C		  FROM WHICH NAME IS TO BE TAKEN
C     VPTR   I*2  POINTER INDICATING WHICH NAME IS TO BE TAKEN FROM
C		  VLIST
C
C------------------------------
C
C
      INTEGER*4 NVEC(8,NVP),VLIST(8,75),VPTR,BLANK
      DATA BLANK/' '/
C
      DO 1 I=1,8
 1    NVEC(I,NVP)=VLIST(I,VPTR)
C
      DO 2 I=1,8
      IF(NVEC(8,NVP).NE.BLANK)GO TO 3
      DO 4 J=1,7
      K=8-J
 4    NVEC(K+1,NVP)=NVEC(K,NVP)
      NVEC(1,NVP)=BLANK
 2    CONTINUE
 3    RETURN
      END
      SUBROUTINE INTERP(IDEV,EOS,ECHOF,NVARS,NVN,NKW,NIVAR,NFVAR,
     1 NFLAGS,NERR,MAXC,NPC,NCG,VLIST,VNUM,GVNPTR,VSPTR,ETYPE,
     2 ENTERD,KWR,KWLIST,KVAL,IVARS,SPECI,DIVARS,DSPECI,FVARS,
     3 SPECF,DFVARS,DSPECF,FLAGS,DFLAGS,*,*)
C
C     INPUT INTERPRETER SUBROUTINE
C
C     DETERMINES VARIABLE SETTINGS FROM FREE FORM INPUT
C
C     IF MULTIPLE VALUES ARE GIVEN FOR ONE OR MORE VARIABLES, A COM-
C     PLETE FACTORIAL DESIGN WILL BE GENERATED
C
C     ALL OUTPUT GOES ON ODEV1
C
C     RETURN 1 MEANS END-OF-FILE WAS ENCOUNTERED
C     RETURN 2 MEANS SUPERVISOR COMMAND PREFIX CHARACTERS WERE FOUND
C     IN THE CURRENT INPUT LINE.  ENTRY INTP1 SHOULD BE CALLED TO
C     RESUME INPUT ANALYSIS WITH THE NEXT LINE FROM IDEV.
C
C
C     ARGUMENTS:
C     IDEV    I*4  INPUT DEVICE NUMBER BETWEEN 0 AND 99
C     EOS     I*4  EOS DETERMINES WHAT CONDITIONS REPRESENT END-OF-
C		   STRING.  IF EOS IS NEGATIVE OR 0, "@END" AND EOF ARE
C		   THE ONLY END-OF-STRING INDICATORS; IF EOS=K>0, INPUT
C                  WILL ALSO BE TERMINATED AFTER K LINES EXCLUDING
C                  SUPERVISIOR COMMAND LINES.
C     ECHOF   L*1  ECHOF=.TRUE. CAUSES INPUT LINES TO BE ECHOED ON
C		   ODEV1
C     NVARS   I*4  TOTAL NUMBER OF VARIABLES OF ALL TYPES
C     NVN     I*4  NUMBER OF VARIABLE NAMES AND ABBREVIATIONS
C     NKW     I*4  NUMBER OF KEYWORDS AND KEYWORD ABBREVIATIONS
C     NIVAR   I*4  NUMBER OF INTEGER VARIABLES
C     NFVAR   I*4  NUMBER OF FLOATING POINT VARIABLES
C     NFLAGS  I*4  NUMBER OF FLAG VARIABLES
C     NERR    I*4  ON RETURN, CONTAINS NUMBER OF SERIOUS ERRORS
C		   DETECTED IN THE INPUT STRING
C     MAXC    I*4  MAXIMUM NUMBER OF CONDITIONS ALLOWED
C     NPC     I*4  NUMBER OF CONDITIONS GENERATED BY PRIOR STRINGS
C     NCG     I*4  CONTAINS NUMBER OF NEW CONDITIONS GENERATED BY NEW
C		   INPUT ON RETURN
C     VLIST   I*2  LIST OF VARIABLE NAMES AND ABBREVIATIONS.
C		   DIM(8,NVN)
C     VNUM    I*2  LIST OF VARIABLE NUMBERS CORRESPONDING TO THE NAMES
C		   AND ABBREVIATIONS IN VLIST.	DIM(NVN)
C     GVNPTR  I*2  LIST OF POINTERS INTO VLIST TO THE GENERIC NAME OF
C		   EACH VARIABLE.  DIM(NVARS)
C     VSPTR   I*2  LIST OF POINTERS INTO IVARS/FVARS (AND SPECI/SPECF)
C		   OR FLAGS WHICH POINT TO THE SLOT FOR EACH VARIABLE
C		   IN THE APPROPRIATE ARRAY.  DIM(NVARS)
C     ETYPE   I*2  LIST CONTAINING EXTERNAL VARIABLE TYPE FOR EACH VAR-
C		   IABLE.  DIM(NVARS).	TYPES ARE:
C		     1	NUMI
C		     2	NUMF
C		     3	KWD
C		     4	KI
C		     5	KF
C		     6	IND
C		     7	FLAG
C		     8	INTERNAL-USE-ONLY INTEGER
C		     9	INTERNAL FLOATING POINT
C		    10	INTERNAL FLAG
C     ENTERD  L*1  VECTOR WHICH CONTAINS ON RETURN THE VALUE .TRUE.
C		   FOR EVERY VARIABLE EXPLICITLY GIVEN A VALUE BY THE
C		   USER, .FALSE. FOR ALL OTHERS.  DIM(NVARS)
C     KWR     I*2  ARRAY OF POINTERS INTO KWLIST.  DIM(2,NVARS).
C		   KWR(1,VARNO) POINTS TO THE 1ST KWD VALUE FOR THE
C		   GIVEN VARIABLE, KWR(2,VARNO) TO THE LAST
C     KWLIST  I*2  LIST OF ALL KEYWORDS AND KWD ABBREVIATIONS.
C		   DIM(8,NKW)
C     KVAL    I*4  VECTOR CONTAINING THE NUMERIC VALUE ASSOCIATED WITH
C		   EACH KWD/ABBREV IN KWLIST.  DIM(NKW)
C     IVARS   I*4  TABLE CONTAINING VALUES OF INTEGER VARIABLES FOR
C		   EACH CONDITION ON RETURN.  DIM(NIVAR,MAXC)
C     SPECI   I*2  TABLE WHICH CONTAINS NONZERO VALUE ON RETURN IF A
C		   KEYWORD IS GIVEN AS A VALUE FOR A KI VARIBLE.
C		   DIM(NIVAR,MAXC)
C     DIVARS  I*4  VECTOR CONTAINING DEFAULT VALUES TO GO IN IVARS
C		   FOR EACH INTEGER VARIABLE.  DIM(NIVAR)
C     DSPECI  I*2  VECTOR CONTAINING DEFAULT VALUES TO GO IN SPECI.
C		   DIM(NIVAR)
C     FVARS   R*4  TABLE CONTAINING VALUES OF FLOATING POINT VAR-
C		   IABLES FOR EVERY CONDITION ON RETURN.
C		   DIM(NFVAR,MAXC)
C     SPECF   I*2  TABLE WHICH CONTAINS NONZERO VALUE ON RETURN
C		   IF A KEYWORD IS GIVEN AS A VALUE FOR A KF
C		   VARIABLE.  DIM(NFVAR,MAXC)
C     DFVARS  R*4  VECTOR CONTAINING DEFAULT VALUES TO GO IN FVARS.
C		   DIM(NFVAR)
C     DSPECF  I*2  VECTOR CONTAINING DEFAULT VALUES TO GO IN
C		   SPECF.  DIM(NFVAR)
C     FLAGS   L*1  TABLE IN WHICH VALUES OF FLAG VARIABLES ARE RET-
C		   URNED.  DIM(NFLAGS,MAXC)
C     DFLAGS  L*1  VECTOR OF DEFAULT VALUES FOR FLAGS.	DIM(NFLAGS)
C
C--------------------
C--------------------
C
      IMPLICIT INTEGER*4 (O), INTEGER*4 (A-H,P-Z)
      INTEGER*4 LINE
      COMMON /IO/IDV(4),ODEV1,ODV(3),LINE(80)
C
C     TYPE SPECS AND DIMENSIONS FOR ARGUMENTS
      INTEGER*4 VLIST(8,NVN), VNUM(NVN), GVNPTR(NVARS), VSPTR(NVARS),
     1 ETYPE(NVARS), KWR(2,NVARS)
      INTEGER*4 KWLIST(8,NKW), SPECI(NIVAR,MAXC), DSPECI(NIVAR),
     1 SPECF(NFVAR,MAXC), DSPECF(NFVAR)
C
      INTEGER*4 EOS, IVARS(NIVAR,MAXC), DIVARS(NIVAR), KVAL(NKW)
C
      REAL*4 FVARS(NFVAR,MAXC), DFVARS(NFVAR)
C
      LOGICAL*4 ECHOF, ENTERD(NVARS), FLAGS(NFLAGS,MAXC),
     1 DFLAGS(NFLAGS)
C
C     TYPE SPECS FOR LOCAL VARIABLES
      INTEGER*4 NAME(8), TRASH(8), VNAME(8)
      INTEGER*4 TYPE
      REAL*4 FNUM
      LOGICAL*4 SPF, EOF, F
C
      DATA EQUALS/'='/
      DATA F/.FALSE./
C
C--------------------
C--------------------
C
C     MAKE SURE NOTHING FUNNY IS GOING ON
      IF(IDEV.LT.0 .OR. IDEV.GT.99)GO TO 1
      IF(NVARS.LE.0)GO TO 1
      IF(NVN.LT.NVARS)GO TO 1
      IF(NKW.LT.0)GO TO 1
      IF(NIVAR.LT.0 .OR. NFVAR.LT.0)GO TO 1
      IF((NIVAR+NFLAGS+NFVAR).LT.NVARS)GO TO 1
      IF(NFLAGS.LT.0)GO TO 1
      IF(NPC.LT.0)GO TO 1
      IF(MAXC.GE.1)GO TO 2
C
C     ERROR
 1    WRITE(ODEV1,1001)IDEV,NVARS,NVN,NKW,NIVAR,NFVAR,NFLAGS,NPC,
     1 MAXC
      STOP 7401
C
C--------------------
C--------------------
C
C     INITIALIZATION
 2    NERR=0
      NCG=1
      NL=0
      NV=0
      NVALS=0
      EOF=.FALSE.
      DO 3 I=1,NVARS
 3    ENTERD(I)=.FALSE.
C
C     STORE DEFAULTS
      I1=NPC+1
      IF(NIVAR.LE.0)GO TO 950
      DO 951 I=1,NIVAR
      IVARS(I,I1)=DIVARS(I)
 951  SPECI(I,I1)=DSPECI(I)
 950  IF(NFVAR.LE.0)GO TO 952
      DO 953 I=1,NFVAR
      FVARS(I,I1)=DFVARS(I)
 953  SPECF(I,I1)=DSPECF(I)
 952  IF(NFLAGS.LE.0)GO TO 954
      DO 955 I=1,NFLAGS
 955  FLAGS(I,I1)=DFLAGS(I)
 954  CONTINUE
      GO TO 6
C
C
C--------------------
C--------------------
C
C     ENTRY TO RESUME INPUT AFTER SUPERVISOR COMMAND
      ENTRY INTP1(NERR,NCG,*,*)
C
C     SEE IF SHOULD ACCEPT MORE LINES
 6    IF(EOS.GT.0 .AND. NL.GE.EOS)GO TO 900
C
C     READ LINE
      CALL INPUT(IDEV,ECHOF,&901,&902)
      NL=NL+1
      LPTR=1
C
C
C     GET NEXT INTERESTING THING FROM LINE
 7    SPF=.FALSE.
      CALL NEXT(LPTR,NAME,BC,NUM,FNUM,F,&8,&9,&6)
C
C     HAVE VAR NAME, INDICATOR KWD, OR KWD VALUE FOR PRECEDING VAR-
C     IABLE; WHAT FOLLOWS MAY DICTATE WHICH
      L1=LPTR
      CALL NEXT(L1,TRASH,BC,NUM,FNUM,F,&11,&12,&11)
C
C     TEST HYPOTHESIS THAT NAME(8) IS VALUE FOR PRECEDING VARIABLE
 11   IF(NV.EQ.0)GO TO 13
C     BRANCH ON TYPE OF PRECEDING VARIABLE
      GO TO (13,13,14,14,14,14,13,13,13,13),TYPE
C
C     SCAN KWD RANGE FOR PRECEDING VARIABLE
 14   K1=KWR(1,NV)
      K2=KWR(2,NV)
      CALL NSRCH(NAME,KWLIST,K1,K2,K,&13)
C
C     HAVE MATCH--NAME IS VALUE FOR PRECEDING VARIABLE
      NUM=KVAL(K)
      IF(TYPE.EQ.4 .OR. TYPE.EQ.5)SPF=.TRUE.
      IF(TYPE.EQ.5)GO TO 200
      GO TO 100
C
C--------------------
C
C     MIGHT HAVE INDICATOR KWD
 13   DO 16 I=1,NVARS
      IF(ETYPE(I).NE.6)GO TO 16
C
C     SCAN KWD RANGE FOR THIS INDICATOR VARIABLE
      K1=KWR(1,I)
      K2=KWR(2,I)
      CALL NSRCH(NAME,KWLIST,K1,K2,K,&16)
C
C     HAVE MATCH--FINISH OFF PRECEDING VARIABLE
      ASSIGN 25 TO IX
      GO TO 800
 25   IF(NV.NE.I)NVALS=0
      NV=I
      J1=GVNPTR(NV)
      DO 17 J=1,8
 17   VNAME(J)=VLIST(J,J1)
      TYPE=ETYPE(NV)
      IF(TYPE.LE.0 .OR. TYPE.GT.10)GO TO 701
      NUM=KVAL(K)
      GO TO 100
 16   CONTINUE
      GO TO 18
C
C--------------------
C
C     SEE IF BREAK CHARACTER IS '='
 12   IF(BC.NE.EQUALS)GO TO 11
C
C
C     PROBABLY HAVE VARIABLE NAME
 18   K=0
 22   J=K+1
      CALL NSRCH(NAME,VLIST,J,NVN,K,&19)
C
C     WE SEEM TO HAVE A VARIABLE NAME
      I=VNUM(K)
      IF(I.LE.0 .OR. I.GT.NVARS)GO TO 700
      IF(ETYPE(I).GT.7)GO TO 22
C
C     SEE IF VARIABLE NAME HAS ALREADY BEEN USED
      IF(.NOT.ENTERD(I))GO TO 23
C
C     VAR NAME HAS ALREADY APPEARED IN THE LIST
      J=GVNPTR(I)
      WRITE(ODEV1,1006)(VLIST(I1,J),I1=1,8),NAME
      GO TO 99
C
C     FINISH PROCESSING PREVIOUS VARIABLE
 23   ASSIGN 24 TO IX
      GO TO 800
 24   NV=I
      TYPE=ETYPE(NV)
      IF(TYPE.LE.0 .OR.TYPE.GT.10)GO TO 701
      DO 20 I=1,8
 20   VNAME(I)=NAME(I)
C
C     FLAG VARIABLES ARE HANDLED HERE
      IF(TYPE.NE.7)GO TO 7
      I=VSPTR(NV)
      DO 27 J=1,NCG
      K=NPC+J
 27   FLAGS(I,K)=.TRUE.
      ENTERD(NV)=.TRUE.
      NVALS=1
      GO TO 7
C
C--------------------
C--------------------
C
C     HAVE NUMBER WHICH BETTER BE VALUE FOR PRECEDING VARIABLE
 8    IF(NV.EQ.0)GO TO 21
      GO TO (100,200,21,100,200),TYPE
C
C     MEANINGLESS NUMBER
 21   WRITE(ODEV1,1003)FNUM
      GO TO 99
C
C
C     INVESTIGATE BREAK CHARACTER
 9    CALL ATQ(BC,LPTR,&900)
      GO TO 7
C
C
C     UNRECOGNIZABLE NAME
 19   WRITE(ODEV1,1005)NAME
C
C
C     COME HERE AFTER ERROR
 99   NERR=NERR+1
C     FINISH PROCESSING OF PREVIOUS VARIABLE
      ASSIGN 7 TO IX
      GO TO 800
C
C--------------------
C--------------------
C
C     ENTER VALUE FOR INTEGER VARIABLE
 100  NVALS=NVALS+1
      ENTERD(NV)=.TRUE.
      LOC=VSPTR(NV)
      IF(NVALS.GT.1)GO TO 101
C
C     FIRST VALUE FOR VARIABLE DOES NOT GENERATE ANY NEW CONDITIONS
      IC=NPC
 106  IF(SPF)GO TO 102
C
C     STORE INTEGER VALUE
      DO 103 I=1,NCG
      I1=IC+I
      SPECI(LOC,I1)=0
 103  IVARS(LOC,I1)=NUM
      GO TO 7
C
C     STORE 'SPECIAL' VALUE
 102  DO 104 I=1,NCG
      I1=IC+I
      IVARS(LOC,I1)=0
 104  SPECI(LOC,I1)=NUM
      GO TO 7
C
C     MUST GENERATE NEW CONDITIONS
 101  NX=NPC+NCG*NVALS
      IF(NX.LE.MAXC)GO TO 105
C
C     TOO MANY CONDITIONS GENERATED
 110  J1=GVNPTR(NV)
      WRITE(ODEV1,1007)MAXC,(VLIST(J2,J1),J2=1,8),VNAME
      GO TO 99
C
C
 105  IC=NPC+(NVALS-1)*NCG
C     COPY VALUES FOR ALL VARIABLES INTO NEW CONDITIONS
      ASSIGN 106 TO IY
      GO TO 600
C
C--------------------
C--------------------
C
C     ENTER VALUE FOR FLOATING POINT VARIABLE
 200  NVALS=NVALS+1
      ENTERD(NV)=.TRUE.
      LOC=VSPTR(NV)
      IF(NVALS.GT.1)GO TO 201
C
C     FIRST VALUE FOR VARIABLE DOES NOT GENERATE ANY NEW CONDITIONS
      IC=NPC
C
 206  IF(SPF)GO TO 202
C
C     STORE FLOATING POINT VALUE
      DO 203 I=1,NCG
      I1=IC+I
      SPECF(LOC,I1)=0
 203  FVARS(LOC,I1)=FNUM
      GO TO 7
C
C     STORE 'SPECIAL' VALUE
 202  DO 204 I=1,NCG
      I1=IC+I
      FVARS(LOC,I1)=0.
 204  SPECF(LOC,I1)=NUM
      GO TO 7
C
C     MUST GENERATE NEW CONDITIONS
 201  NX=NPC+NCG*NVALS
      IF(NX.GT.MAXC)GO TO 110
      IC=NPC+(NVALS-1)*NCG
C
C     COPY VALUES FOR ALL VARIABLES INTO NEW CONDITIONS
      ASSIGN 206 TO IY
      GO TO 600
C
C--------------------
C--------------------
C
C     COPY VALUES FOR ALL VARIABLES INTO NEW CONDITIONS
 600  DO 601 I=1,NCG
      I1=NPC+I
      I2=IC+I
      IF(NIVAR.LE.0)GO TO 603
      DO 602 J=1,NIVAR
      IVARS(J,I2)=IVARS(J,I1)
 602  SPECI(J,I2)=SPECI(J,I1)
C
 603  IF(NFVAR.LE.0)GO TO 605
      DO 604 J=1,NFVAR
      FVARS(J,I2)=FVARS(J,I1)
 604  SPECF(J,I2)=SPECF(J,I1)
C
 605  IF(NFLAGS.LE.0)GO TO 601
      DO 606 J=1,NFLAGS
 606  FLAGS(J,I2)=FLAGS(J,I1)
 601  CONTINUE
      GO TO IY,(106,206)
C
C--------------------
C--------------------
C
C     END-OF-FILE
 901  EOF=.TRUE.
C
C     END-OF-STRING
C     FINISH PROCESSING OF LAST VARIABLE
 900  ASSIGN 903 TO IX
      GO TO 800
 903  IF(EOF)RETURN1
      RETURN
C
C     COMMAND ENCOUNTERED
 902  RETURN 2
C
C
C     BAD VALUE IN VNUM
 700  WRITE(ODEV1,1008)K,NAME,I
      STOP 7402
C
C     BAD VALUE IN ETYPE
 701  WRITE(ODEV1,1009)NV,NAME,TYPE
      STOP 7403
C
C--------------------
C--------------------
C
C     THIS SECTION FINISHES PROCESSING OF A VARIABLE AFTER NEW
C     VARIABLE ENCOUNTERED, EOS, OR ERROR
 800  IF(NV.EQ.0)GO TO 801
      IF(NVALS.GT.0)GO TO 802
C
C     NO VALUES ENTERED FOR LAST VARIABLE
      NERR=NERR+1
      J1=GVNPTR(NV)
      WRITE(ODEV1,1004)(VLIST(J2,J1),J2=1,8),VNAME
      NVALS=1
C
C     ADD IN EXTRA CONDITIONS, IF ANY
 802  NCG=NCG*NVALS
      NV=0
 801  NVALS=0
      GO TO IX,(25,24,7,903)
C
C--------------------
C--------------------
C
 1001 FORMAT(1H0,'BAD CALL TO INTERP'/1X,'   IDEV  NVARS',
     1 4X,'NVN	  NKW  NIVAR  NFVAR NFLAGS    NPC   MAXC'
     2 /1X,10I7)
 1003 FORMAT(1H0,'MEANINGLESS NUMBER: ', G12.5)
 1004 FORMAT(1H0,'NO VALUE GIVEN FOR VARIABLE ',8A1,' ('
     1 , 8A1, ')')
 1005 FORMAT(1H0,'UNRECOGNIZABLE NAME: ',8A1/)
 1006 FORMAT(1H0,'VARIABLE NAME ',8A1,' (',8A1,') APPEARS',
     1	' TWICE')
 1007 FORMAT(1H0,'MORE THAN',I4,' CONDITIONS GENERATED BY'/
     1 ' EXTRA VALUES FOR VARIABLE ',8A1,' (',8A1,')')
 1008 FORMAT(1H0,'BAD VALUE IN VNUM FOR NAME ',I4,2X,8A1,
     1	'  VALUE=',I8)
 1009 FORMAT(1H0,'BAD VALUE IN ETYPE FOR VARIABLE',I4,2X,8A1,
     1	'  TYPE=',I8)
      END
      SUBROUTINE JIGGLE
C     SYSTEM-DEPENDENT SUBROUTINE TO 'JIGGLE' THE RANDOM NUMBERS
C     GENERATOR BEFORE EACH RUN TO PREVENT EXACT REPEATABILITY OF
C     RESULTS.
C
      INTEGER*4 SINK
      INTEGER*4 OVNAM
      LOGICAL*4 COSTPT,DATOUT
      COMMON /IO1/NOV,SINK,OVNAM(8,6),COSTPT,DATOUT
      DIMENSION JX(2)
	CALL TIME (X,Y)
	DECODE (5,1,Y) A
  1     FORMAT (1X,F4.1)
	I=10*(A-INT(A))*A
	CALL SETRAN(I)
	RETURN
C
C
C     SYSTEM-DEPENDENT SUBROUTINE TO PRINT A LINE CONTAINING TIME
C     AND DATE IN OUTPUT HEADING FOR EACH EXPTL GROUP.
C     PURELY OPTIONAL.
C
      ENTRY TANDD
!!      INTEGER*4 SINK
!!      INTEGER*4 OVNAM
!!      LOGICAL*4 COSTPT,DATOUT
!!      COMMON /IO1/NOV,SINK,OVNAM(8,6),COSTPT,DATOUT
!!      DIMENSION JX(2)
C
	CALL TIME (IX,IY)
	CALL DATE (JX)
	WRITE (SINK,1011) JX,IX,IY
C  
 1011    FORMAT (1X,2A5,3X,2A5)
      END
      INTEGER FUNCTION JMAXX(ILINE,JJ)
C	FUNCTION RETURNS THE POSTION NUMBER(PLACE) OF THE LAST 
C	NON-BLANK CHARACTER IN THE PASSED VECTOR OR MATRIX,ILINE.
C	JJ=MAX DIMENSION OF PASSED VECTOR OR MATRIX(I*J*ETC).
      DIMENSION ILINE(1)
      DATA IBLANK/' '/
      JMAX=1
      DO 2 K=1,JJ
      IF(ILINE(K).NE.IBLANK) GO TO 1
      GO TO 2
 1    JMAX=K
 2    CONTINUE
      JMAXX=JMAX
      RETURN
      END
      SUBROUTINE BATCHK(ODEV1)
C	PROGRAM TO ASSIGN A VALUE OF 3 TO ODEV1 IF A BATCH
C	JOB IS BEING RUN, OR A VALUE OF 5 IF JOB IS AT THE TERMINAL
C	USES DEC-10 SPECFIC MACRO FUNCTION SUBROUTINE 'BATCH'
C	'BATCH' MAY OR MAY NOT BE INSTELATION SPECFIC
      INTEGER BATCH, ODEV1
      I=BATCH(I)
      IF(I.EQ.0) ODEV1=3
      IF(I.NE.0) ODEV1=5
      RETURN
      END