!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                           C
!                                     DECOMP                                C
!                                                                           C
!          THIS PROGRAM CARRIES OUT MULTIPLE STANDARDIZATION ANALYSIS       C
!          AND DEMOGRAPHIC DECOMPOSITION USING DAS GUPTA'S METHOD.          C
!          THE PROGRAM IS CONTROLLED BY A COMMAND FILE CONTAINING           C
!          COMMANDS SIMILAR TO THOSE USED IN SPSS-PC.                       C
!                                                                           C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!       THIS PROGRAM WAS DESIGNED AND TESTED USING MICROSOFT FORTRAN 4.10
!       BETWEEN 1986 AND 1989. 
!
!       POTENTIAL PROBLEMS IN USING THIS PROGRAM ON ANOTHER COMPILER ARE
!       DISCUSSED IN COMMENTS CONTAINING A STRING OF SIX ASTERISKS, AS
!       FOLLOWS: ******.  OF COURSE, DEPENDING ON THE COMPILER, YOU MAY
!       RUN INTO ADDITIONAL PROBLEMS.
! 
!       IN 2021, THE PROGRAM WAS REVISED FOR COMPLIANCE WITH FORTRAN 90. 
!       THIS REQUIRED THOUSANDS OF CHANGES, MOST OF WHICH WERE IMPLEMENTED 
!       AUTOMATICALLY BY THE NAG FORTRAN BUILDER CODE CLEANER, AND MANY OTHERS
!       IMPLEMENTED MANUALLY.
!    
!       THIS IS THE MAIN PROGRAM.  IT READS COMMANDS FROM THE COMMAND
!       FILE, CALLS SUBROUTINES TO INTERPRET THEM, AND THEN CALLS THE
!       SUBROUTINES TO EXECUTE INSTRUCTIONS.
!
    Program deco
      Implicit Integer (A-Z)
      Integer argc
      Character *3 scrnam
      Character *80 setcmd
      Character *10 varnam
      Character *35 files
      Character *5 command
      Character *35 comnam
      Character *500 list
!
!  ****** WARNING: SOME COMPILERS WILL NOT PERMIT CHARACTER VARIABLES
!         OVER 127 BYTES
!
      Real table
      Dimension setcmd(10)
      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension files(4)
      Dimension recode(200, 4)
      Dimension tratab(300, 2)
      Dimension seldes(50, 3)
      Dimension comdes(10, 5)
      Dimension brkloc(5)
      Dimension stdloc(5)
      Dimension facloc(5)
      Common /tabs/table(640000, 2), maxvars(5)
! 
!     THE ARRAYS ABOVE WILL CONTAIN THE EXECUTION INSTRUCTIONS, AFTER THEY
!     HAVE BEEN INTERPRETED.  THE INFORMATION IS ORGANIZED AS FOLLOWS:
!
!     VARDES: VARIABLE DESCRIPTIONS
!             ROWS  - VARIABLE NUMBER, IN ORDER SPECIFIED
!             COL 1 - BEGINS COLUMN
!             COL 2 - ENDS COLUMN
!             COL 3 - MINIMUM VALUE
!             COL 4 - MAXIMUM VALUE
!             COL 5 - PURPOSE
!                      1 - DEPENDENT VARIABLE
!                      2 - INDEPENDENT VARIABLE 1
!                      3 - INDEPENDENT VARIABLE 2
!                      4 - INDEPENDENT VARIABLE 3
!                      5 - INDEPENDENT VARIABLE 4
!                      6 - INDEPENDENT VARIABLE 5
!                      7 - WEIGHT
!                      8 - OTHER
!
!      VARNAM: VARIABLE NAMES, BY NUMBER
!
!      FILES:          1 - COMMAND FILE
!                      2 - DATA FILE
!                      3 - OUTPUT FILE
!                      4 - WRITE TABLE FILE
!
!      RECODE: RECODE DESCRIPTIONS
!              ROWS  - RECODE NUMBER
!              COL 1 - VARIABLE NUMBER
!              COL 2 - GE VALUE
!              COL 3 - LE VALUE
!              COL 4 - NEW VALUE
!
!      SELDES: SELECT IF DESCRIPTIONS
!              ROWS  - SELECT IF NUMBER
!              COL 1 - VARIABLE NUMBER
!              COL 2 - OPERATOR
!                      1 - EQ
!                      2 - GT
!                      3 - LT
!                      4 - GE
!                      5 - LE
!              COL 3 - CONSTANT FOR COMPARISON
!
!       VARNUM - NUMBER OF VARIABLES
!       RECNUM - NUMBER OF RECODES
!       SELNUM - NUMBER OF SELECT IFS
!       POP1D  - DEFINITION OF POPULATION 1
!       POP2D  - DEFINITION OF POPULATION 2
!       WTDIG  - NUMBER OF DIGITS THAT SHOULD BE TO THE
!                RIGHT OF THE DECIMAL POINT ON VARIABLE WEIGHT
!
!       TRATAB: TRANSFORMATION TABLE
!                ROWS  - TANSFORMATION NUMBER
!                COL 1 - TRANSFORMATION TYPE
!                        1 - SELECT IF
!                        2 - RECODE
!                        3 - COMBINE
!                COL 2 - NUMBER WITHIN TRANSFORMATION TYPE
!                         (E.G. RECNUM, SELNUM)
!
!        FIRST, CALL SUBROUTINE SETUP, WHICH CHECKS TO SEE
!        IF A 'SETUP.CMD' FILE EXISTS, AND INTERPRETS THE
!        SET PROMPT AND SET PAGENUMS COMMANDS
!
      optyp = 1
      Call setup(optyp, setcmd)
!
!   SET DEFAULTS
!
      prnum = 0
      recnum = 0
      selnum = 0
      tranum = 0
      comnum = 0
      wtdig = 0
      depdig = 0
      chkdat = 0
      chkmak = 0
      comnam = '                                   '
      files(3) = 'declis'
      files(4) = 'dectab'
!
!       OPEN OUTPUT FILE AND INTIALIZE PAGE AND LINE NUMBERS.
!
      Open (23, File=files(3))
      pg = 0
      ln = 100
      Call page(pg, ln)
      ln = ln + 1
!
!   ****** WARNING: MANY SYSTEMS WILL NOT ACCEPT FILENAME EXTENSIONS
!
!
!       WRITE INITIAL MESSAGE AND GET COMMAND FILE NAME
!
      Inquire (Unit=6, Name=scrnam)
!
!   ****** THE FOLLOWING ASSUMES THAT UNIT 6 IS SET TO FILENAME 'CON'
!          BY DEFAULT.  THIS IS USED THROUGHOUT THE PROGRAM TO SUPRESS
!          SCREEN OUTPUT IF THE 'SET SCREEN=OFF' OPTION HAS BEEN
!          INVOKED IN A 'SETUP.CMD' FILE; WHEN THE SCREEN IS SET OFF,
!          THE FILENAME ASSOCIATED WITH UNIT 6 (SCRNAM) IS RENAMED.
!          TO GET RID OF THIS FEATURE, SEARCH THIS SOURCE FILE FOR
!          'CON'.
!
!cc alternate method ... set SCRNAM to CON manually
      scrnam = 'CON'
!*                                      ** check for command line argument

      Call get_command_argument(1, comnam, argc)
      If (argc>=1) Then
        If (optyp==1 .And. scrnam=='CON') Then
          Write (*, 100)
100       Format (/, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, &
            /, /, ' _______________________________________', &
            '_______________________________________', /, /, 10X, &
            '                          DECOMP', /, &
            ' _______________________________________', &
            '_______________________________________', /, /, /, 10X, &
            '               Components Analysis Program', /, 10X, &
            '                       Version 2.0', /, /, /, 10X, &
            '        Institute for Social Research and Data Innovation', /, 10X, &
            '                  University of Minnesota', /, 10X, &
            '                   Minneapolis, MN 55455', /, /, /, /, /, 10X, &
            '  (C) COPYRIGHT 1986-2021 by Steven Ruggles;', /, 10X, &
            '      ALL RIGHTS RESERVED', /, /, /, /, /)

        End If
      Else If (optyp==1 .And. scrnam=='CON') Then
        Write (*, 110)
110     Format (/, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, /, &
          ' _______________________________________', &
          '_______________________________________', /, /, 10X, &
          '                          DECOMP', /, &
          ' _______________________________________', &
          '_______________________________________', /, /, /, 10X, &
          '               Components Analysis Program', /, 10X, &
          '                       Version 2.0', /, /, /, 10X, &
          '        Institute for Social Research and Data Innovation', /, 10X, &
          '                  University of Minnesota', /, 10X, &
          '                   Minneapolis, MN 55455', /, /, /, /, /, 10X, &
          '  (C) COPYRIGHT 1986-2021 by Steven Ruggles;', /, 10X, &
          '      ALL RIGHTS RESERVED', /, /, /, /, /)
          write(*,'(a41)', advance='no') "Enter the name of your command file --->  "
!
!  ****** WARNING: THE BACKSLASH EDIT DESCRIPTOR IS NOT ALLOWED ON
!         SOME COMPILERS.  IF YOU INTEND TO RUN THE PROGRAM IN BATCH
!         MODE, YOU MAY WANT TO GET RID OF THIS WRITE STATEMENT.
!
!  note: changed to $ edit descriptor for Sun
!        changed to advance='no' for F90 compatibility
!
        Read (*, 120) comnam
120     Format (A35)
      End If
!
!          CLEAR SCREEN AND OPEN COMMAND FILE
!
!cc          IF(SCRNAM.EQ.'CON')WRITE(*,150)
!
!   IF THE COMMAND FILE WAS NOT GIVEN, IT WILL BE TAKEN FROM
!   THE COMMAND LINE
!
      Open (3, File=comnam, Status='OLD', Iostat=iocheck)
!
!   IF THE COMMAND FILE WAS NOT FOUND, WRITE ERROR MESSAGE.
!
      If (iocheck/=0) Call error(704)
!
!   NOW GET A COMMAND.
!
130   Call getlist(command, list, long, done, pg, ln, setcmd, prnum)
!
!   CALL THE APPROPRIATE COMMAND INTERPRETER.  IF THERE IS AN ERROR,
!   EXIT THE PROGRAM.  IF THE COMMAND IS AN EXECUTABLE PROCEDURE,
!   CALL THE APPROPRIATE EXECUTION ROUTINE.
!
!   START WITH THE DATA LIST COMMAND
!
!
      If (command=='DATA ') Then
!
!   INITIALIZE RECODES, ETC. (THE DATA LIST COMMAND ERASES ANY
!           PREVIOUS DATA TRANSFORMATIONS AND DATA TABLES)
!
        recnum = 0
        selnum = 0
        tranum = 0
        comnum = 0
        wtdig = 0
        depdig = 0
        Call datlist(list, long, files, varnam, vardes, varnum, maxno)
!
!   INSERT A PERIOD FOR THE DELIMITER IN THE DATA FILE NAME, IF NECESSARY
!   (IT WAS PREVIOUSLY REMOVED TO AVOID CONFUSION WITH THE END-OF-COMMAND
!   PERIOD)
!
        perloc = index(files(2), '^')
        If (perloc>0) files(2)(perloc:perloc) = '.'
        chkdat = 1
        chkmak = 0
!
!   CALL THE RECODE COMMAND INTERPRETER IF NEEDED
!
      Else If (command=='RECOD') Then
        Call recdat(list, long, recode, varnam, recnum, tratab, tranum)
!
!   CALL THE SELECT IF COMMAND INTERPRETER IF NEEDED
!
      Else If (command=='SELEC') Then
        Call seldat(list, long, seldes, varnam, selnum, tratab, tranum)
!
!   CALL THE STANDARDIZE COMMAND INTERPRETER,
!   FIRST CHECKING FOR CORRECT COMMAND SEQUENCE.
!   THEN CALL STANDAR TO EXECUTE THE STANDARDIZATION.
!
      Else If (command=='STAND') Then
        Call check(chkdat, chkmak)
        Call stddat(list, long, varnam, vardes, varnum, numbrk, brkloc, &
          numstd, stdloc, stpop, stcat, form, wrex)
        Call standar(varnam, vardes, numbrk, brkloc, numstd, stdloc, stpop, &
          stcat, form, wrex, files, wtdig, wtdiv, pg, ln)
!cc            IF(SCRNAM.EQ.'CON')WRITE(*,150)
!
!   CALL THE DECOMPOSE COMMAND INTERPRETER,
!   FIRST CHECKING FOR CORRECT COMMAND SEQUENCE.
!   THEN CALL GUPTA TO EXECUTE THE STANDARDIZATION.
!
      Else If (command=='DECOM') Then
        Call check(chkdat, chkmak)
        Call decomp(list, long, varnam, vardes, varnum, poploc, pop1d, pop2d, &
          numfac, facloc, wrex)
        Call gupta(varnam, vardes, pop1d, pop2d, poploc, numfac, facloc, &
          files, wtdig, wtdiv, wrex, pg, ln)
!cc            IF(SCRNAM.EQ.'CON')WRITE(*,150)
!
!   CALL THE WEIGHT COMMAND INTERPRETER IF NEEDED
!
      Else If (command=='WEIGH') Then
        Call weight(list, long, varnam, vardes, varnum, wtdig)
!
!   CALL THE COMBINE COMMAND INTERPRETER IF NEEDED
!
      Else If (command=='COMBI') Then
        Call combdat(list, long, comdes, varnam, vardes, varnum, comnum, &
          tratab, tranum, pg, ln)
!
!   CALL THE MAKETAB COMMAND INTERPRETER, FIRST CHECKING FOR
!   COMMAND SEQUENCE.  THEN CREATE THE TABLE WITH SUBROUTINE TAB.
!
      Else If (command=='MAKET') Then
        chkmak = 1
        Call check(chkdat, chkmak)
        Call make(list, long, varnam, vardes, varnum, depdig, wrtab)
        Call tab(vardes, recode, seldes, comdes, tratab, files, maxno, varnam, &
          tranum, varnum, wtdig, depdig, wrtab, wtdiv, pg, ln)
!cc            IF(SCRNAM.EQ.'CON')WRITE(*,150)
!
!   CALL THE SET LISTING COMMAND INTERPRETER IF NEEDED, THEN
!   CHANGE THE OUTPUT FILE AND RESET PAGE NUMBER.
!
      Else If (command=='SET L') Then
        Call set(command, list, files, pg, ln)
        Close (23)
        Open (23, File=files(3))
        ln = 100
        pg = 0
        Call page(pg, ln)
!
!   CALL THE SET RESULTS COMMAND INTERPRETER IF NEEDED.
!
      Else If (command=='SET R') Then
        Call set(command, list, files, pg, ln)
!
!   IF SET COMMAND INVALID OR NO COMMAND FOUND, WRITE AN ERROR MESSAGE.
!
      Else If (command(1:3)=='SET') Then
        Call page(pg, ln)
!
!  ****** WARNING: THE FOLLOWING STATEMENT WRITES THE BELL CHARACTER.
!                  THIS WILL NOT WORK ON ALL SYSTEMS.
!
!cc          WRITE(6,*)CHAR(007)
        Write (23, 140)
        Write (*, 140)
140     Format (/, 1X, 'ERROR 4 in SET command: Command ignored', /)
      Else If (command/='NULL ') Then
        If (command(1:1)==' ') Then
          Call error(705)
        Else
          Write (23, 150) command
          If (scrnam=='CON') Write (*, 150) command
150       Format (/, /, ' ERROR 3: Command ', A5, ' not recognized')
          Call faterr
        End If
      End If
!
!  IF THERE ARE MORE COMMANDS, GO GET THEM AND EXECUTE THEM.
!
      If (done/=1) Go To 130
!
!       WRITE TERMINATION MESSAGE, DELETE WORKING TABLE, AND END.
!
      If (scrnam=='CON') Then
!cc          WRITE(*,150)
        Write (*, 160)
      End If
      Close (22, Status='DELETE')
!cc        CLOSE(6,STATUS='DELETE')
      Endfile 23
!cc 150    FORMAT(12(1H0/))
160   Format (' No errors detected')
    End Program
!
!
!
    Subroutine check(chkdat, chkmak)
!
!       THIS SUBROUTINE CHECKS FOR VALID STATEMENT ORDER.
!
      Implicit Integer (A-Z)

      If (chkdat==0) Then
        Write (*, 100)
        Write (23, 100)
        Call faterr
      End If
      If (chkmak==0) Then
        Write (*, 110)
        Write (23, 110)
        Call faterr
      End If
100   Format (/, /, 1X, 'ERROR 1: DATA LIST statement not found')
110   Format (/, /, 1X, 'ERROR 2: MAKETAB statement not found')
      Return
    End Subroutine
!
!
!
    Subroutine standar(varnam, vardes, numbrk, brkloc, numstd, stdloc, stpop, &
      stcat, form, wrex, files, wtdig, wtdiv, pg, ln)
!
!      THIS SUBROUTINE EXECUTES THE STANDARDIZE PROCEDURE
!
      Implicit Integer (A-Z)
      Integer indval(6)
      Real table, total, results, tnum, tden, oldtot
      Character *10 varnam
      Character *35 files
      Dimension files(4)
      Dimension results(250, 3)
      Dimension table(640000, 2)
      Dimension maxvars(5), brkloc(5), stdloc(5)
      Dimension varnam(30)
      Dimension vardes(30, 5)
      Dimension newpos(5)
      Dimension wrform(7, 2)
      Common /tabs/table, maxvars
!
!      RUN THE PROGRAM SEPARATELY FOR EACH BREAKDOWN VARIABLE
!
      Do brk = 1, numbrk
!
!      INITIALIZE WORKING TABLE AND RESULTS TABLE
!
        Do i = 1, 640000
          table(i, 1) = 0.0
          table(i, 2) = 0.0
        End Do

        Do i = 1, 250
          results(i, 1) = 0.0
          results(i, 2) = 0.0
          results(i, 3) = 0.0
        End Do
!
!      DETERMINE WHICH VARIABLE IS THE CURRENT BREAKDOWN VARIABLE
!      THE LOCATION OF EACH BREAKDOWN VARIABLE IS CONTAINED IN ARRAY
!      BRKLOC.
!
        brkpos = vardes(brkloc(brk), 5) - 1
        If (brkpos<1 .Or. brkpos>5) Then
          Call error(303)
        End If
!
!      CALCULATE MAXIMUM VALUE FOR BREAKDOWN VARIABLE AND INITIALIZE
!      THE MAXIMUM OF CONTROL VARIABLES TO ONE.
!
!      NEWPOS WILL RELATE THE ORIGINAL POSITION OF THE VARIABLES TO THEIR
!      POSITION IN THE WORKING TABLE.
!
        Do i = 1, 5
          maxvars(i) = 1
          newpos(i) = 6
        End Do
        total = 0.0
        indval(6) = 1
        newpos(1) = brkpos
        maxvars(1) = vardes(brkloc(brk), 4) - vardes(brkloc(brk), 3) + 2
!
!      IF A PARTICULAR BREAKDOWN CATEGORY HAS BEEN SPECIFIED AS THE
!      STANDARD, ADJUST ITS VALUE.  IF THE CATEGORY FALLS OUTSIDE
!      THE RANGE OF THE VARIABLE, WRITE AN ERROR MESSAGE.
!
        If (stpop==3) Then
          stval = stcat - vardes(brkloc(brk), 3) + 2
          If (stval<2 .Or. stval>maxvars(1)) Then
            Call error(308)
          End If
        End If
!
!      THE FOLLOWING SECTION INSERTS THE POSITION FOR EACH OF THE
!      CONTROL VARIABLES (WHOSE LOCATION HAS BEEN TRANSFERRED IN
!      ARRAY STDLOC), CALCULATES THEIR MAXIMUM VALUE, AND MAKES SURE
!      THAT THE CONTROL VARIABLE IS NOT THE SAME AS THE BREAKDOWN VARIABLE.
!      (IF IT IS, WE JUST SKIP IT).
!
        next = 1
        nstnow = numstd
        Do i = 1, numstd
          If (stdloc(i)==brkloc(brk)) Then
            nstnow = nstnow - 1
            If (nstnow==0) Go To 360
            Go To 100
          End If
          next = next + 1
          newpos(next) = vardes(stdloc(i), 5) - 1
          If (newpos(next)<1 .Or. newpos(next)>5) Then
            Call error(303)
          End If
          maxvars(next) = vardes(stdloc(i), 4) - vardes(stdloc(i), 3) + 1
100     End Do
!
!       NOW WE ARE READY TO READ IN THE PRETABULATED DATA, AND STICK IT
!       IN THE APPROPRIATE PLACE.  WE USE INDVAL AND NEWPOS TO GET THEM
!       THERE.  FIRST, REWIND THE FILE WITH THE PRETABULATED DATA.
!
        Rewind 22
        Do i = 1, 640000
          Read (22, End=110) indval(1), indval(2), indval(3), indval(4), &
            indval(5), tnum, tden
!
!       THE FIRST POSITION OF OUR WORKING TABLE WILL CONTAIN THE
!       CURRENT BREAKDOWN VARIABLE.  LEAVE AN EXTRA SPACE AT THE
!       BEGINNING TO CONTAIN MARGINALS.
!
          a = indval(newpos(1)) + 1
          b = indval(newpos(2))
          c = indval(newpos(3))
          d = indval(newpos(4))
          e = indval(newpos(5))
!
!       CALCULATE THE INDEX AND INCREMENT THE TABLE.
!
          index = ind(a, b, c, d, e)
          table(index, 1) = table(index, 1) + tnum
          table(index, 2) = table(index, 2) + tden
!
!       IF STPOP IS EQUAL TO 2, TABULATE STANDARD POPULATION AS THE
!       COMBINATION OF ALL BREAKDOWN CATEGORIES.
!
          If (stpop==2) Then
            index = ind(1, b, c, d, e)
            table(index, 2) = table(index, 2) + tden
!
!         IF STPOP IS SPECIFIED AS 3, CALCULATE STANDARD POPULATION
!         AS SPECIFIED BREAKDOWN CATEGORY
!
          Else If (stpop==3 .And. stval==a) Then
            index = ind(1, b, c, d, e)
            table(index, 2) = table(index, 2) + tden
          End If
!
!       TEMPORARILY PUT THE TOTAL CASES OF EACH BREAKDOWN CATEGORY
!       IN ARRAY RESULTS, AND ADD UP TOTAL CASES.
!
          results(a, 1) = results(a, 1) + tden
          total = total + tden
        End Do
110     Continue
!
!       NOW EXCLUDE ANY MISSING CELLS. SAVE THE ORIGINAL TOTAL CASES.
!
        oldtot = total
        exclu = 0
!
!       IF THE "WRITE EXCLUDED" SUBCOMMAND WAS ISSUED, PREPARE TO DO SO.
!
        If (wrex>0) Then
          Call exprep(vardes, wrform, files, brkloc(brk), wtdig, stdloc, &
            numstd, nvars)
        End If
!
!       NOW IT IS TIME TO WRITE THE HEADER FOR THE TABLE.
!       CHECK TO SEE IF THERE IS SUFFICIENT ROOM ON THE PAGE TO FIT
!       A SIGNIFICANT CHUNK OF THE DECOMP TABLE.  IF NOT, START A NEW
!       PAGE. (ONLY DO THIS FOR THE FIRST BREAKDOWN VARIABLE).
!
        If (brk==1) Then
          If (ln>20) ln = 100
          Call page(pg, ln)
          ln = ln + 9
          If (ln>10) Then
            ln = ln + 1
            Write (23, 350)
          End If
          Write (23, 120)
120       Format (/, 7X, 'M U L T I P L E   S T A N D A R D I Z A T I O N   ', &
            'A N A L Y S I S')
        End If
!
!       THE FOLLOWING SECTION DEALS WITH THE PROBLEM OF EXCLUDED
!       CELLS.  A CELL IS EXCLUDED IF IT IS EMPTY IN ANY OF THE
!       BREAKDOWN CATEGORIES.
!
!       FOR EACH COMBINATION OF CONTROL CHARACTERISTICS, COUNT
!       THE NUMBER OF CELLS WITH ZERO CASES.
!
        Do b = 1, maxvars(2)
          Do c = 1, maxvars(3)
            Do d = 1, maxvars(4)
              Do e = 1, maxvars(5)
                zeros = 0
                Do a = 2, maxvars(1)
                  If (results(a,1)<=0) Go To 130
                  index = ind(a, b, c, d, e)
                  If (table(index,2)<=0) zeros = zeros + 1
130             End Do
!
!       IF THERE ARE CELLS WITH ZEROS, TAKE CARE OF IT.
!
                If (zeros>0) Then
                  Do a = 2, maxvars(1)
                    index = ind(a, b, c, d, e)
!
!       IF THE "WRITE EXCLUDED" OPTION WAS SPECIFIED, DETERMINE THE
!       UNADJUSTED VALUES OF EACH VARIABLE AND CALL EXDATA TO
!       ACTUALLY WRITE THE EXCLUDED CELLS.  ONLY DO THIS, HOWEVER,
!       IF THERE ARE CASES IN THIS PARTICULAR CELL.
!
                    If (wrex>0 .And. table(index,2)>0) Then
                      i = a + vardes(brkloc(brk), 3) - 2
                      j = b + vardes(stdloc(1), 3) - 1
                      k = c + vardes(stdloc(2), 3) - 1
                      l = d + vardes(stdloc(3), 3) - 1
                      m = e + vardes(stdloc(4), 3) - 1
                      Call exdata(wrform, i, j, k, l, m, nvars, wtdig, wtdiv, &
                        table(index,1), table(index,2))
                    End If
!
!       NOW COUNT THE NUMBER OF EXCLUDED CASES, ADJUST THE TOTAL
!       NUMBER OF CASES, AND ZERO OUT THE OFFENDING CELLS.
!
! INT4 WAS HERE
                    exclu = exclu + int(table(index,2))
                    total = total - table(index, 2)
                    results(a, 1) = results(a, 1) - table(index, 2)
                    table(index, 1) = 0.0
                    table(index, 2) = 0.0
                  End Do
!
!       FINALLY, ZERO OUT THE MARGINALS.
!
                  index1 = ind(1, b, c, d, e)
                  table(index1, 1) = 0.0
                  table(index1, 2) = 0.0
                End If
              End Do
            End Do
          End Do
        End Do
!       CALCULATE THE STANDARD POPULATION, DEPENDING ON WHICH
!       STANDARD HAS BEEN SPECIFIED.
!
!       IF STDPOP IS 1, THEN THE STANDARD IS THE AVERAGE OF THE
!       PROPORTIONS OF THE POPULATION IN THE CELL OVER ALL
!       BREAKDOWN CATEGORIES.  THIS IS CALCULATED AS THE SUM OF THE
!       PROPORTIONS IN EACH CATEGORY DIVIDED BY THE NUMBER OF CATEGORIES.
!
        If (stpop==1) Then
!
!       FIRST WE NEED TO KNOW HOW MANY CATEGORIES OF THE BREAKDOWN
!       VARIABLE ARE ACTUALLY BEING USED (THAT IS, HOW MANY OF THEM HAVE
!       CASES).
!
          ncats = 0
          Do a = 2, maxvars(1)
            If (results(a,1)>0) ncats = ncats + 1
          End Do
!
          Do b = 1, maxvars(2)
            Do c = 1, maxvars(3)
              Do d = 1, maxvars(4)
                Do e = 1, maxvars(5)
                  index1 = ind(1, b, c, d, e)
                  Do a = 2, maxvars(1)
                    index2 = ind(a, b, c, d, e)
                    If (results(a,1)<=0 .Or. maxvars(1)<2) Go To 140
                    table(index1, 2) = table(index1, 2) + &
                      (table(index2,2)/results(a,1))/float(ncats)
140               End Do
                End Do
              End Do
            End Do
          End Do
!
!       IF STPOP IS SPECIFIED AS 2, CALCULATE STANDARD POPULATION
!       AS THE COMBINATION OF ALL BREAKDOWN CATEGORIES
!
        Else If (stpop==2) Then
          Do b = 1, maxvars(2)
            Do c = 1, maxvars(3)
              Do d = 1, maxvars(4)
                Do e = 1, maxvars(5)
                  index1 = ind(1, b, c, d, e)
                  If (total<=0) Go To 150
                  table(index1, 2) = table(index1, 2)/total
150             End Do
              End Do
            End Do
          End Do
!
!         IF STPOP IS SPECIFIED AS 3, CALCULATE STANDARD POPULATION
!         AS SPECIFIED BREAKDOWN CATEGORY
!
        Else If (stpop==3) Then
          If (results(stval,1)<=0) Then
            Call error(309)
          End If
          Do b = 1, maxvars(2)
            Do c = 1, maxvars(3)
              Do d = 1, maxvars(4)
                Do e = 1, maxvars(5)
                  index1 = ind(1, b, c, d, e)
                  If (total<=0) Go To 160
                  table(index1, 2) = table(index1, 2)/results(stval, 1)
160             End Do
              End Do
            End Do
          End Do
        End If
!
!       RE-INITIALIZE RESULTS TABLE AND TOTAL.
!       TOTAL NOW WILL BE USED TO REPRESENT THE SUM OF THE STANDARD
!       POPULATION.
!       RESULTS WILL BE USED, LOGICALLY ENOUGH, TO HOLD THE RESULTS,
!       AS FOLLOWS:
!
!       COLUMN 1 OF RESULTS: STANDARDIZED
!       COLUMN 2             RAW MEANS
!       COLUMN 3             NUMBER OF CASES
!
        Do a = 1, 250
          results(a, 1) = 0.0
        End Do
        total = 0.0
!
!       THE GUTS ARE IN THE FOLLOWING LOOP.
!
        Do b = 1, maxvars(2)
          Do c = 1, maxvars(3)
            Do d = 1, maxvars(4)
              Do e = 1, maxvars(5)
!
!       INDEX1 REFERS TO THE STANDARD POPULATION
!
                index1 = ind(1, b, c, d, e)
                total = total + table(index1, 2)
                Do a = 2, maxvars(1)
!
!       INDEX2 REFERS TO THE POPULATIONS TO BE STANDARDIZED
!
                  index2 = ind(a, b, c, d, e)
!
!       CONVERT SUM OF NUMERATORS TO RATES OR MEANS
!
                  If (table(index2,2)>0.01) Then
                    table(index2, 1) = table(index2, 1)/table(index2, 2)
                  End If
!
!       THE STANDARDIZED RESULTS ARE THE SUM OF CELLS, WEIGHTED BY
!       THE STANDARD POPULATION.
!       THE UNSTANDARDIZED RESULTS ARE WEIGHTED BY EACH BREAKDOWN
!       POPULATION.
!
                  results(a, 1) = results(a, 1) + table(index2, 1)*table( &
                    index1, 2)
                  results(a, 2) = results(a, 2) + table(index2, 1)*table( &
                    index2, 2)
                  results(a, 3) = results(a, 3) + table(index2, 2)
                End Do
              End Do
            End Do
          End Do
        End Do
!
!      ADD UP THE TOTALS
!
        Do a = 2, maxvars(1)
          results(1, 1) = results(1, 1) + results(a, 1)
          results(1, 2) = results(1, 2) + results(a, 2)
          results(1, 3) = results(1, 3) + results(a, 3)
        End Do
!
!      FIND THE DEPENDENT VARIABLE LOCATION
!
        Do a = 1, 30
          If (vardes(a,5)==1) deploc = a
        End Do
!
!      WRITE THE SPECIFICATIONS OF THIS RUN.
!
        Write (23, 350)
        If (brk==1) Then
          Write (23, 170) varnam(deploc)
170       Format (/, ' Standardization of ', A10)
          If (stpop==2) Then
            Write (23, 180)
180         Format (/, ' Standard is the total combined population')
          Else If (stpop==1) Then
            Write (23, 190)
190         Format (/, ' Standard is the average of breakdown categories')
          Else If (stpop==3) Then
            Write (23, 200) stcat
200         Format (/, ' Standard is breakdown category', I3)
          End If
          If (form==1) Then
            Write (23, 210)
210         Format (' Results expressed as means of dependent variable')
          Else If (form==2) Then
            Write (23, 220)
220         Format (' Results expressed as percentages')
          Else If (form==3) Then
            Write (23, 230)
230         Format (' Results expressed as deviations from grand mean')
          End If
        End If
!
!      WRITE BREAKDOWN VARIABLE.
!
240     Format (/, ' Breakdown by ', A10)
        size = 13 + numstd + maxvars(1)
        oldtot = float(exclu)/oldtot*100.0
!
!      FIGURE OUT HOW MUCH SPACE WILL BE NEEDED, DEPENDING ON ANY
!      EXCLUDED DATA MESSAGES THAT MAY BE WRITTEN.
!
        exlin = 1
        If (oldtot>0.0 .And. oldtot<5.0) Then
          exlin = 2
        Else
          exlin = 4
        End If
        size = size + exlin
        If (size<50) Then
          ln = ln + size - 1
        Else
          ln = ln + 9 + numstd
        End If
!
!      IF NECESSARY, WE WILL START A FRESH PAGE.
!
        Call page(pg, ln)
        Write (23, 240) varnam(brkloc(brk))
        Write (23, 260)
250     Format (14X, A10)
260     Format (/, ' Controlling for:')
        Do i = 1, numstd
          If (brkloc(brk)/=stdloc(i)) Then
            Write (23, 250) varnam(stdloc(i))
          End If
        End Do
!
!      CALCULATE OVERALL (GRAND) UNSTANDARDIZED MEAN
!
        If (results(1,3)>0.0) results(1, 2) = results(1, 2)/results(1, 3)
!
!      IF RESULTS ARE SUPPOSED TO BE EXPRESSED AS PERCENTAGES, ADJUST.
!
        If (form==2) Then
          results(1, 2) = results(1, 2)*100.0
        End If
! INT4 WAS HERE
        num = int(results(1,3))
!
!      WRITE GRAND MEAN AND NUMBER OF CASES.
!
        Write (23, 270) results(1, 2), num
270     Format (/, ' Grand Mean:', F10.4, '        Total cases included:', I8)
!
!      WRITE SUBHEAD
!
        Write (23, 280)
280     Format (/, '   Category of      Unstandardized      Standardized')
        Write (23, 290) varnam(brkloc(brk))
290     Format ('   ', A10, '           Mean                Mean         ', &
          '       N', /)
!
!      CHECK TO MAKE SURE WE HAD VALID CASES
!
        If (total<=0) Then
          Call error(310)
        End If
!
!      FIGURE OUT NECESSARY ADJUSTMENT OF BREAKDOWN CATEGORY VALUES
!
        minbrk = vardes(brkloc(brk), 3) - 1
!
!       WRITE RESULTS, CALCULATING THEM ON THE FLY
!
        Do i = 2, maxvars(1)
          results(i, 1) = results(i, 1)/total
          If (results(i,3)>0.0) Then
            results(i, 2) = results(i, 2)/results(i, 3)
          End If
! INT4 WAS HERE
          num = int(results(i,3))
!
!       ADJUST THE OUTPUT FORMAT IF NEEDED.
!
          If (form==2) Then
            results(i, 2) = results(i, 2)*100.0
            results(i, 1) = results(i, 1)*100.0
          Else If (form==3) Then
            results(i, 2) = results(i, 2) - results(1, 2)
            results(i, 1) = results(i, 1) - results(1, 2)
          End If
!
!       CHECK THE SPACE LEFT ON THE PAGE, AND START A NEW ONE IF NEEDED.
!
          If (size>=50) Then
            Call page(pg, ln)
            If (ln==1) Then
              ln = ln + 6
              Write (23, 280)
              Write (23, 290) varnam(brkloc(brk))
            End If
          End If
          Write (23, 300) i - 1 + minbrk, results(i, 2), results(i, 1), num
        End Do
300     Format (7X, I3, 9X, 2(F10.4,10X), I8)
!
!       CHECK AGAIN.
!
        If (size>=50) Then
          ln = ln + exlin
          Call page(pg, ln)
        End If
        Write (23, 350)
!
!       NOW WRITE TRAILING MESSAGES AND HINTS ON EXCLUDED CELLS.
!
        If (oldtot<=0) Then
          Write (23, 310)
310       Format (/, 7X, 'No cases excluded on grounds of empty cells')
        Else If (oldtot<5.0) Then
          Write (23, 320) exclu, oldtot
320       Format (/, 7X, 'WARNING:', I9, ' cases excluded because of empty', &
            ' Cells ', /, 26X, '(', F4.1, ' percent of total)')
        Else If (oldtot<33.3) Then
          Write (23, 330) exclu, oldtot
330       Format (/, 7X, 'DANGER:', I9, ' cases excluded because of empty', &
            ' cells ', /, 26X, '(', F4.1, ' percent of total)', /, /, 7X, &
            'SUGGESTION: recode variables into fewer categories,', /, 7X, &
            '            eliminating', ' small populations.')
        Else
          Write (23, 340) exclu, oldtot
340       Format (/, 7X, 'DISASTER:', I9, ' cases excluded because of empty', &
            ' cells ', /, 26X, '(', F4.1, ' percent of total)', /, /, 7X, &
            'reconceptualize your problem or get more data.')
        End If
!
!       WRITE CODEBOOK FOR EXCLUDED CASES.
!
        If (wrex>0 .And. oldtot>0.0) Then
          Call excode(varnam, wrform, files, brkloc(brk), stdloc, numstd, pg, &
            ln)
        End If
        If (numbrk>1 .And. brk==numbrk) Write (23, 350)
350     Format ('________________________________________', &
          '_______________________________________')

360   End Do
      ln = 100
      Return
    End Subroutine

    Subroutine gupta(varnam, vardes, pop1d, pop2d, poploc, numfac, facloc, &
      files, wtdig, wtdiv, wrex, pg, ln)
!
!      THIS SUBROUTINE CONTROLS THE DECOMP PROCEDURE.
!      MUCH OF THE ACTUAL WORK IS CARRIED OUT BY SUBROUTINES;
!      E.G., THE ACTUAL GUTS OF THE PROGRAM (CARRIED OUT BY
!      THE SUBROUTINE GUTS) A LITTLE BIT OF PRELIMINARY SETUP
!      (CARRIED OUT BY THE SUBROUTINE GUPSET), AND THE DESCRIPTIVE
!      TABULATION (CARRIED OUT BY SUBROUTINE DESTAB).
!
      Implicit Integer (A-Z)
      Character *35 files
      Character *3 scrnam
      Integer indval(6)
      Real effect(4), tot1, tot2, comb, rate, dif, ix
      Real miss1, miss2, num1, num2
      Real table, n1, n2, n1t, n2t, t1, t2
      Real tnum, tden
      Character *10 varnam
      Dimension files(4)
      Dimension table(640000, 2)
      Dimension maxvars(5), facloc(5)
      Dimension varnam(30)
      Dimension vardes(30, 5)
      Dimension newpos(5)
      Dimension wrform(7, 2)
      Common /tabs/table, maxvars
!
!      INITIALIZE STUFF
!
!cc        INQUIRE(UNIT=6,NAME=SCRNAM)
      scrnam = 'CON'
      tot1 = 0.0
      tot2 = 0.0
      comb = 0.0
      rate = 0.0
      Do i = 1, 640000
        table(i, 1) = 0.0
        table(i, 2) = 0.0
      End Do
      Do i = 1, 4
        effect(i) = 0.0
      End Do
      indval(6) = 1
!
!       CALL SETUP ROUTINE
!
      Call gupset(vardes, pop1d, pop2d, poploc, numfac, newpos, facloc, &
        poppos, popa, popb)
!
!       REWIND THE FILE CONTAINING THE TABULATED DATA.
!
      Rewind 22
!
!       READ THE DATA FROM THE FILE.
!       AS WE GO, WE WILL REORGANIZE IT AND TABULATE THE MARGINALS.
!
!       FIRST, STICK THE VARIABLE VALUES FOR EACH CELL IN THE ARRAY
!       INDVAL AND THE NUMERATOR OF THE CELL IN TNUM AND THE
!       DENOMINATOR IN TDEN.
!
      Do i = 1, 640000
        Read (22, End=110) indval(1), indval(2), indval(3), indval(4), &
          indval(5), tnum, tden
!
!       IN OUR NEW TABLE, THE COMPARISON VARIABLE WILL BE LOCATED
!       IN THE FIRST POSITION.  THE ARRAY NEWPOS (CONSTRUCTED IN
!       SUBROUTINE GUPSET) RELATES THE POSITION OF THE VARIABLE IN
!       THE MAKETAB TABLE TO ITS NEW POSITION IN THE DECOMP TABLE.
!
        a = 0
        pop = indval(newpos(1))
!
!       POPA IS THE VALUE OF THE FIRST COMPARISON CATEGORY SPECIFIED
!       IN THE DECOMP COMMAND, AND POPB IS THE VALUE OF THE SECOND.
!
        If (pop==popa) a = 1
        If (pop==popb) a = 2
!
!       SKIP THIS CELL IF IT DOES NOT REFER TO EITHER COMPARISON CATEGORY.
!
        If (a==0) Go To 100
!
!       THE FACTORS GO IN VARIABLES B THROUGH E.  ADD ONE TO EACH VALUE
!       TO MAKE ROOM FOR THE MARGINALS.
!
        b = indval(newpos(2)) + 1
        c = indval(newpos(3)) + 1
        d = indval(newpos(4)) + 1
        e = indval(newpos(5)) + 1
!
!       GET THE INDEX, AND TABULATE.
!
        index = ind(a, b, c, d, e)
        table(index, 1) = table(index, 1) + tnum
        table(index, 2) = table(index, 2) + tden
!
!       THE MARGINALS WILL BE TABULATED BY SUBROUTINE MARGINS.
!
        Call margins(a, b, c, d, e, tden, tnum)
100   End Do
110   Continue

!
!       STICK THE TOTALS OF EACH POPULATION IN N1T AND N2T
!
      n1t = table(ind(1,1,1,1,1), 2)
      n2t = table(ind(2,1,1,1,1), 2)
!
!       IF THE "WRITE EXCLUDED" SUBCOMMAND WAS ISSUED, GET READY
!       FOR IT BY CALLING THE SUBROUTINE EXPREP.
!
      If (wrex>0) Then
        Call exprep(vardes, wrform, files, poploc, wtdig, facloc, numfac, &
          nvars)
      End If
!
!       WRITE HEADER.  START A NEW PAGE IF THERE ISN'T MUCH ROOM
!       LEFT ON THE CURRENT ONE.
!
      If (ln>20) ln = 100
      Call page(pg, ln)
      ln = ln + 14
!
!       IF WE DID NOT JUST START A NEW PAGE, WRITE A LINE ACROSS
!       THE PAGE.
!
      If (ln>15) Then
        Write (23, 150)
        ln = ln + 2
      End If
!
!       NOW WRITE TITLE.
!
      Write (23, 120)
120   Format (/, 18X, 'D E C O M P O S I T I O N   ', 'A N A L Y S I S')
!
!       INITIALIZE MISSING (EXCLUDED) NUMBERS AND INCLUDED NUMBERS.
!
      miss1 = 0
      miss2 = 0
      num1 = 0
      num2 = 0
!
!       NOW TAKE CARE OF MISSING CELLS PROBLEMS.  ANY CELL NOT
!       REPRESENTED IN BOTH POPULATIONS WILL BE IGNORED.
!
      Do i = 2, maxvars(2)
        Do j = 2, maxvars(3)
          Do k = 2, maxvars(4)
            Do l = 2, maxvars(5)
              n1 = table(ind(1,i,j,k,l), 2)
              n2 = table(ind(2,i,j,k,l), 2)
              t1 = table(ind(1,i,j,k,l), 1)
              t2 = table(ind(2,i,j,k,l), 1)
!
!       FIRST, WE DEAL WITH THE CELLS WHERE THERE ARE SOME CASES IN
!       POPULATION A BUT NO CASES IN POPULATION B.
!
              If (n1>0 .And. n2<=0) Then
                If (wrex>0) Then
!
!         IF THE "WRITE EXCLUDED" COMMAND WAS SPECIFIED, DO SO.
!         FIRST FIGURE OUT WHAT THE ORIGINAL VALUES WERE.
!
                  a = i + vardes(facloc(1), 3) - 2
                  b = j + vardes(facloc(2), 3) - 2
                  c = k + vardes(facloc(3), 3) - 2
                  d = l + vardes(facloc(4), 3) - 2
                  Call exdata(wrform, pop1d, a, b, c, d, nvars, wtdig, wtdiv, &
                    t1, n1)
                End If
!
!         COUNT THE NUMBER EXCLUDED FROM POPULATION.
!
                miss1 = miss1 + n1
!
!         ADJUST THE TOTAL N OF THE POPULATION
!
                n1t = n1t - n1
!
!         NOW WE ADJUST ALL OF THE MARGINAL N'S
!
                Call adjust(1, i, j, k, l, n1, t1)
!
!       THE PROCESS IS REPEATED FOR THE CELLS THAT HAVE CASES IN
!       POPULATION B AND NO COMPARABLE CASES IN POPULATION A.
!
              Else If (n2>0 .And. n1<=0) Then
                If (wrex>0) Then
                  a = i + vardes(facloc(1), 3) - 2
                  b = j + vardes(facloc(2), 3) - 2
                  c = k + vardes(facloc(3), 3) - 2
                  d = l + vardes(facloc(4), 3) - 2
                  Call exdata(wrform, pop2d, a, b, c, d, nvars, wtdig, wtdiv, &
                    t2, n2)
                End If
                n2t = n2t - n2
                miss2 = miss2 + n2
                Call adjust(2, i, j, k, l, n2, t2)
              End If
            End Do
          End Do
        End Do
      End Do
!
!       WE ARE READY TO CALCULATE MEANS AND AVERAGES OF DEPENDENT
!       FIRST WRITE A MESSAGE TO THE SCREEN TO LET THE USER KNOW
!       THAT IT'S GOING TO BE A WHILE.
!
      If (scrnam=='CON') Write (*, 130)
130   Format (/, /, 3X, ' CALCULATING DECOMPOSITION . . . ', /, /)
!
!       PASS THROUGH ALL CELLS
!
      Do j = 2, maxvars(2)
        Do k = 2, maxvars(3)
          Do l = 2, maxvars(4)
            Do m = 2, maxvars(5)
!
!       CALCULATE THE MEAN OF THE DEPENDENT
!
              Do i = 1, 2
                index = ind(i, j, k, l, m)

                If (table(index,2)>0.0001) Then
                  table(index, 1) = table(index, 1)/table(index, 2)
                End If
              End Do
!
!       CALCULATE THE AVERAGE OF THE TWO POPULATIONS, AND STICK
!       IT IN THE SPACE PREVIOUSLY RESERVED BY SUBROUTINE GUPSET.
!
              pop1 = ind(1, j, k, l, m)
              pop2 = ind(2, j, k, l, m)
              aver = ind(3, j, k, l, m)
              table(aver, 1) = (table(pop1,1)+table(pop2,1))/2.0
            End Do
          End Do
        End Do
      End Do
!
!      THE SUBROUTINE GUTS DOES THE ACTUAL DECOMPOSITION, NOW THAT
!      WE HAVE EVERYTHING NEATLY SET UP.  THE RESULTS COME OUT IN
!      COMB, RATE, AND EFFECT.
!
      Call guts(n1t, n2t, num1, num2, tot1, tot2, comb, rate, effect)

!
!      FIND THE LOCATION OF THE DEPENDENT VARIABLE.
!
      Do a = 1, 30
        If (vardes(a,5)==1) deploc = a
      End Do
!
!      CALCULATE TOTAL POPULATION DIFFERENCE.
!
      dif = tot1 - tot2
!
!       CONVERT N'S TO INTEGERS.
!
      nu1 = int(num1)
      nu2 = int(num2)
!
!       IF EITHER POPULATION HAS NO CASES, WE'VE GOT A PROBLEM
!
      If (nu1==0 .Or. nu2==0) Then
        Call error(407)
      End If
!
!       CONVERT MISSING N'S TO INTEGERS.
!
      mi1 = int(miss1)
      mi2 = int(miss2)
!
!       WRITE HEADER DESCRIBING DECOMPOSITION.
!
      Write (23, 140) varnam(deploc), varnam(poploc), pop1d, pop2d
140   Format (/, 1X, ' Decomposition of ', A10, /, 1X, ' comparing ', A10, &
        'categories', I3, ' and', I3)
      Write (23, 150)
150   Format ('________________________________________', &
        '_______________________________________')
!
!       CALL SUBROUTINE DESTAB TO GENERATE DESCRIPTIVE TABLE
!
      Call destab(vardes, varnam, facloc, numfac, n1t, n2t, mi1, mi2, pop1d, &
        pop2d, deploc, nu1, nu2, tot1, tot2, pg, ln)
!
!       IF THE "WRITE EXCLUDED" SUBCOMMAND WAS ISSUED AND SOME CELLS
!       WERE ACTUALLY EXCLUDED, THIS IS THE PLACE TO WRITE THE
!       CODEBOOK FOR THE EXCLUDED DATA TO THE LISTING FILE.
!       THIS TASK IS HANDLED BY SUBROUTINE EXCODE.
!
      exclu = mi1 + mi2
      If (wrex>0 .And. exclu>0) Then
        Call excode(varnam, wrform, files, poploc, facloc, numfac, pg, ln)
      End If
!
!       NOW WE ARE READY TO WRITE THE DECOMPOSITION RESULTS.
!       FIRST CALCULATE WETHER OR NOT WE WILL BE ABLE TO SQUEEZE
!       THEM ONTO THE CURRENT PAGE; IF NOT, MAKE A NEW PAGE.
!
      space = 20
      If (wrex>0) space = 21
      room = 63 - ln
      If (space>room) ln = 100
      Call page(pg, ln)
      If (ln>1) Then
        Write (23, 170)
        If (wrex>0) Write (23, 150)
      End If
!
!       WRITE THE HEADER FOR THE DECOMP TABLE.
!
      Write (23, 160)
160   Format (/, ' B. Decomposition Table')
170   Format (/)
      Write (23, 180)
180   Format (/, 8X, '                          Components of        Index of' &
        , /, 8X, '                           Difference         Difference', &
        /, 8X, '                          _____________       __________', /)
!
!       AS WE WRITE THE RESULTS, WE WILL CALCULATE THE INDEX OF
!       DIFFERENCE FOR EACH ITEM.
!       FIRST WRITE DIFFERENCE BETWEEN POPULATIONS.  THE INDEX
!       OF THIS IS 100.0 BY DEFINITION.
!
      ix = 100.0
      Write (23, 190) dif, ix
190   Format (' Total population difference   ', 3X, F10.4, 12X, F5.1, /, /)
!
!       WRITE EFFECTS OF EACH FACTOR
!
      Do fact = 1, numfac
        effect(fact) = effect(fact)/24.0
        If (dif>0) ix = effect(fact)/dif*100.0
        Write (23, 200) varnam(facloc(fact)), effect(fact), ix
      End Do
200   Format ('    Effect of factor ', A10, 3X, F10.4, 11X, F6.1)
!
!       WRITE COMBINED EFFECT OF FACTORS
!
      If (dif>0) ix = comb/dif*100.0
      Write (23, 210) comb, ix
210   Format (/, /, ' Combined effect of factors  ', 5X, F10.4, 11X, F6.1)
!
!       WRITE RATE EFFECT
!
      If (dif>0) ix = rate/dif*100.0
      Write (23, 220) rate, ix
220   Format (/, ' Rate effect                 ', 5X, F10.4, 11X, F6.1, /)
      Write (23, 150)
!
!       SET LINE NUMBER HIGH ENOUGH THAT A NEW PAGE WILL BE GENERATED
!       UNLESS THIS IS THE END OF THE JOB.
!
      ln = 100
      Return
    End Subroutine


    Subroutine destab(vardes, varnam, facloc, numfac, n1t, n2t, mi1, mi2, &
      pop1d, pop2d, deploc, nu1, nu2, tot1, tot2, pg, ln)
!
!       THIS SUBROUTINE CALCULATES AND WRITES A GENERAL DESCRIPTIVE
!       TABLE WHICH IS USEFUL FOR INTERPRETING THE DECOMPOSITION RESULTS.
!
      Implicit Integer (A-Z)
      Real mean1, mean2, pct1, pct2
      Real table, n1t, n2t, tot1, tot2
      Character *10 varnam
      Dimension table(640000, 2)
      Dimension maxvars(5), facloc(5)
      Dimension varnam(30)
      Dimension vardes(30, 5)
      Common /tabs/table, maxvars


      Write (23, 100)
100   Format (/, ' A. Descriptive Table')
!
!       BEFORE WRITING THE DECOMPOSITION TABLE, WE WILL CALCULATE
!       AND WRITE A SIMPLE DESCRIPTIVE TABLE GIVING THE MEAN OF THE
!       DEPENDENT FOR EACH CATEGORY OF EACH FACTOR AND THE FREQUENCY
!       DISTRIBUTION OF EACH FACTOR IN EACH POPULATION.  SUCH A TABLE
!       IS USUALLY CRITICAL TO THE INTERPRETATION OF THE DECOMPOSITION.
!
!       FIRST WRITE THE HEADER TO THE DESCRIPTIVE SECTION
!
      Write (23, 110) varnam(deploc), pop1d, pop2d, pop1d, pop2d
110   Format (/, 23X, 'Mean of ', A10, 11X, 'Population distribution', /, 19X, &
        '_______________________', 10X, '_______________________', /, /, 21X, &
        I3, 9X, I3, 19X, I3, 9X, I3, /, ' Factor/category', &
        '   ________________________________________________________')
!
!       WRITE THE NAME OF THE FIRST FACTOR.
!
      ln = ln + 3
      Write (23, 120) varnam(facloc(1))
120   Format (/, 4X, A10)
!
!       MINBRK IS THE DIFFERENCE BETWEEN THE FACTOR'S INDEX AND
!       THE ACTUAL VALUE OF THE VARIABLE
!
      minbrk = vardes(facloc(1), 3) - 1
!
!       LOOP GOES THROUGH EACH VALUE OF FIRST FACTOR
!
      Do i = 2, maxvars(2)
!
!       INITIALIZE MEANS AND PERCENTS
!
        mean1 = 0.0
        mean2 = 0.0
        pct1 = 0.0
        pct2 = 0.0
!
!       CALCULATE MEANS AND PERCENTS
!
        index = ind(1, i, 1, 1, 1)
        If (table(index,2)>0.0) Then
          mean1 = table(index, 1)/table(index, 2)
        End If
        If (n1t>0.0) pct1 = table(index, 2)/n1t*100.0
        index = ind(2, i, 1, 1, 1)
        If (table(index,2)>0.0) Then
          mean2 = table(index, 1)/table(index, 2)
        End If
        If (n2t>0.0) pct2 = table(index, 2)/n2t*100.0
!
!       SEE IF WE ARE AT THE END OF THE PAGE.
!
        Call page(pg, ln)
!
!       IF THIS IS A NEW PAGE, WRITE A NEW HEADER.
!
        If (ln==1) Then
          Write (23, 130)
130       Format (/, ' A. Descriptive Table (Continued)')
          Write (23, 110) varnam(deploc), pop1d, pop2d, pop1d, pop2d
          Write (23, 120) varnam(facloc(1))
          ln = ln + 13
        End If
!
!       NOW WRITE THE FACTOR CATEGORY AND THE MEANS AND PERCENTS.
!
        Write (23, 140) i - 1 + minbrk, mean1, mean2, pct1, pct2
140     Format (6X, I3, 8X, F9.4, 3X, F9.4, 16X, F5.1, 7X, F5.1)
      End Do
!
!       THE FOLLOWING SECTION DOES THE SAME THING FOR THE SECOND
!       THROUGH THE FOURTH FACTORS. CHECK FIRST TO MAKE SURE THAT
!       THESE FACTORS EXIST.
!
      If (numfac>=2) Then
!
!       DO THE LOOP FOR EACH FACTOR.
!
        Do nf = 2, numfac
!
!       IF THERE IS NOT GOING TO BE ROOM FOR THE ENTIRE FACTOR ON THIS
!       PAGE, WE WANT TO START A NEW PAGE (UNLESS THE FACTOR ISN'T GOING
!       TO FIT ON A SINGLE PAGE ANYWAY). THE FOLLOWING STATEMENTS
!       CALCULATE THE AMOUNT OF SPACE THE FACTOR WILL REQUIRE AND THE
!       AMOUNT OF ROOM LEFT ON THE PAGE, AND COMPARE THE TWO; IF THERE
!       ISN'T ENOUGH ROOM, WE CALL A NEW PAGE.
!
          trail = 0
          If (numfac==nf) trail = 5
          space = maxvars(nf+1) + trail + 2
          room = 63 - ln
          If (space>room .And. space<50) ln = 100
          Call page(pg, ln)
!
!         IF WE DID START A NEW PAGE, WE NEED TO WRITE A NEW HEADER.
!
          If (ln==1) Then
            Write (23, 130)
            Write (23, 110) varnam(deploc), pop1d, pop2d, pop1d, pop2d
            ln = ln + 9
          End If
          ln = ln + 2
!
!         WRITE THE NAME OF THE FACTOR
!
          Write (23, 120) varnam(facloc(nf))
!
!         CALCULATE MINBRK FOR THIS FACTOR
!
          minbrk = vardes(facloc(nf), 3) - 1
!
!         PASS THROUGH EACH CATEGORY OF THE CURRENT FACTOR
!
          Do i = 2, maxvars(nf+1)
!
!         INITIALIZE MEANS AND PERCENTS
!
            mean1 = 0.0
            mean2 = 0.0
            pct1 = 0.0
            pct2 = 0.0
!
!         CALCULATE MEANS AND PERCENTS, FIRST SETTING THE INDEX
!         TO THE APPROPRIATE MARGINALS FOR THIS FACTOR.
!
            If (nf==2) Then
              index = ind(1, 1, i, 1, 1)
            Else If (nf==3) Then
              index = ind(1, 1, 1, i, 1)
            Else If (nf==4) Then
              index = ind(1, 1, 1, 1, i)
            End If
            If (table(index,2)>0.0) Then
              mean1 = table(index, 1)/table(index, 2)
            End If
            If (n1t>0.0) pct1 = table(index, 2)/n1t*100.0
!
!         DITTO FOR POPULATION B
!
            If (nf==2) Then
              index = ind(2, 1, i, 1, 1)
            Else If (nf==3) Then
              index = ind(2, 1, 1, i, 1)
            Else If (nf==4) Then
              index = ind(2, 1, 1, 1, i)
            End If
            If (table(index,2)>0.0) Then
              mean2 = table(index, 1)/table(index, 2)
            End If
            If (n2t>0.0) pct2 = table(index, 2)/n2t*100.0
!
!         CHECK OUT PAGE POSITION AND WRITE NEW HEADER IF
!         THIS IS A NEW PAGE
!
            Call page(pg, ln)
            If (ln==1) Then
              Write (23, 130)
              Write (23, 110) varnam(deploc), pop1d, pop2d, pop1d, pop2d
              Write (23, 120) varnam(facloc(nf))
              ln = ln + 13
            End If
!
!         FINALLY, WRITE THE MEANS AND PERCENTS.
!
            Write (23, 140) i - 1 + minbrk, mean1, mean2, pct1, pct2
          End Do
        End Do
      End If
!
!       NOW WRITE TOTALS AND N'S AT THE END OF THE DESCRIPTIVE TABLE.
!
      pct1 = 100.0
      ln = ln + 4
      Call page(pg, ln)
      Write (23, 150) tot1, tot2, pct1, pct1
150   Format (/, '    Total', 8X, F9.4, 3X, F9.4, 16X, F5.1, 7X, F5.1)
      Write (23, 160) nu1, nu2, nu1, nu2
160   Format ('              N', 2X, I9, 3X, I9, 12X, I9, 3X, I9)
      Write (23, 170) mi1, mi2, mi1, mi2
170   Format (' Excluded cases', 2X, I9, 3X, I9, 12X, I9, 3X, I9)
      Write (23, 180)
180   Format ('________________________________________', &
        '_______________________________________')
      Return
    End Subroutine


    Subroutine margins(a, b, c, d, e, tden, tnum)
      Implicit Integer (A-Z)
      Real table, tden, tnum
      Dimension table(640000, 2)
      Dimension maxvars(5)
      Common /tabs/table, maxvars
!
!       THIS SUBROUTINE TABULATES EVERY POSSIBLE TYPE
!       OF MARGINAL.  FOR THE THREE-WAY AND TWO-WAY MARGINALS,
!       WE WILL ONLY NEED THE NUMBER OF CASES.
!
      index = ind(a, b, c, d, 1)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, b, c, 1, e)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, b, 1, d, e)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, 1, c, d, e)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, b, c, 1, 1)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, b, 1, d, 1)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, 1, c, d, 1)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, b, 1, 1, e)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, 1, c, 1, e)
      table(index, 2) = table(index, 2) + tden
      index = ind(a, 1, 1, d, e)
      table(index, 2) = table(index, 2) + tden
!
!       FOR THE ONE-WAY MARGINALS, WE NEED THE NUMERATORS AS
!       WELL, BECAUSE WE WILL NEED THEM FOR THE DESCRIPTIVE TABLE.
!
      index = ind(a, b, 1, 1, 1)
      table(index, 2) = table(index, 2) + tden
      table(index, 1) = table(index, 1) + tnum
      index = ind(a, 1, 1, d, 1)
      table(index, 2) = table(index, 2) + tden
      table(index, 1) = table(index, 1) + tnum
      index = ind(a, 1, c, 1, 1)
      table(index, 2) = table(index, 2) + tden
      table(index, 1) = table(index, 1) + tnum
      index = ind(a, 1, 1, 1, e)
      table(index, 2) = table(index, 2) + tden
      table(index, 1) = table(index, 1) + tnum
!
!       FINALLY, THE NO-WAY MARGINAL (THE TOTAL CASES IN EACH POPULATION)
!
      index = ind(a, 1, 1, 1, 1)
      table(index, 2) = table(index, 2) + tden
      Return
    End Subroutine


    Subroutine adjust(a, i, j, k, l, n, t)
!
!         THIS ROUTINE ADJUSTS THE MARGINALS TO ACCOUNT FOR EXCLUDED
!         CASES.
!
      Implicit Integer (A-Z)
      Real table, n, t
      Dimension table(640000, 2)
      Dimension maxvars(5)
      Common /tabs/table, maxvars

      table(ind(a,i,j,k,1), 2) = table(ind(a,i,j,k,1), 2) - n
      table(ind(a,i,j,1,l), 2) = table(ind(a,i,j,1,l), 2) - n
      table(ind(a,i,1,k,l), 2) = table(ind(a,i,1,k,l), 2) - n
      table(ind(a,1,j,k,l), 2) = table(ind(a,1,j,k,l), 2) - n
      table(ind(a,1,1,k,l), 2) = table(ind(a,1,1,k,l), 2) - n
      table(ind(a,1,j,1,l), 2) = table(ind(a,1,j,1,l), 2) - n
      table(ind(a,1,j,k,1), 2) = table(ind(a,1,j,k,1), 2) - n
      table(ind(a,i,1,1,l), 2) = table(ind(a,i,1,1,l), 2) - n
      table(ind(a,i,j,1,1), 2) = table(ind(a,i,j,1,1), 2) - n
      table(ind(a,i,1,k,1), 2) = table(ind(a,i,1,k,1), 2) - n
      index = ind(a, i, 1, 1, 1)
      table(index, 2) = table(index, 2) - n
      table(index, 1) = table(index, 1) - t
      index = ind(a, 1, j, 1, 1)
      table(index, 2) = table(index, 2) - n
      table(index, 1) = table(index, 1) - t
      index = ind(a, 1, 1, k, 1)
      table(index, 2) = table(index, 2) - n
      table(index, 1) = table(index, 1) - t
      index = ind(a, 1, 1, 1, l)
      table(index, 2) = table(index, 2) - n
      table(index, 1) = table(index, 1) - t
      Return
    End Subroutine

    Subroutine gupset(vardes, pop1d, pop2d, poploc, numfac, newpos, facloc, &
      poppos, popa, popb)
!
!      THIS SUBROUTINE DOES A LITTLE PRELIMINARY SETUP NECESSARY
!      FOR THE DECOMP PROCEDURE.
!
      Implicit Integer (A-Z)
      Dimension maxvars(5), facloc(5)
      Dimension vardes(30, 5)
      Dimension newpos(5)
      Dimension table(640000, 2)
      Common /tabs/table, maxvars
!
!      CHECK TO MAKE SURE THE COMPARISON VARIABLE HAS BEEN
!      SPECIFIED IN THE PRECEEDING MAKETAB COMMAND.
!
      poppos = vardes(poploc, 5) - 1
      If (poppos<1 .Or. poppos>5) Then
        Call error(402)
      End If
!
!      INITIALIZE MAXVARS TO TWO.  THIS LEAVES ROOM FOR THE AVERAGE OF
!      THE POPULATIONS IN POSITION THREE (WE WILL NEED TO HAVE THE
!      AVERAGE LATER).
!
!      THE ARRAY NEWPOS WILL CONTAIN THE NEW POSITION OF EACH OF
!      THE VARIABLES PREVIOUSLY TABULATED BY MAKETAB. INITIALIZE.
!
      Do i = 1, 5
        maxvars(i) = 2
        newpos(i) = 6
      End Do
!
!      THE COMPARISON VARIBLE, LOCATED IN POPPOS, WILL BE STUCK IN
!      THE FIRST POSITION IN THE NEW TABLE.
!
      newpos(1) = poppos
!
!      FIGURE THE RANGE OF THE COMPARISON VARIABLE, AND PUT IT IN
!      MAXVARS(1)
!
      maxvars(1) = vardes(poploc, 4) - vardes(poploc, 3) + 1
!
!      MAKE SURE THE SPECIFIC COMPARISON CATEGORIES FALL WITHIN THE
!      SPECIFIED RANGE OF THE COMPARISON VARIABLE.
!
      popa = pop1d - vardes(poploc, 3) + 1
      popb = pop2d - vardes(poploc, 3) + 1
      If (popa>maxvars(1) .Or. popb>maxvars(1)) Then
        Call error(405)
      End If
!
!      THE FOLLOWING LOOP DOES THE SAME PROCEDURES FOR EACH OF
!      THE FACTORS.
!
      next = 1
      ndenow = numfac
      Do i = 1, numfac
!
!      IF THE COMPARISON VARIABLE WAS SPECIFIED AS A FACTOR,
!      WE WILL SIMPLY THROW IT OUT AS A FACTOR.
!      MAKE SURE, IF IT IS THROWN OUT, THAT WE STILL HAVE AT
!      LEAST ONE FACTOR TO WORK WITH.
!
        If (facloc(i)==poploc) Then
          ndenow = ndenow - 1
          If (ndenow==0) Then
            Call error(402)
          End If
          Go To 100
        End If
!
!      THE POSITION OF EACH OF THE FACTORS ON THE NEW TABLE IS
!      INDICATED BY NEXT.  THEIR PREVIOUS POSITION, IN THE TABLE
!      CREATED BY VARDES, CAN BE DETERMINED FROM COLUMN 5 OF VARDES.
!
        next = next + 1
        newpos(next) = vardes(facloc(i), 5) - 1
!
!      MAKE SURE WE HAVE THE RIGHT VARIABLES; IF NOT, WRITE AN ERROR.
!
        If (newpos(next)<1 .Or. newpos(next)>5) Then
          Call error(402)
        End If
!
!      CALCULATE MAXIMUM RANGE FOR EACH VARIABLE.  LEAVE ROOM FOR
!      THE MARGINALS.
!
        maxvars(next) = vardes(facloc(i), 4) - vardes(facloc(i), 3) + 2
100   End Do
!
!      NOW WE CAN RESET THE MAXVARS FOR THE COMPARISON VARIABLE;
!      WE KNOW THAT ONLY TWO CATEGORIES WILL BE COMPARED, SO WE
!      SET IT TO THREE, LEAVING ROOM FOR MARGINALS.
!
      maxvars(1) = 3
!
!      NOW MAKE SURE THAT THE NEW TABLE WILL FIT IN THE MAIN ARRAY.
!
      biggest = maxvars(1)*maxvars(2)*maxvars(3)*maxvars(4)*maxvars(5)
      If (biggest>640000) Then
        Call error(406)
      End If
      Return
    End Subroutine
!
!
!
    Subroutine wrprep(varnum, vardes, wrform, files, numind, indnam, wtdig, &
      pg, ln)
!
!       THIS ROUTINE DOES THE SETUP FOR THE "WRITE TABLE" SUBCOMMAND
!       OF THE MAKETAB COMMAND.  IT ALSO WRITES A CODEBOOK FOR THE
!       TABLE TO THE LISTING FILE.
!
      Implicit Integer (A-Z)
      Real table
      Character *10 indnam
      Character *35 files
      Dimension vardes(30, 5)
      Dimension files(4)
      Dimension indnam(5)
      Dimension wrform(7, 2)
      Common /tabs/table(640000, 2), maxvars(5)
!
!       OPEN THE OUTPUT FILE.  THE NAME OF THE FILE IS STORED IN
!       FILES(4)
!
      Open (66, File=files(4))
!
!       THE FOLLOWING LOOP FILLS AN ARRAY WRFORM.
!       THE FIRST SUBSCRIPT INDICATES THE OUTPUT VARIABLE NUMBER.
!       WRFORM (N,1) WILL CONTAIN THE STARTING COLUMN OF THE VARIABLE;
!       WRFORM (N,2) WILL CONTAIN THE ENDING COLUMN OF THE VARIABLE.
!
!       THE PROGRAM IS SET UP TO SKIP A SPACE BETWEEN EACH OUTPUT
!       VARIABLE, IN ORDER TO ALLOW USERS TO INPUT THE TABLE INTO
!       PROGRAMS REQUIRING FREE-FORM INPUT.
!
!       START AT CHARACTER -1, AND FIRST LOOK FOR THE INDEPENDENT
!       VARIABLES.  THEY ARE IDENTIFIED IN VARDES(N,5) WITH VALUES
!       BETWEEN 2 AND 6.
!
      chr = -1
      Do pos = 2, 6
        Do i = 1, varnum
!
!       IF THE CURRENT VARIABLE IS AN INDEPENDENT VARIABLE, THE
!       BEGINNING COLUMN IS SET TWO COLUMNS BEYOND THE LAST COLUMN
!       SPECIFIED (NOTE: TO ELIMINATE THE SPACE BETWEEN VARIABLES,
!       CHANGE THIS INCREMENT TO ONE)
!       SINCE CHR WAS INITIALIZED AT COLUMN -1, THE FIRST VARIABLE
!       WILL BEGIN IN COLUMN 1.
!
          If (vardes(i,5)==pos) Then
            wrform(pos-1, 1) = chr + 2
!
!       NOW DETERMINE THE ENDING COLUMN BY LOOKING AT THE MAXIMUM
!       SIZE OF THE VARIABLE.
!
            If (vardes(i,4)<10) Then
              wrform(pos-1, 2) = chr + 2
            Else If (vardes(i,4)<100) Then
              wrform(pos-1, 2) = chr + 3
            Else If (vardes(i,4)<1000) Then
              wrform(pos-1, 2) = chr + 4
            Else If (vardes(i,4)<10000) Then
              wrform(pos-1, 2) = chr + 5
            Else
              wrform(pos-1, 2) = chr + 6
            End If
            chr = wrform(pos-1, 2)
          End If
        End Do
      End Do
!
!       THAT TAKES CARE OF THE INDEPENDENT VARIABLES.
!       WRITE THE CODES TO THE LISTING FILE.
!       FIRST, SEE IF WE WILL NEED TO START A NEW PAGE AND
!       WRITE A HEADER.
!
      ln = ln + numind + 10
      Call page(pg, ln)
      Write (23, 100) files(4)
100   Format (/, ' The table is written to file ', A35, /, &
        ' using the following format:', /, /, 22X, 'Variable', /, 22X, &
        'Name                   Columns', /)
!
!       NOW WRITE THE FIRST PART OF THE CODEBOOK.
!
      Do i = 1, numind
        Write (23, 110) indnam(i), wrform(i, 1), wrform(i, 2)
      End Do
110   Format (22X, A10, 12X, I2, '-', I2)
!
!       NOW SEARCH THE MAIN TABLE FOR THE LARGEST MEAN VALUE OF THE
!       DEPENDENT VARIABLE (BIGNUM) AND THE LARGEST NUMBER OF
!       CASES (BIGDEN).
!
      bignum = 0
      bigden = 0
      Do i = 1, 640000
        If (table(i,2)>0.0) Then
          If (table(i,1)/table(i,2)>float(bignum)) bignum = int(table(i,1)/ &
            table(i,2))
        End If
        If (table(i,2)>float(bigden)) bigden = int(table(i,2))
      End Do
!
!       THE MEAN OF THE DEPENDENT WILL BE WRITTEN BEGINNING TWO
!       COLUMNS AFTER THE LAST DEPENDENT VARIABLE.
!
      wrform(numind+1, 1) = wrform(numind, 2) + 2
!
!       THE END OF IT WILL DEPEND ON BIGNUM.  LEAVE FIVE ADDITIONAL
!       PLACES FOR A DECIMAL POINT AND DIGITS.
!
      If (bignum<10) Then
        wrform(numind+1, 2) = wrform(numind+1, 1) + 5
      Else If (bignum<100) Then
        wrform(numind+1, 2) = wrform(numind+1, 1) + 6
      Else If (bignum<1000) Then
        wrform(numind+1, 2) = wrform(numind+1, 1) + 7
      Else If (bignum<10000) Then
        wrform(numind+1, 2) = wrform(numind+1, 1) + 8
      Else
        wrform(numind+1, 2) = wrform(numind+1, 1) + 9
      End If
!
!       NOW FIGURE THE COLUMNS FOR THE NUMBER OF CASES IN THE CELL.
!       IN CASE THE WEIGHT VARIABLE HAS BEEN READ IN AS A REAL NUMBER,
!       WE WILL HAVE TO FIGURE HOW MANY PLACES TO THE RIGHT OF THE
!       DECIMAL ARE NEEDED.
!
      places = 0
      If (wtdig>0) places = wtdig + 1
      wrform(numind+2, 1) = wrform(numind+1, 2) + 2
      If (bigden<10) Then
        wrform(numind+2, 2) = wrform(numind+2, 1) + places
      Else If (bigden<100) Then
        wrform(numind+2, 2) = wrform(numind+2, 1) + places + 1
      Else If (bigden<1000) Then
        wrform(numind+2, 2) = wrform(numind+2, 1) + places + 2
      Else If (bigden<10000) Then
        wrform(numind+2, 2) = wrform(numind+2, 1) + places + 3
      Else If (bigden<100000) Then
        wrform(numind+2, 2) = wrform(numind+2, 1) + places + 4
      Else
        wrform(numind+2, 2) = wrform(numind+2, 1) + places + 5
      End If
!
!       NOW WRITE THE REST OF THE TABLE: THE COLUMNS FOR THE MEAN OF THE
!       DEPENDENT AND NUMBER OF CASES IN THE CELL.
!
      Write (23, 120) wrform(numind+1, 1), wrform(numind+1, 2), &
        wrform(numind+2, 1), wrform(numind+2, 2)
120   Format (22X, 'Mean of dependent', 5X, I2, '-', I2, /, 22X, &
        'Number of cases  ', 5X, I2, '-', I2, /)
      Return
    End Subroutine


    Subroutine wrdata(varnum, vardes, wrform, i, j, k, l, m, numind, wtdig, &
      wtdiv, tnum, tden)
!
!       THIS SUBROUTINE ACTUALLY WRITES THE TABLE WHEN THE "WRITE TABLE"
!       OPTION OF MAKETAB IS SPECIFIED.
!
      Implicit Integer (A-Z)
      Integer i, j, k, l, m
      Character *10 convert
      Character *80 wrline
      Real tnum, tden
      Dimension vardes(30, 5)
      Dimension wrform(7, 2)
!
!       THE DATA WILL BE CONVERTED TO A CHARACTER STRING BEFORE IT IS
!       WRITTEN, BOTH TO SIMPLIFY THE WRITE FORMAT AND TO SPEED THE
!       PROCESS.
!
! ****** WARNING: SOME COMPILERS DO NOT ALLOW CHARACTER STRINGS THIS LONG.
!
!       INITIALIZE THE LINE WITH BLANKS.
!
      wrline(1:40) = '                                        '
      wrline(41:80) = '                                        '
!
!       THE FOLLOWING LOOP INSERTS THE DATA INTO THE STRING WRLINE.
!       THE POSITION OF EACH VARIABLE HAS BEEN PREVIOUSLY DETERMINED
!       BY THE SUBROUTINE WRPREP, AND STORED IN ARRAY WRFORM.
!       THE VALUES OF EACH VARIABLE ARE CONVERTED TO CHARACTERS BY
!       THE FUNCTION "CONVERT".
!       (THE VALUES OF THE VARIABLES HAVE BEEN TRANSFERRED TO THIS
!       SUBROUTINE AS I,J,K,L, AND M.
!
      Do f = 1, varnum
        pos = vardes(f, 5) - 1
        If (pos==1) Then
          wrline(wrform(pos,1):wrform(pos,2)) = convert(i+vardes(f,3)-1, &
            wrform(pos,2)-wrform(pos,1)+1)
        Else If (pos==2) Then
          wrline(wrform(pos,1):wrform(pos,2)) = convert(j+vardes(f,3)-1, &
            wrform(pos,2)-wrform(pos,1)+1)
        Else If (pos==3) Then
          wrline(wrform(pos,1):wrform(pos,2)) = convert(k+vardes(f,3)-1, &
            wrform(pos,2)-wrform(pos,1)+1)
        Else If (pos==4) Then
          wrline(wrform(pos,1):wrform(pos,2)) = convert(l+vardes(f,3)-1, &
            wrform(pos,2)-wrform(pos,1)+1)
        Else If (pos==5) Then
          wrline(wrform(pos,1):wrform(pos,2)) = convert(m+vardes(f,3)-1, &
            wrform(pos,2)-wrform(pos,1)+1)
        End If
      End Do
!
!       NOW WE HAVE TO INSERT THE MEAN OF THE DEPENDENT VARIABLE.
!       TO DO THIS WE MUST SPLIT IT INT TWO PARTS.
!       DIGITS TO THE LEFT OF THE DECIMAL ARE PLACED IN THE
!       VARIABLE "INTEG."  DIGITS TO THE RIGHT OF THE DECIMAL
!       ARE PLACED IN THE VARIABLE "FRACT."
!
      integ = int((tnum/tden))
! INT4 WAS HERE
      fract = int(((tnum/tden+.00005)-(tnum/tden))*10000.0)
      If (fract==10000) Then
        fract = 0
        integ = integ + 1
      End If
!
!       NOW INSERT THE INTEG PART, FIRST CONVERTING IT WITH "CONVERT"
!
      beg = wrform(numind+1, 1)
      end = wrform(numind+1, 2) - 5
      wrline(beg:end) = convert(integ, end-beg+1)
!
!       NOW INSERT THE DECIMAL POINT IN THE APPROPRIATE SPOT.
!
      wrline(end+1:end+1) = '.'
!
!       AND FINALLY INSERT THE FRACT PART.
!
      beg = end + 2
      end = wrform(numind+1, 2)
      wrline(beg:end) = convert(fract, end-beg+1)
!
!       THE LAST ITEM TO BE WRITTEN IS THE NUMBER OF CASES IN THE CELL.
!       IF IT IS A REAL NUMBER, WE MUST AGAIN SPLIT IT INTO INTEG AND
!       FRACT, BUT IF NOT (INDICATED BY WTDIG=0) WE CAN SIMPLY CONVERT
!       IT AS AN INTEGER.
!
      If (wtdig==0) Then
        integ = int(tden)
        beg = wrform(numind+2, 1)
        end = wrform(numind+2, 2)
!
!       INSERT INTEGER VALUE
!
        wrline(beg:end) = convert(integ, end-beg+1)
      Else
!
!       OR SPLIT IT UP, AND INSERT INTEG, DECIMAL POINT, AND FRACT.
!
        integ = int(tden)
        beg = wrform(numind+2, 1)
        end = wrform(numind+2, 2) - wtdig - 1
        wrline(end+1:end+1) = '.'
        wrline(beg:end) = convert(integ, end-beg+1)
! INT4 WAS HERE
        fract = int((tden-(tden))*wtdiv)
        beg = end + 2
        end = wrform(numind+2, 2)
        wrline(beg:end) = convert(fract, end-beg+1)
      End If
!
!       NOW THE STRING IS COMPLETE.  WRITE IT.
!
      Write (66, 100) wrline(1:wrform(numind+2,2))
100   Format (A)
      Return
    End Subroutine


    Subroutine exprep(vardes, wrform, files, deploc, wtdig, locate, numind, &
      nvars)
!
!       THIS ROUTINE DOES THE SETUP FOR THE "WRITE EXCLUDED" SUBCOMMAND
!       OF A PROCEDURE COMMAND.
!
      Implicit Integer (A-Z)
      Real table
      Character *35 files
      Dimension locate(5)
      Dimension vardes(30, 5)
      Dimension files(4)
      Dimension wrform(7, 2)
      Common /tabs/table(640000, 2), maxvars(5)
!
!       OPEN THE OUTPUT FILE.  THE NAME OF THE FILE IS STORED IN
!       FILES(4)
!
      Open (66, File=files(4))
!
!       THE FOLLOWING SECTION FILLS AN ARRAY WRFORM.
!       THE FIRST SUBSCRIPT INDICATES THE OUTPUT VARIABLE NUMBER.
!       WRFORM (N,1) WILL CONTAIN THE STARTING COLUMN OF THE VARIABLE;
!       WRFORM (N,2) WILL CONTAIN THE ENDING COLUMN OF THE VARIABLE.
!
!       THE PROGRAM IS SET UP TO SKIP A SPACE BETWEEN EACH OUTPUT
!       VARIABLE, IN ORDER TO ALLOW USERS TO INPUT THE TABLE INTO
!       PROGRAMS REQUIRING FREE-FORM INPUT.
!
!       FIRST LOOK FOR THE COMPARISON OR BREAKDOWN
!       VARIABLE, LOCATED IN POSITION DEPLOC OF VARDES.
!       THE NUMBER OF COLUMNS REQUIRED IS DETERMINED BY THE MAXIMUM
!       VALUE OF THE VARIABLE.
!

      pos = 1
      wrform(pos, 1) = 1
      If (vardes(deploc,4)<10) Then
        wrform(pos, 2) = 1
      Else If (vardes(deploc,4)<100) Then
        wrform(pos, 2) = 2
      Else If (vardes(deploc,4)<1000) Then
        wrform(pos, 2) = 3
      Else If (vardes(deploc,4)<10000) Then
        wrform(pos, 2) = 4
      Else
        wrform(pos, 2) = 5
      End If
      chr = wrform(pos, 2)
!
!       NOW DO THE INDEPENDENT VARIABLES; THEIR NUMBER IS NUMIND,
!       AND THEIR LOCATION IS GIVEN IN ARRAY LOCATE.
!
      Do i = 1, numind
        If (locate(i)/=deploc) Then
          pos = pos + 1
          wrform(pos, 1) = chr + 2
          If (vardes(locate(i),4)<10) Then
            wrform(pos, 2) = chr + 2
          Else If (vardes(locate(i),4)<100) Then
            wrform(pos, 2) = chr + 3
          Else If (vardes(locate(i),4)<1000) Then
            wrform(pos, 2) = chr + 4
          Else If (vardes(locate(i),4)<10000) Then
            wrform(pos, 2) = chr + 5
          Else
            wrform(pos, 2) = chr + 6
          End If
          chr = wrform(pos, 2)
        End If
      End Do
!
!       NOW DETERMINE THE POSITIONS OF THE MEAN VALUE OF THE DEPENDENT
!       VARIABLE AND THE NUMBER OF CASES IN THE CELL.
!       TO FIGURE OUT HOW MANY COLUMNS WILL BE REQUIRED, WE SEARCH THE
!       MAIN TABLE FOR THE MAXIMUM VALUES, WHICH ARE STUCK IN BIGNUM
!       AND BIGDEN.
!
      bignum = 0
      bigden = 0
      Do i = 1, 640000
        If (table(i,2)>0.0) Then
          If (table(i,1)/table(i,2)>float(bignum)) bignum = int(table(i,1)/ &
            table(i,2))
        End If
        If (table(i,2)>float(bigden)) bigden = int(table(i,2))
      End Do
!
!       FIGURE OUT COLUMNS FOR MEAN
!
      wrform(pos+1, 1) = wrform(pos, 2) + 2
      If (bignum<10) Then
        wrform(pos+1, 2) = wrform(pos+1, 1) + 5
      Else If (bignum<100) Then
        wrform(pos+1, 2) = wrform(pos+1, 1) + 6
      Else If (bignum<1000) Then
        wrform(pos+1, 2) = wrform(pos+1, 1) + 7
      Else If (bignum<10000) Then
        wrform(pos+1, 2) = wrform(pos+1, 1) + 8
      Else
        wrform(pos+1, 2) = wrform(pos+1, 1) + 9
      End If
!
!       FIGURE OUT COLMNS FOR NUMBER OF CASES, LEAVING ROOM FOR
!       DECIMAL PLACES IF THE CASES ARE EXPRESSED AS REAL NUMBERS.
!
      places = 0
      If (wtdig>0) places = wtdig + 1
      wrform(pos+2, 1) = wrform(pos+1, 2) + 2
      If (bigden<10) Then
        wrform(pos+2, 2) = wrform(pos+2, 1) + places
      Else If (bigden<100) Then
        wrform(pos+2, 2) = wrform(pos+2, 1) + places + 1
      Else If (bigden<1000) Then
        wrform(pos+2, 2) = wrform(pos+2, 1) + places + 2
      Else If (bigden<10000) Then
        wrform(pos+2, 2) = wrform(pos+2, 1) + places + 3
      Else If (bigden<100000) Then
        wrform(pos+2, 2) = wrform(pos+2, 1) + places + 4
      Else
        wrform(pos+2, 2) = wrform(pos+2, 1) + places + 5
      End If
      nvars = pos
      Return
    End Subroutine


    Subroutine excode(varnam, wrform, files, deploc, locate, numind, pg, ln)
!
!       THIS SUBROUTINE WRITES A CODEBOOK TO THE LISTING FILE
!       WHEN A "WRITE EXCLUDED" SUBCOMMAND IS ISSUED.
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *35 files
      Dimension locate(5)
      Dimension files(4)
      Dimension varnam(30)
      Dimension wrform(7, 2)
!
!       FIRST SEE IF THERE WILL BE ENOUGH ROOM ON THE PAGE,
!       AND START A NEW PAGE IF NECESSARY.
!
      ln = ln + numind + 10
      Call page(pg, ln)
!
!       WRITE A HEADER.
!
      Write (23, 100) files(4)
100   Format (/, ' The excluded data are written to file ', A35, /, &
        ' using the following format:', /, /, 22X, 'Variable', /, 22X, &
        'Name                   Columns', /)
!
!       THE COUMNS FOR EACH VARIABLE ARE STORED IN ARRAY WRFORM.
!       WRITE THEM.
!
      Write (23, 110) varnam(deploc), wrform(1, 1), wrform(1, 2)
      pos = 1
      Do i = 1, numind
        If (locate(i)/=deploc) Then
          pos = pos + 1
          Write (23, 110) varnam(locate(i)), wrform(pos, 1), wrform(pos, 2)
        End If
      End Do
110   Format (22X, A10, 12X, I2, '-', I2)
!
!       THE MEAN OF DEPENDENT AND NUMBER OF CASES COME LAST.
!
      Write (23, 120) wrform(pos+1, 1), wrform(pos+1, 2), wrform(pos+2, 1), &
        wrform(pos+2, 2)
120   Format (22X, 'Mean of dependent', 5X, I2, '-', I2, /, 22X, &
        'Number of cases  ', 5X, I2, '-', I2, /)
      Return
    End Subroutine



    Subroutine exdata(wrform, i, j, k, l, m, nvars, wtdig, wtdiv, tnum, tden)
!
!       THIS SUBROUTINE WRITES THE EXCLUDED CASES WHEN THE
!       "WRITE EXCLUDED" SUBCOMMAND HAS BEEN SPECIFIED.
!
      Implicit Integer (A-Z)
      Character *10 convert
      Character *80 wrline
      Real tnum, tden
      Dimension wrform(7, 2)
      Dimension indv(5)
!
!       THE DATA WILL BE CONVERTED TO A CHARACTER STRING BEFORE IT IS
!       WRITTEN, BOTH TO SIMPLIFY THE WRITE FORMAT AND TO SPEED THE
!       PROCESS.
!
! ****** WARNING: SOME COMPILERS DO NOT ALLOW CHARACTER STRINGS THIS LONG.
!
!       INITIALIZE THE LINE WITH BLANKS.
!
      wrline(1:40) = '                                        '
      wrline(41:80) = '                                        '
!
!       THE FOLLOWING LOOP INSERTS THE DATA INTO THE STRING WRLINE.
!       THE POSITION OF EACH VARIABLE HAS BEEN PREVIOUSLY DETERMINED
!       BY THE SUBROUTINE EXPREP, AND STORED IN ARRAY WRFORM.
!       THE VALUES OF EACH VARIABLE ARE CONVERTED TO CHARACTERS BY
!       THE FUNCTION "CONVERT".
!       (THE VALUES OF THE VARIABLES HAVE BEEN TRANSFERRED TO THIS
!       SUBROUTINE AS I,J,K,L, AND M. FIRST STICK THEM IN THE ARRAY
!       INDV TO SIMPLFY PROCESSING).
!
      indv(1) = i
      indv(2) = j
      indv(3) = k
      indv(4) = l
      indv(5) = m
      Do n = 1, nvars
        wrline(wrform(n,1):wrform(n,2)) = convert(indv(n), &
          wrform(n,2)-wrform(n,1)+1)
      End Do
!
!       NOW WE HAVE TO INSERT THE MEAN OF THE DEPENDENT VARIABLE.
!       TO DO THIS WE MUST SPLIT IT INT TWO PARTS.
!       DIGITS TO THE LEFT OF THE DECIMAL ARE PLACED IN THE
!       VARIABLE "INTEG."  DIGITS TO THE RIGHT OF THE DECIMAL
!       ARE PLACED IN THE VARIABLE "FRACT."
!
      pos = nvars
      integ = int(tnum/tden)
! INT4 WAS HERE
      fract = int(((tnum/tden+.00005)-(tnum/tden))*10000.0)
      If (fract==10000) Then
        fract = 0
        integ = integ + 1
      End If
!
!       NOW INSERT THE INTEG PART.
!
      beg = wrform(pos+1, 1)
      end = wrform(pos+1, 2) - 5
      wrline(beg:end) = convert(integ, end-beg+1)
!
!       NOW STICK IN THE DECIMAL POINT.
!
      wrline(end+1:end+1) = '.'
!
!       NOW INSERT THE FRACT PART.
!
      beg = end + 2
      end = wrform(pos+1, 2)
      wrline(beg:end) = convert(fract, end-beg+1)
!
!       THE LAST ITEM TO BE WRITTEN IS THE NUMBER OF CASES IN THE CELL.
!       IF IT IS A REAL NUMBER, WE MUST AGAIN SPLIT IT INTO INTEG AND
!       FRACT, BUT IF NOT (INDICATED BY WTDIG=0) WE CAN SIMPLY CONVERT
!       IT AS AN INTEGER.
!
      If (wtdig==0) Then
        integ = int(tden)
        beg = wrform(pos+2, 1)
        end = wrform(pos+2, 2)
        wrline(beg:end) = convert(integ, end-beg+1)
      Else
!
!       OTHERWISE, WE DO THE INTEG AND FRACT BUSINESS.
!       THE NUMBER OF PLACES TO THE RIGHT OF THE DECIMAL POINT
!       IS GIVEN IN WTDIG.
!
        integ = int(tden)
        beg = wrform(pos+2, 1)
        end = wrform(pos+2, 2) - wtdig - 1
        wrline(end+1:end+1) = '.'
        wrline(beg:end) = convert(integ, end-beg+1)
! INT4 WAS HERE
        fract = int((tden-(tden))*wtdiv)
        beg = end + 2
        end = wrform(pos+2, 2)
        wrline(beg:end) = convert(fract, end-beg+1)
      End If
!
!      THE STRING IS COMPLETE, AND CAN BE WRITTEN.
!
      Write (66, 100) wrline(1:wrform(pos+2,2))
100   Format (A)
      Return
    End Subroutine





    Character *10 Function convert(integ, len)
!
!       THIS IS THE FUNCTION TO CONVERT INTEGERS TO CHARACTERS.
!
! ****** WARNING: THIS FUNCTION ASSUMES AN EXTENDED ASCII CHARACTER SET.
!                 IT WILL HAVE TO BE MODIFIED FOR COMPUTERS USING OTHER
!                 CHARACTER CODES.
!
      Implicit Integer (A-Z)

      int = integ
      Do i = 1, len
        int2 = (int/10)*10
        digit = int - int2
        place = len - i + 1
        convert(place:place) = char(digit+48)
        int = int/10
      End Do
      Return
    End Function

    Subroutine page(pg, ln)
!
!      THIS SUBROUTINE KEEPS TRACK OF LINE AND PAGE NUMBERS.
!      AT THE BEGINNING OF EACH PAGE, A HEADER IS WRITTEN
!      CONTAINING THE PAGE NUMBER AND DATE.
!
      Implicit Integer (A-Z)
!cc       INTEGER*4 YR,MON,DAY
        character *19 datetime
        character*8 date
        character*10 time

      Save skip
!
!      FIRST CHECK TO SEE IF PAGE NUMBERING HAS BEEN SUPPRESSED
!      SUPRESSION IS INDICATED IF PG IS EVER EQUAL TO 9999
!      SAVE THIS IN VARIABLE SKIP
!
      If (pg==9999) skip = pg
!
      If (skip/=9999) Then
!
!      LN IS THE LINE NUMBER.  IF IT IS OVER PGSIZE, START A NEW PAGE.
!
        ln = ln + 1
        If (ln>=61) Then
          ln = 1
          pg = pg + 1
!
!  ****** WARNING: THE FOLLOWING LIBRARY ROUTINE, WHICH GETS THE DATE
!                  WILL NOT WORK ON MOST SYSTEMS.
!cc note function changed for Sun
!cc changed again for GNU
!cc changed again for F90
!             I=TIME8()
!              CALL ctime(I,DATE)
!
!cc           CALL IDATE(DAY,MON,YR)
!
!  ****** WARNING: THE FOLLOWING STATEMENT WRITES THE FORMFEED CHARACTER.
!                  THIS WILL NOT WORK ON ALL SYSTEMS.
!
          If (pg>1) Write (23, *) char(012)
!
!        NOW WRITE THE PAGE HEADER.
!        DEPENDING ON THE DAY OF THE MONTH, THE WRITE FORMAT CAN
!        VARY SLIGHTLY.
!
!           IF (DAY.LT.10)THEN
!             WRITE(23,100)PG,MON,DAY,YR
!           ELSE
!             WRITE(23,101)PG,MON,DAY,YR
!           ENDIF
!            WRITE(23,100)PG,DATE
        call date_and_time(date, time)
        datetime(1:2)=date(5:6)
        datetime(3:3)='/'
        datetime(4:5)=date(7:8)
        datetime(6:6)='/'
        datetime(7:10)=date(1:4)
        datetime(11:11)=' '
        datetime(12:13)=time(1:2)
        datetime(14:14)=':'  
        datetime(15:16)=time(3:4)
        datetime(17:17)=':'  
        datetime(18:19)=time(5:6)

          Write (23, 100) pg, datetime

!        WRITE A LINE ACROSS THE TOP OF THE PAGE.
!
          Write (23, 110)
        End If
      End If
100   Format ('Page', I3, 15X, 'DECOMP Version 2.0', 7X, A31)

! 100   FORMAT('Page',I3,23X,'DECOMP Version 0.51',21x,
!     x        i2,'/',i1,'/',i4)
! 101   FORMAT('Page',I3,23X,'DECOMP Version 0.51',20x,
!     x        i2,'/',i2,'/',i4)
110   Format ('________________________________________', &
        '_______________________________________')
      Return
    End Subroutine
!
!
!
    Subroutine error(errnum)
!
!         THIS ROUTINE WRITES ORDINARY ERROR MESSAGES.
!         THE MESSAGES ARE STUCK IN ARRAYS AND THEN REFERRED TO
!         BY THEIR NUMBER.
!
!         THERE ARE EIGHT CHARACTER ARRAYS TO HOLD THE MESSAGES;
!         SEVEN OF THESE REFER TO ERRORS IN SPECIFIC COMMANDS,
!         AND THE OTHER ONE HOLDS MISCELLANEOUS MESSAGES.
!
      Implicit Integer (A-Z)
      Character *35 daterr
      Character *35 makerr
      Character *55 sterr
      Character *55 decerr
      Character *35 recerr
      Character *35 selerr
      Character *35 generr
      Character *35 comerr
      Dimension daterr(8)
      Dimension makerr(11)
      Dimension sterr(15)
      Dimension decerr(10)
      Dimension recerr(6)
      Dimension selerr(3)
      Dimension generr(6)
      Dimension comerr(6)
      Data daterr/'Too many variables listed', 'Quote mark missing', &
        'Slash not found', 'Command sequence problem', 'Too few variables', &
        'Columns illegal', 'Cannot read beyond column 200', &
        'Data file not found'/
      Data makerr/'Illegal use of WEIGHT variable', 'Check punctuation', &
        'Too many variables', 'Independent variable(s) not found', &
        'Missing DEPENDENT subcommand', 'Missing VARIABLES subcommand', &
        'Subcommands out of sequence', 'Missing equals sign', &
        'Incomplete statement', 'Invalid minimum or maximum values', &
        'Too many categories of variables'/
      Data sterr/'Illegal use of WEIGHT variable', &
        'Too many BREAKDOWN variables listed', &
        'All variables must be from MAKETAB command', &
        'Too many CONTROL variables listed', &
        'Syntax of STANDARD=CATEGORY subcommand', &
        'Cannot understand STANDARD subcommand', &
        'BREAKDOWN variable(s) must be specified', &
        'Specified STANDARD category out of range', &
        'Specified STANDARD category has no cases', &
        'All cases excluded from the analysis', 'Equals sign missing', &
        'BREAKDOWN subcommand missing', 'CONTROL subcommand missing', &
        'Option permits only one breakdown category', &
        'CONTROL variable(s) must be specified'/
      Data decerr/'Illegal use of WEIGHT variable', &
        'All variables must be from MAKETAB command', &
        'Too many FACTORS listed', 'Problem with comparison categories', &
        'Comparison category out of variable range', &
        'Too many categories of variables for DECOMP', &
        'No cases found in at least one population', 'Equals sign missing', &
        'FACTOR subcommand not found', 'COMPARE subcommand not found'/
      Data recerr/'Too many recodes', 'Punctuation problem', &
        'Range of values invalid', 'ELSE operator illegal', &
        'Incomplete statement', 'Check syntax'/
      Data selerr/'Operator illegal or missing', &
        'Too many SELECT IF statements', 'Incomplete statement'/
      Data generr/': Command too long', ' in WEIGHT command', &
        ': Array index exceeds 640000', ': Command file not found', &
        ': Command must begin in column 1', ': Warning count exceeded'/
      Data comerr/'Punctuation problem', 'Too many variables', &
        'Invalid minimum or maximum values', 'Incomplete statement', &
        'You must have two or more variables', &
        'Too many categories of variables'/
!
!         DEPENDING ON THE ERROR NUMBER, WRITE THE MESSAGE FROM THE
!         APPROPRIATE ARRAY.
!
      If (errnum<200) Then
        index = errnum - 100
        Write (*, 100) errnum, daterr(index)
        Write (23, 100) errnum, daterr(index)
100     Format (/, /, /, ' ERROR ', I3, 1X, 'in DATA LIST command: ', A)
      Else If (errnum<300) Then
        index = errnum - 200
        Write (*, 110) errnum, makerr(index)
        Write (23, 110) errnum, makerr(index)
110     Format (/, /, /, ' ERROR ', I3, 1X, 'in MAKETAB command: ', A)
      Else If (errnum<400) Then
        index = errnum - 300
        Write (*, 120) errnum, sterr(index)
        Write (23, 120) errnum, sterr(index)
120     Format (/, /, /, ' ERROR ', I3, 1X, 'in STANDARDIZE command: ', A)
      Else If (errnum<500) Then
        index = errnum - 400
        Write (*, 130) errnum, decerr(index)
        Write (23, 130) errnum, decerr(index)
130     Format (/, /, /, ' ERROR ', I3, 1X, 'in DECOMPOSE command: ', A)
      Else If (errnum<600) Then
        index = errnum - 500
        Write (*, 140) errnum, recerr(index)
        Write (23, 140) errnum, recerr(index)
140     Format (/, /, /, ' ERROR ', I3, 1X, 'in RECODE command: ', A)
      Else If (errnum<700) Then
        index = errnum - 600
        Write (*, 150) errnum, selerr(index)
        Write (23, 150) errnum, selerr(index)
150     Format (/, /, /, ' ERROR ', I3, 1X, 'in SELECT IF command: ', A)
      Else If (errnum<800) Then
        index = errnum - 700
        Write (*, 160) errnum, generr(index)
        Write (23, 160) errnum, generr(index)
160     Format (/, /, /, ' ERROR ', I3, A)
      Else If (errnum<900) Then
        index = errnum - 800
        Write (*, 170) errnum, comerr(index)
        Write (23, 170) errnum, comerr(index)
170     Format (/, /, /, ' ERROR ', I3, 1X, 'in COMBINE command: ', A)
      End If
!
!         FATERR MEANS FATAL ERROR
!
      Call faterr
    End Subroutine



    Subroutine errnam(errnum, varname)
!
!         THIS SUBROUTINE HANDLES THE ERROR MESSAGES IN WHICH
!         A VARIABLE NAME MUST BE WRITTEN.
!
      Implicit Integer (A-Z)
      Character *10 varname
      Character *11 command
      Dimension command(7)
      Data command/'  MAKETAB ', 'STANDARDIZE', ' DECOMPOSE', '   RECODE  ', &
        ' SELECT IF ', '  WEIGHT  ', '  COMBINE '/
!
!         MOST OF THE ERRORS (NUMBERS 1-7) ARE SIMPLY THAT A
!         VARIABLE NAME WAS NOT FOUND.
!
      If (errnum<100) Then
        num = errnum*100 + 150
!
!         NUM MAKES THE ERROR NUMBERS CONSISTENT WITH THOSE USED IN
!         SUBROUTINE ERROR
!
        Write (*, 100) num, command(errnum), varname
        Write (23, 100) num, command(errnum), varname
100     Format (/, /, 1X, 'ERROR ', I3, ' in ', A, ' Command: Variable ', A, &
          ' not found')
!
!         THE ONLY OTHER MESSAGE PERTAINS TO MIN AND MAX IN MAKETAB.
!
      Else If (errnum==251) Then
        Write (*, 110) errnum, varname
        Write (23, 110) errnum, varname
110     Format (/, /, 1X, 'ERROR ', I3, ' in MAKETAB command: Minimum', &
          ' and maximum for ', A10, 'must be specified')
      End If
      Call faterr
    End Subroutine

    Subroutine faterr
!
!         THIS LITTLE ROUTINE JUST WRITES A TERMINAL MESSAGE AND DIES.
!
      Write (*, 100)
      Write (23, 100)
!
!  ****** WARNING: THE FOLLOWING STATEMENT WRITES THE BELL CHARACTER.
!                  THIS WILL NOT WORK ON ALL SYSTEMS.
!
!cc        WRITE(6,*)CHAR(007)
100   Format (/, 1X, 'Error is fatal - program terminated')
      Close (22, Status='DELETE')
!cc        CLOSE(6,STATUS='DELETE')
      Stop ' '
    End Subroutine

    Integer Function select(rec, beg, end)
!
!  THIS FUNCTION CONVERTS SELECTED COLUMNS OF CHARACTER STRINGS INTO INTEGERS
!
      Implicit Integer (A-Z)
      Character *500 rec
      Dimension digit(8)

      numdig = 0
      Do i = 1, 8
        digit(i) = 0
      End Do
!
! ****** WARNING: THE FOLLOWING ASSUMES THAT THE SYSTEM IS RUNNING ON
!                 THE EXTENDED ASCII CHARACTER CODES.  THIS WILL HAVE
!                 TO BE MODIFIED FOR COMPUTERS USING OTHER CODES.
!
!     CALCULATE THE VALUE OF EACH DIGIT SEPARATELY.
!
      Do i = beg, end
        c = ichar(rec(i:i))
!
!     CONVERT BLANKS TO ZEROS.
!
        If (c==32) c = 48
!
!     CONVERT SLASH CHARACTER TO NULL
!
        If (c==47) c = 0
!
!     IF THE CODE IS LESS THAN A PERIOD OR GREATER THAN A NINE,
!     WRITE AN ERROR MESSAGE.
!
        If (c<46 .Or. c>57) Then
          Write (23, 100) rec(i:i), ichar(rec(i:i)), i
          Write (*, 100) rec(i:i), ichar(rec(i:i)), i
100       Format (/, 1X, 'ERROR 7: Non-integer found where integer &
            &expected in preceding command.', /, 10X, 'The character is "', &
            A1, '", ASCII decimal code is ', I3, ' at column', I4, /)
!
!     OTHERWISE, IF THE CHARACTER IS NOT A PERIOD, CONVERT IT TO
!     THE APPROPRIATE INTEGER. IGNORE PERIODS, AND WRITE A WARNING
!     IF THE INTEGER HAS MORE THAN EIGHT CHARACTERS.
!
        Else If (c/=46) Then
          numdig = numdig + 1
          If (numdig<=8) Then
            digit(numdig) = c - 48
          Else
            Write (23, 110)
            Write (*, 110)
110         Format (3X, 'WARNING: Value over 8 characters truncated')
            numdig = numdig - 1
          End If
        End If
      End Do
!
!     NOW ADD UP THE DIGITS INTO A SINGLE INTEGER.
!
      select = 0
      Do i = 1, numdig
        place = numdig - i
        select = select + digit(i)*(10**place)
      End Do
      Return
    End Function

    Integer Function sdata(rec, beg, end)
!
!  THIS FUNCTION CONVERTS SELECTED COLUMNS OF CHARACTER STRINGS INTO INTEGERS
!  IT IS STREAMLINED FOR SPEED
!
      Implicit Integer (A-Z)
      Character *200 rec

      sdata = 0
      factor = 100000000
!
! ****** WARNING: THE FOLLOWING ASSUMES THAT THE SYSTEM IS RUNNING ON
!                 THE EXTENDED ASCII CHARACTER CODES.  THIS WILL HAVE
!                 TO BE MODIFIED FOR COMPUTERS USING OTHER CODES.
!
!     CALCULATE THE VALUE OF EACH DIGIT SEPARATELY.
!
      Do i = beg, end
        c = ichar(rec(i:i))
!
!     CONVERT BLANKS TO ZEROS.
!
        If (c==32) c = 48
!
!     CONVERT SLASH CHARACTER TO NULL
!
        If (c==47) c = 0
!
!     IF THE CODE IS LESS THAN A PERIOD OR GREATER THAN A NINE,
!     WRITE AN ERROR MESSAGE.
!
        If (c<46 .Or. c>57) Then
          Write (23, 100) rec(i:i), ichar(rec(i:i)), i
          Write (*, 100) rec(i:i), ichar(rec(i:i)), i
100       Format (/, 1X, 'ERROR 6: Illegal character found in data.', /, 10X, &
            'The character is "', A1, '", ASCII decimal Code is ', I3, &
            ' at column', I4, /)
!
!     OTHERWISE, IF THE CHARACTER IS NOT A PERIOD, CONVERT IT TO
!     THE APPROPRIATE INTEGER. IGNORE PERIODS.
!
        Else If (c/=46) Then
          factor = factor/10
          sdata = sdata + (c-48)*factor
        End If
      End Do
!
!     NOW ADD UP THE DIGITS INTO A SINGLE INTEGER.
!
      sdata = sdata/factor
      Return
    End Function

    Subroutine datlist(list, len, files, varnam, vardes, varnum, maxno)
!
!   THIS SUBROUTINE INTERPRETS THE DATA LIST COMMAND.  DETAILS ON VARIABLES
!   ARE STORED IN ARRAY VARDES (VARIABLE DESCRIPTIONS) AS FOLLOWS:
!
!     VARDES: VARIABLE DESCRIPTIONS
!             ROWS  - VARIABLE NUMBER, IN ORDER SPECIFIED
!             COL 1 - BEGINS COLUMN
!             COL 2 - ENDS COLUMN
!
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *35 files
      Character *500 list
      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension files(4)

      chr = 0
      quote1 = 0
      quote2 = 0
!
!       FIND DATA FILE NAME
!
      beg = 0
      end = 0
!
!       LOOK FOR FIRST QUOTE MARK -- IF WE GET TO THE END OF THE LINE,
!       GO TO END
!
100   chr = chr + 1
      If (chr>len+1) Go To 170
      If (ichar(list(chr:chr))/=39) Go To 100
      chr = chr + 1
      If (chr>len+1) Go To 170
      beg = chr
      quote1 = chr
!
!       LOOK FOR THE SECOND QUOTE THE SAME WAY
!
110   chr = chr + 1
      If (chr>len+1) Go To 170
      If (ichar(list(chr:chr))/=39) Go To 110
      end = chr - 1
      quote2 = chr
!
!       NOW STICK THE FILE NAME IN POSITION 2 OF ARRAY FILES
!
      files(2) = list(beg:end)

!
!       FIND SLASH
!
120   chr = chr + 1
      If (chr>len+1) Go To 170
      If (list(chr:chr)/='/') Go To 120
!
!       FIND EACH VARIABLE
!
      varnum = 0
130   Continue
      varnum = varnum + 1
!
!       NO MORE THAN 30 VARIABLES MAY BE DEFINED
!
      If (varnum>30) Then
        Call error(101)
      End If
!
!       THE CHARACTER | SERVES AS A DELIMITER (SEE SUBROUTINE GETLIST)
!       SET BEG TO NEXT CAHARACTER AFTER THE SECOND QUOTE.
!       IF IT IS A DELIMITER, GO TO THE NEXT CHARACTER.
!
      beg = chr + 1
      chr = chr + 1
      If (chr>len+1) Go To 170
      If (list(chr:chr)=='|') Then
        beg = chr + 1
      End If
140   chr = chr + 1
      If (chr>len+1) Go To 170
!
!       NOW LOOK FOR NEXT DELIMITER.
!
      If (list(chr:chr)/='|') Go To 140
      end = chr - 1
!
!       THIS SHOULD DEFINE A VARIABLE NAME
!
      varnam(varnum) = list(beg:end)
!
!       THE NEXT COLUMNS AFTER THE LAST DELIMITER SHOULD BE A COLUMN NUMBER.
!       IT WILL BE TERMINATED BY A DASH OR A DELIMITER
!
      beg = chr + 1
150   chr = chr + 1
      If (chr>len+1) Go To 170
      If (list(chr:chr)/='-' .And. list(chr:chr)/='|' .And. chr/=len+1) &
        Go To 150
!
!       STICK THE COLUMN NUMBER IN THE APPROPRIATE POSITION OF
!       ARRAY VARDES AFTER CONVERTING IT TO AN INTEGER USING
!       FUNCTION SELECT.
!
      end = chr - 1
      vardes(varnum, 1) = select(list, beg, end)
      If (list(chr:chr)=='|' .Or. chr==len+1) Then
        vardes(varnum, 2) = vardes(varnum, 1)
        Go To 130
      End If
!
!       IF THE COLUMN NUMBER WAS TERMINATED WITH A DASH, LOOK FOR ANOTHER.
!
      beg = chr + 1
160   chr = chr + 1
      If (chr>len+1) Go To 170
      If (list(chr:chr)/='|' .And. chr/=len+1) Go To 160
!
!       STICK IT IN THE APPROPRIATE PLACE.
!
      end = chr - 1
      vardes(varnum, 2) = select(list, beg, end)
!
!       NOW GO UP AND LOOK FOR ANOTHER VARIABLE.
!
      Go To 130
170   Continue
!
!       ADJUST VARNUM; WE'VE COUNTED ONE TOO MANY.
!
      varnum = varnum - 1
!
!       NOW CHECK VALIDITY OF DATA LIST COMMAND.
!       THE FOLLOWING BLOCK CHECKS THE NUMBER OF VARIABLES, THE
!       EXISTANCE OF A SLASH, AND THE EXISTENCE AND ORDER OF
!       QUOTE MARKS.
!
      slash = index(list, '/')
      If (quote2==0) Then
        Call error(102)
      Else If (quote2-quote1>35) Then
        Call error(102)
      Else If (slash==0) Then
        Call error(103)
      Else If (slash<quote2) Then
        Call error(104)
      Else If (varnum<3) Then
        Call error(105)
      End If
!
!      NOW CHECK TO SEE THAT INDICATED COLUMNS ARE LEGIT.
!
      maxno = 0
      Do i = 1, varnum
        dif = vardes(i, 2) - vardes(i, 1)
        If (dif<0 .Or. dif>8 .Or. vardes(i,1)==0 .Or. vardes(i,2)==0) Then
          Call error(106)
        End If
        If (vardes(i,2)>maxno) maxno = vardes(i, 2)
      End Do
      If (maxno>200) Then
        Call error(107)
      End If
      Return
    End Subroutine



    Subroutine caseconv(line)
!
!       THIS SUBROUTINE CONVERTS LOWER CASE LETTERS TO CAPS
!
      Implicit Integer (A-Z)
      Character *80 line

      numquot = 0
      lca = ichar('a')
      lcz = ichar('z')
      difcas = lca - ichar('A')
      Do i = 1, 80
        decval = ichar(line(i:i))
        If (decval==39) numquot = numquot + 1
        If (numquot/=1) Then
          If (decval>=lca .And. decval<=lcz) Then
            decval = decval - difcas
            line(i:i) = char(decval)
          End If
        End If
      End Do
      Return
    End Subroutine


    Subroutine getlist(command, list, len, done, pg, ln, setcmd, prnum)
      Implicit Integer (A-Z)
!
!  THIS SUBROUTINE GETS A COMMAND, REMOVES THE SPACES, STICKS IT ALL IN ONE
!  BIG STRING, AND STICKS IN ALTERNATE DELIMITERS FOR SPACES.
!
      Character *3 scrnam
      Character *5 command
      Character *80 line
      Character *80 setcmd
      Character *500 list
      Dimension setcmd(10)
!cc        INQUIRE(UNIT=6,NAME=SCRNAM)
      scrnam = 'CON'
!
      Do i = 1, 500
        list(i:i) = ' '
      End Do
      done = 0
100   len = 0
      command = 'NULL '
!
!         READ FIRST LINE OF COMMAND
!
      If (prnum<10) Then
        prnum = prnum + 1
        line = setcmd(prnum)
      Else
        Read (3, 110, End=160) line
      End If
110   Format (A80)
!
!         CONVERT LOWER CASE TO CAPS ON FIRST LINE OF COMMAND
!
      Call caseconv(line)
!
!           SKIP THE LINE IF IT IS A COMMENT OR BLANK, BUT
!           ECHO IT TO THE FILE IF IT BEGINS WITH AN ASTERISK
!
      If (line(1:1)=='*') Then
        Call page(pg, ln)
        Write (23, 130) line(1:79)
        Go To 100
      End If
      If (line(1:1)=='[' .Or. line(1:50)== &
        '                                                  ') Go To 100
!
!           ECHO LINE TO OUTPUT FILE AND SCREEN.  IF THIS IS THE
!           FIRST LINE ON A PAGE, SKIP A SPACE.
!
!           SUBROUTINE PAGE CONTROLS PAGE NUMBERS AND HEADERS.
!
      Call page(pg, ln)
      If (ln==1) Then
        ln = ln + 1
        Write (23, *)
      End If
      Write (23, 130) line(1:79)
      If (scrnam=='CON') Write (*, 120) line(1:79)
120   Format (1X, A79)
130   Format (A79)
!
!           CHECK FOR SOME ALTERNATE COMMAND FORMS
!
      command = line(1:5)
      If (command=='DATAL') Then
        command = 'DATA '
      Else If (command=='MAKE ') Then
        command = 'MAKET'
      Else If (command=='SETLI') Then
        command = 'SET L'
      Else If (command=='SETRE') Then
        command = 'SET R'
      End If
!
!           DO NO ALTERATIONS FOR SET LISTING, SET RESULTS,
!           OR WEIGHT COMMANDS
!
      If (command=='SET L' .Or. command=='SET R' .Or. command=='WEIGH') Then
        list(1:80) = line
        Go To 170
      End If
!
!          CHANGE SPACE DELIMITERS TO ~ IN SELECT IF COMMAND
!
!  ****** WARNING: THE ~ CHARACTER IS NOT AVAILABLE ON SOME SYSTEMS.
!         YOU MAY HAVE TO CHANGE IT FOR SOMETHING ELSE
!
      If (command=='SELEC') Then
        beg = index(line, ' GE ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~GE~'
        beg = index(line, ' GT ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~GT~'
        beg = index(line, ' LE ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~LE~'
        beg = index(line, ' LT ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~LT~'
        beg = index(line, ' EQ ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~EQ~'
        beg = index(line, ' NE ')
        end = beg + 3
        If (beg>0) line(beg:end) = '~NE~'
      End If
!
!        CHANGE SPACE DELIMITERS IN DATA LIST COMMAND TO |
!        ALSO CHANGE PERIOD IN FILENAME IN DATA LIST COMMAND TO A ^,
!        IN ORDER TO AVOID CONFUSION WITH THE END-OF-COMMAND PERIOD.
!
!  ****** WARNING: THE | CHARACTER IS NOT AVAILABLE ON SOME SYSTEMS.
!         YOU MAY HAVE TO CHANGE IT FOR SOMETHING ELSE
!
!
      quote1 = 0
      quote2 = 0
140   Do i = 1, 80
        If (command=='DATA ' .And. i<80) Then
          If (i==1) Then
            len = len + 1
            list(len:len) = '|'
          End If
          If (line(i:i)==' ' .And. line(i+1:i+1)/=' ') Then
            line(i:i) = '|'
          End If
        End If
        If (line(i:i)=='''' .And. quote1==0) Then
          quote1 = i
        Else If (line(i:i)=='''' .And. quote1/=0) Then
          quote2 = i
        Else If (line(i:i)=='.' .And. quote1/=0 .And. quote2==0) Then
          line(i:i) = '^'
        Else If (line(i:i)=='|' .And. quote1/=0 .And. quote2==0) Then
          line(i:i) = ' '
        End If
!
!   NOW MOVE THE COLUMN FROM STRING LINE TO THE APPROPRIATE POSITION IN
!   THE STRING LIST.  MAXIMIMUM COMMAND LENGTH IS 500 CHARACTERS, ONCE
!   THE SPACES ARE REMOVED.
!
        If (line(i:i)=='.') Then
          Go To 170
        Else If (line(i:i)/=' ') Then
          len = len + 1
          If (len>499) Then
            Call error(701)
          End If
          list(len:len) = line(i:i)
        End If
      End Do
!
!         READ NEXT LINE OF THE COMMAND.
!
150   If (prnum<10) Then
        prnum = prnum + 1
        line = setcmd(prnum)
      Else
        Read (3, 110, End=160) line
      End If
!
!         CONVERT LOWER CASE TO CAPS FOR SUBSEQUENT LINES
!
      Call caseconv(line)
!
!         SKIP COMMENTS
!
      If (line(1:1)=='*') Then
        Call page(pg, ln)
        Write (23, 130) line(1:79)
        Go To 150
      End If
      If (line(1:1)=='[' .Or. line(1:50)== &
        '                                                  ') Go To 150
!
!           ECHO LINE TO OUTPUT FILE AND SCREEN
!
      Call page(pg, ln)
      Write (23, 130) line(1:79)
      If (scrnam=='CON') Write (*, 120) line(1:79)

      Go To 140
160   done = 1
170   Continue
      Return
    End Subroutine



    Subroutine get(line, maxno, done)
      Implicit Integer (A-Z)
!
!  THIS SUBROUTINE GETS DATA STRINGS.
!
      Character *200 line

      Read (8, 100, End=110) line(1:maxno)
100   Format (A)
      Return
!
!   IF THE END OF THE DATA FILE IS ENCOUNTERED, SO INDICATE.
!
110   done = 1
      Return
    End Subroutine


    Subroutine tab(vardes, recode, seldes, comdes, tratab, files, maxno, &
      varnam, tranum, varnum, wtdig, depdig, wrtab, wtdiv, pg, ln)
!
!  THIS IS THE MAIN TABULATION ROUTINE.
!
      Implicit Integer (A-Z)
      Integer i, j, k, l, m
      Real table, wt, tnum, tden
      Character *3 scrnam
      Character *10 varnam
      Character *10 indnam
      Character *10 badvar
      Character *35 files
      Character *200 record

      Dimension varval(30)
      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension files(4)
      Dimension recode(200, 4)
      Dimension seldes(50, 3)
      Dimension comdes(10, 5)
      Dimension tratab(300, 2)
      Dimension indval(5)
      Dimension indnam(5)
      Dimension wrform(7, 2)
      Common /tabs/table(640000, 2), maxvars(5)

!       INITIALIZE WORKING TABLES

      Do i = 1, 640000
        table(i, 1) = 0.0
        table(i, 2) = 0.0
      End Do
      Do i = 1, 5
        maxvars(i) = 1
      End Do
!
!       SET WARNING COUNT TO ZERO
!
      warn = 0
!
!       CALCULATE MAXIMUM RANGE OF EACH INDEPENDENT VARIABLE,
!       AND STICK IT IN ARRAY MAXVARS.
!       COUNT INDEPENDENT VARIABLES AND STORE THEIR NAMES IN
!       ARRAY INDNAM, AND LOCATE THE POSITION OF THE DEPENDENT
!       VARIABLE AND THE WEIGHT VARIABLE, IF IT EXISTS.
!
!
      wtpos = 999
      deppos = 999
      numind = 0
      Do i = 1, varnum
        If (vardes(i,5)>1 .And. vardes(i,5)<=6) Then
          maxvars(vardes(i,5)-1) = vardes(i, 4) - vardes(i, 3) + 1
          indnam(vardes(i,5)-1) = varnam(i)
          numind = numind + 1
        Else If (vardes(i,5)==1) Then
          deppos = i
        Else If (vardes(i,5)==7) Then
          wtpos = i
        End If
      End Do
!
!       MAKE SURE THE DEPENDENT VARIABLE IS PROPERLY
!       IDENTIFIED
!
      If (deppos==999) Call error(212)
!
!       SEE IF THE TABLE WILL EXCEED THE SIZE LIMIT.
!
      biggest = maxvars(1)*maxvars(2)*maxvars(3)*maxvars(4)*maxvars(5)
      If (biggest>640000) Call error(211)
!
!       CALCULATE DIVISION FACTOR RELEVANT TO WTDIG (WEIGHT),
!       BASED ON THE NUMBER OF DIGITS TO THE RIGHT OF THE DECIMAL
!       POINT.
!
      wtdiv = 1
      If (wtdig==1) Then
        wtdiv = 10
      Else If (wtdig>1) Then
        wtdiv = 10**wtdig
      End If
!
!       DO THE SAME FOR DEPDIG (DEPENDENT VARIABLE DIVISION FACTOR)
!
      depdiv = 1
      If (depdig==1) Then
        depdiv = 10
      Else If (depdig>1) Then
        depdiv = 10**depdig
      End If
!
!       OPEN RAW DATA FILE AND INITIALIZE READ SUBROUTINE
!
!  ****** WARNING: MOST COMPILERS WILL NOT ACCEPT THE BLOCKSIZE
!         OPTION ON THE OPEN COMMAND.
!
      Open (8, File=files(2), Status='OLD', Err=100)
      done = 0
!
!       IF THE RAW DATA FILE WAS NOT FOUND, WRITE ERROR MESSAGE
!
      Go To 110
100   Continue
      Call error(108)
110   Continue
!
!       BEGIN PROCESSING
!       BEGIN MAIN LOOP.  FIRST SKIP A COUPLE LINES ON THE SCREEN,
!       AND WRITE THE CASE NUMBER EVERY 100 LINES READ.
!
      scrnam = 'CON'
!cc        INQUIRE(UNIT=6,NAME=SCRNAM)
      If (scrnam=='CON') Write (*, 120)
120   Format (/)
      Rewind 8
  !  generate carriage return
      Do case = 1, 50000
        If (scrnam=='CON')then
           Write (*, '(a,a4, I10)', advance='no')char(13), 'Case', case*100000 - 100000
        end if   
!        If (scrnam=='CON') Write (*, 130) case*100000 - 100000
        
        
!         write(*,'(a41)', advance='no') "Enter the name of your command file --->  "
        
!130     Format ('Case', I10)
        Do num = 1, 100000
!
!       SUBROUTINE GET PUTS THE DATA INTO THE STRING  "RECORD"
!       WHEN DONE, JUMP OUT OF LOOP.
!
          Call get(record, maxno, done)
          If (done==1) Go To 160
!
!       CONVERT VARIABLES TO INTEGERS, USING FUNCTION "SDATA"
!
          Do i = 1, varnum
            If (vardes(i,1)/=9999) Then
              If (vardes(i,3)==0 .And. vardes(i,4)==0) vardes(i, 3) = 1
              varval(i) = sdata(record, vardes(i,1), vardes(i,2))
            End If
          End Do
!
!       DO TRANSFORMATIONS
!
          Call trans(tratab, tranum, seldes, recode, comdes, varval, vardes, &
            varnum, range, quit)
!
!       IF THIS CASE SHOULD BE SKIPPED BECAUSE OF A SELECT IF
!       COMMAND, JUMP TO THE END OF THE LOOP.
!
          If (quit==1) Go To 150
!
!       FIND DEPENDENT VARIABLE, INDEPENDENT VARIABLES,
!       WEIGHT, AND POPULATION DEFINITION VARIABLE.
!       PUT THE VALUE OF THE INDEPENDENT VARIABLES IN
!       ARRAY INDVAL, AND THE OTHERS IN DEPVAL AND WTVAL.
!
          wtval = 1
          Do i = 1, 5
            indval(i) = 1
          End Do
!
!
          depval = varval(deppos)
          If (wtpos<=varnum) wtval = varval(wtpos)
          Do i = 1, varnum
            If (vardes(i,5)>=2 .And. vardes(i,5)<=6) indval(vardes(i,5)-1) &
              = varval(i)
          End Do
!
!       INSURE THAT THE VALUES OF THE INDEPENDENT VARIABLES FALL
!       WITHIN SPECIFIED RANGE, AND WRITE WARNING IF NOT.
!
          Do i = 1, numind
            If (indval(i)<1 .Or. indval(i)>maxvars(i) .Or. range>0) Then
              warn = warn + 1
              ln = ln + 4
              Call page(pg, ln)
              badvar = indnam(i)
              If (range>0) badvar = varnam(range)
              Write (23, 140) badvar, case*100 - 100 + num
              Write (*, 140) badvar, case*100 - 100 + num

140           Format (/, 3X, 'WARNING: Variable ', A10, 'outside ', &
                'specified range at case number ', I6, /, &
                '            Case ignored', /, /)
              If (warn>25) Call error(706)
              Go To 150
            End If
          End Do
!
!       CALCULATE WEIGHT
!
          wt = float(wtval)/float(wtdiv)
!
!       CALCULATE TABULATION INDEX AND TABULATE.  FUNCTION IND WILL
!       GENERATE THE INDEX.
!
          a = indval(1)
          b = indval(2)
          c = indval(3)
          d = indval(4)
          e = indval(5)
          index = ind(a, b, c, d, e)
          table(index, 2) = table(index, 2) + wt
          table(index, 1) = table(index, 1) + (float(depval)/float(depdiv))*wt
150     End Do
      End Do
160   Continue
!
!       IF THE USER HAS SPECIFIED THE WRITE TABLE OPTION, CALL SUBROUTINE
!       TO PREPARE FOR IT.
!
      If (wrtab>0) Call wrprep(varnum, vardes, wrform, files, numind, indnam, &
        wtdig, pg, ln)
!
!       WRITE TABULATED DATA TO TEMPORARY WORK FILE FOR LATER USE
!
! ****** WARNING: MOST COMPILERS WILL NOT ACCEPT THE BLOCKSIZE OPTION,
!        AND SOME DON'T LIKE FILENAME EXTENSIONS OR BINARY FILES.
!        IF BINARY FILES ARE NOT ALLOWED, USE UNFORMATTED.
!
      Open (22, File='WRKFILE.TMP', Form='UNFORMATTED')
      Rewind 22
      Do i = 1, maxvars(1)
        Do j = 1, maxvars(2)
          Do k = 1, maxvars(3)
            Do l = 1, maxvars(4)
              Do m = 1, maxvars(5)
                a = i
                b = j
                c = k
                d = l
                e = m
                index = ind(a, b, c, d, e)
                If (table(index,2)>0.0) Then
                  Write (22) i, j, k, l, m, table(index, 1), table(index, 2)
!
!    WRITE ASCII FILE IF WRITE TABLE OPTION IS SPECIFIED
!    (FIRST CONVERT INDEPENDENT VALUES)
!
                  If (wrtab>0) Then
                    tnum = table(index, 1)
                    tden = table(index, 2)
                    Call wrdata(varnum, vardes, wrform, i, j, k, l, m, numind, &
                      wtdig, wtdiv, tnum, tden)
                  End If
                End If

              End Do
            End Do
          End Do
        End Do
      End Do

      Return
    End Subroutine



    Subroutine trans(tratab, tranum, seldes, recode, comdes, varval, vardes, &
      varnum, range, quit)
      Implicit Integer (A-Z)
!
!   THIS SUBROUTINE CARRIES OUT THE DATA TRANSFORMATIONS.
!
      Dimension varval(30)
      Dimension recode(200, 4)
      Dimension seldes(50, 3)
      Dimension comdes(10, 5)
      Dimension tratab(300, 2)
      Dimension vartem(30)
      Dimension vardes(30, 5)
      Dimension comval(4, 2)
!
!       MOVE VALUES TO ARRAY VARTEM TEMPORARILY
!
      Do i = 1, varnum
        vartem(i) = varval(i)
      End Do
!
!       INITIALIZE FLAGS.
!
      range = 0
      quit = 0
      Do t = 1, tranum
!
!       FOR EACH TRANSFORMATION, FIND OUT WHAT TYPE IT IS AND DO IT.
!
        If (tratab(t,1)==1) Then
!
!       DO SELECT IFS (SEE SELDES DESCRIPTIONS)
!
          s = tratab(t, 2)
          If (seldes(s,2)==1) Then
            If (vartem(seldes(s,1))/=seldes(s,3)) quit = 1
          Else If (seldes(s,2)==2) Then
            If (vartem(seldes(s,1))<=seldes(s,3)) quit = 1
          Else If (seldes(s,2)==3) Then
            If (vartem(seldes(s,1))>=seldes(s,3)) quit = 1
          Else If (seldes(s,2)==4) Then
            If (vartem(seldes(s,1))<seldes(s,3)) quit = 1
          Else If (seldes(s,2)==5) Then
            If (vartem(seldes(s,1))>seldes(s,3)) quit = 1
          Else If (seldes(s,2)==6) Then
            If (vartem(seldes(s,1))==seldes(s,3)) quit = 1
          End If

        Else If (tratab(t,1)==2 .And. quit==0) Then
!
!       DO RECODES (SEE RECODE DESCRIPTIONS)
!
          r = tratab(t, 2)
          If (varval(recode(r,1))>=recode(r,2) .And. varval(recode(r, &
            1))<=recode(r,3)) vartem(recode(r,1)) = recode(r, 4)
!
!       DO COMBINATIONS.  THIS IS A LITTLE MORE COMPLICATED.
!       A NEW VARIABLE IS CREATED FROM UP TO FOUR OTHERS.
!
        Else If (tratab(t,1)==3 .And. quit==0) Then
          comnum = tratab(t, 2)
!
!       ARRAY COMVAL (4,2) WILL HOLD THE VALUES AND RANGES OF
!       THESE VARIABLES.  THE ROWS INDICATE VARIABLE NUMBERS,
!          COLUMN 1 WILL CONTAIN THE RANGES OF THESE VARIABLES;
!          COLUMN 2 WILL CONTAIN THEIR VALUES.
!
!       VARIABLE LOC INDICATES POSITION IN VARDES OF THE CURRENT
!       VARIABLE TO BE COMBINED
!
!       VARIABLE RANGE WILL BE SET TO THE CURRENT VARIABLE NUMBER
!       IF THE VALUE OF THE VARIABLE IS OUT OF THE SPECIFIED MINIMUM
!       OR MAXIMUM.
!
!
!       THE FOLLOWING LOOP LOCATES RANGES AND VALUES AND STICKS THEM
!       IN ARRAY COMVAL
!
          Do c = 1, 4
            comval(c, 1) = 1
            comval(c, 2) = 1
            loc = 30
            If (comdes(comnum,c)/=99) Then
              loc = comdes(comnum, c)
              comval(c, 1) = vardes(loc, 4) - vardes(loc, 3) + 1
              If (vartem(loc)<vardes(loc,3) .Or. vartem(loc)>vardes(loc,4)) &
                Then
                range = loc
              End If
              comval(c, 2) = vartem(loc) - vardes(loc, 3) + 1
            End If
          End Do
!
!      FUNCTION IND2 WILL CALCULATE THE NEW VALUE OF THE COMBINED
!      VARIABLE, BASED ON THE VALUES AND RANGES DETERMINED ABOVE.
!      THE NEW VALUE IS PLACED IN THE INDICATED POSITION FOR THE NEW
!      VARIABLE.
!
          newloc = comdes(comnum, 5)
          vartem(newloc) = ind2(comval(4,2), comval(3,2), comval(2,2), &
            comval(1,2), comval(4,1), comval(3,1), comval(2,1))
          varval(newloc) = vartem(newloc)
        End If
!
!       TRANSFORMATION LOOP ENDS HERE
!
      End Do
!
!       RETURN VALUES FROM VARTEM TO VARVAL
!
      Do i = 1, varnum
        varval(i) = vartem(i) - vardes(i, 3) + 1
      End Do

      Return
    End Subroutine




    Subroutine set(command, list, files, pg, ln)
!
!   THIS SUBROUTINE INTERPRETS THE SET COMMANDS.
!
      Implicit Integer (A-Z)
      Character *5 command
      Character *35 files
      Character *500 list
      Dimension files(4)
!
!       FIND QUOTATION MARKS
!
      quote1 = 0
      quote2 = 0
      Do i = 1, 80
        If (list(i:i)=='''' .And. quote1==0) Then
          quote1 = i
        Else If (list(i:i)=='''' .And. quote1/=0) Then
          quote2 = i
        End If
      End Do
!
!       IF THEY ARE NOT BOTH FOUND, WRITE ERROR MESSAGE
!
      If (quote1==0 .Or. quote2==0 .Or. quote1+1==quote2) Then
        ln = ln + 2
        Call page(pg, ln)
!
!  ****** WARNING: THE FOLLOWING STATEMENT WRITES THE BELL CHARACTER.
!                  THIS WILL NOT WORK ON ALL SYSTEMS.
!
!cc          WRITE(6,*)CHAR(007)
        Write (23, 100)
        Write (*, 100)
100     Format (/, 1X, 'ERROR 4 in SET command: Command ignored', /)
        Return
      End If
!
!       OTHERWISE, GET NEW NAME AND STICK IT IN THE APPROPRIATE
!       POSITION OF ARRAY FILES
!           POSITION 3 FOR SET LISTING COMMAND
!           POSITION 4 FOR SET RESULTS COMMAND
!
      beg = quote1 + 1
      end = quote2 - 1
      fnum = 3
      If (command=='SET R') fnum = 4
      files(fnum) = list(beg:end)
      Return
    End Subroutine







    Subroutine recdat(list, len, recode, varnam, recnum, tratab, tranum)
!
!   THIS SUBROUTINE INTERPRETS THE RECODE COMMAND.  THE RECODES ARE STORED
!   IN THE ARRAY RECODE, AS FOLLOWS:
!
!      RECODE: RECODE DESCRIPTIONS
!              ROWS  - RECODE NUMBER
!              COL 1 - VARIABLE NUMBER
!              COL 2 - GE VALUE
!              COL 3 - LE VALUE
!              COL 4 - NEW VALUE
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *500 list
      Character *10 thisv
      Dimension recode(200, 4)
      Dimension tratab(300, 2)
      Dimension varnam(30)
!
!       FIND FIRST VARIABLE NAME; IT SHOULD BEGIN RIGHT AFTER THE
!       KEYWORD RECODE, AND END RIGHT BEFORE THE FIRST (.
!
      chr = index(list, 'RECODE') + 6
      beg = chr
100   chr = chr + 1
      If (chr>len+1) Go To 180
      If (list(chr:chr)/='(') Go To 100
      end = chr - 1
!
!       STICK THE VARIABLE NAME IN STRING THISV, AND CHECK TO MAKE
!       SURE IT WAS DEFINED IN THE DATA LIST COMMAND.
!       IF NOT, WRITE AN ERROR MESSAGE.
!
      thisv = list(beg:end)
      varnum = 0
      Do i = 1, 30
        If (thisv==varnam(i)) varnum = i
      End Do
!
!       START FINDING RECODES. INCREMENT RECODE NUMBER
!       AND TRANSFORMATION NUMBER.
!
110   recnum = recnum + 1
      tranum = tranum + 1
!
!       NO MORE THAN 200 RECODES ALLOWED.
!
      If (recnum>200) Then
        Call error(501)
      End If
!
!       FILL IN THE RECODE VARIABLE NUMBER AND SET RANGES TO -9999999
!       TEMPORARILY. FILL IN THE TRANSFORMATION TYPE AND RECODE NUMBER
!       IN THE TRANSFORMATION TABLE.
!
      recode(recnum, 1) = varnum
      recode(recnum, 3) = -99999999
      recode(recnum, 4) = -99999999
      tratab(tranum, 1) = 2
      tratab(tranum, 2) = recnum
!
!       NOW MOVE TO THE COLUMN AFTER THE (. THIS IS WHERE THE FIRST
!       VALUE SHOULD APPEAR.
!
120   chr = chr + 1
      If (chr>len+1) Go To 180
      If (list(chr:chr)=='(') Go To 120
      beg = chr
!
!       NOW LOOK FOR A CHARACTER THAT WILL INDICATE THE COLUMN
!       RIGHT AFTER THE END OF THE VALUE.
!
130   chr = chr + 1
      If (chr>len+1) Go To 180
      If (list(chr:chr)/=',' .And. list(chr:chr)/='-' .And. &
        list(chr:chr)/='=' .And. list(chr:chr)/='T' .And. list(chr:chr)/='E') &
        Go To 130
!
!       IF THE USER TRIED TO USE AN "ELSE" OPERATOR,
!       WRITE AN ERROR MESSAGE.
!
      If (list(chr:chr)=='E') Then
        If (list(chr:chr+3)=='ELSE') Then
          Call error(504)
        Else
          Call error(506)
        End If
      End If
!
!       CONVERT THE VALUE TO AN INTEGER WITH FUNCTION SELECT,
!       AND STICK IT IN THE APPROPRIATE POSITION OF ARRAY RECODE.
!
      end = chr - 1
      recode(recnum, 2) = select(list, beg, end)
!
!       IF WE ARE DONE WITH THE VALUES (AS INDICATED BY THE
!       PRESENCE OF AN EQUALS SIGN), JUMP TO 40 TO LOOK FOR NEW
!       VALUE.
!
      If (list(chr:chr)=='=') Go To 160
!
!       IF THIS IS A "THRU", FIND THE "U'.
!
      If (list(chr:chr)=='T') Then
140     chr = chr + 1
        If (chr>len+1) Go To 180
        If (list(chr:chr)/='U') Go To 140
      End If
!
!       IF THE COMMAND INDICATES A RANGE OF VALUES, FIND THE
!       SECOND VALUE. STICK IT IN THE APPROPRIATE POSITION OF
!       RECODE.  IF THE END OF THE VALUE IS DEFINED BY AN EQUALS
!       SIGN, JUMP TO 40
!
      beg = chr + 1
      If (list(chr:chr)=='U' .Or. list(chr:chr)=='-') Then
150     chr = chr + 1
        If (chr>len+1) Go To 180
        If (list(chr:chr)/=',' .And. list(chr:chr)/='=') Go To 150
        end = chr - 1
        recode(recnum, 3) = select(list, beg, end)
        If (list(chr:chr)=='=') Go To 160
      End If
!
!       SINCE WE HAVE NOT YET ENCOUNTERED AN EQUALS SIGN, GO LOOK
!       FOR MORE VALUES.
!
      Go To 110
160   Continue
!
!       NOW FIND THE NEW VALUE. IF WE ENCOUNTER A "(" OR A "="
!       BEFORE WE ENCOUNTER A ")", THERE IS A SYNTAX ERROR;
!       WRITE A MESSAGE.
!
      beg = chr + 1
170   chr = chr + 1
      If (chr>len+1) Go To 180
      If (list(chr:chr)=='(' .Or. list(chr:chr)=='=') Then
        Call error(502)
      End If
!
!       STICK THE NEW VALUE IN ITS APPROPRIATE SPOT.
!
      If (list(chr:chr)/=')') Go To 170
      end = chr - 1
      recode(recnum, 4) = select(list, beg, end)
!
!       WE ARE DONE WITH THIS SET OF PARENTHESES; GO LOOK TO
!       SEE IF THERE ARE MORE, AND PROCESS THEM IF NECESSARY.
!
      Go To 110
180   recnum = recnum - 1
      tranum = tranum - 1
!
!       NOW FILL IN THE CELLS OF RECODE THAT HAVE NOT YET
!       BEEN ASSIGNED.  IF A MAXIMIMUM OF A RANGE WAS NOT
!       ASSIGNED, IT IS THE SAME AS THE MINIMUM.
!       IF A NEW VALUE WAS NOT ASSIGNED, IT IS THE SAME AS
!       THE LAST NEW VALUE ASSIGNED WITHIN ANY GIVEN SET OF
!       PARENTHESES.
!
      Do i = 1, recnum
        p = recnum - i + 1
        If (recode(i,3)==-99999999) recode(i, 3) = recode(i, 2)
        If (recode(p,4)==-99999999) recode(p, 4) = recode(p+1, 4)
      End Do
!
!       CHECK FOR RECODE ERRORS:
!       --THE LOWER BOUND OF A RANGE OF VALUES MAY NOT BE
!         GREATER THAN THE UPPER BOUND
!       --ALL VARIABLES MUST HAVE APPEARED IN DATA LIST
!       --THE STATEMENT MUST END WITH A ")."
!
      Do i = 1, recnum
        If (recode(i,2)>recode(i,3)) Then
          Call error(503)
        End If
        If (recode(i,1)==0) Then
          Call errnam(4, thisv)
        End If
      End Do
      If (list(len:len)/=')') Then
        Call error(505)
      End If
      Return
    End Subroutine



    Subroutine stddat(list, len, varnam, vardes, varnum, numbrk, brkloc, &
      numstd, stdloc, stpop, stcat, form, wrex)
!
!       THIS SUBROUTINE INTERPRETS THE STANDARDIZE COMMAND.
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *10 thisv
      Character *500 list
      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension brkloc(5)
      Dimension stdloc(5)
!
!       INITIALIZE NUMBER OF BREAKDOWN VARIABLES AND NUMBER
!       OF STANDARDIZE VARIABLES
!
      numbrk = 0
      numstd = 0
!
!       LOOK FOR "BREAKDOWN=" OR "BREAK=".  IF THERE IS SOMETHING
!       WRONG WITH THE SYNTAX HERE, WRITE AN ERROR MESSAGE.
!
      chr = index(list, 'BREAKDOWN=') + 10
      If (chr==10) Then
        chr = index(list, 'BREAK=') + 6
        If (chr==6) Then
          If (index(list,'BREAK')>0) Then
            Call error(311)
          Else
            Call error(312)
          End If
        End If
      End If
!
!       THE LIST OF BREAKDOWN VARIABLES SHOULD BEGIN AT THIS COLUMN.
!
      beg = chr
!
!       LOOK FOR A COMMA OR A SLASH, WHICH WILL INDICATE THE END
!       OF THE VARIABLE NAME.
!
100   chr = chr + 1
      If (chr>len+1) Go To 120
      If (list(chr:chr)/=',' .And. list(chr:chr)/='/') Go To 100
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER
!
      brknum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) brknum = i
      End Do
!
!       IF THE VARIABLE WAS NOT FOUND, WRITE AN ERROR MESSAGE.
!
      If (brknum==0) Then
        Call errnam(2, thisv)
      End If
!
!       IF THE VARIABLE WAS PREVIOUSLY SPECIFIED AS THE WEIGHT
!       VARIABLE, WRITE AN ERROR MESSAGE.
!
      If (vardes(brknum,5)==7) Then
        Call error(301)
      End If
!
!       INCREMENT NUMBER OF BREAKDOWN VARIABLES.
!       IF MORE THAN FIVE SPECIFIED, WRITE ERROR MESSAGE.
!
      numbrk = numbrk + 1
      If (numbrk>5) Then
        Call error(302)
      End If
!
!       CHECK TO ENSURE THAT THE VARIABLE WAS MENTIONED IN
!       THE PRECEDING MAKETAB COMMAND.
!
      mloc = vardes(brknum, 5) - 1
      If (mloc<1 .Or. mloc>5) Then
        Call error(303)
      End If
!
!       STICK THE VARIABLE NUMBER IN THE APPROPRIATE POSITION
!       OF ARRAY BRKLOC
!
      brkloc(numbrk) = brknum
      beg = chr + 1
!
!       IF THIS IS NOT THE END OF THE BREAKDOWN VARIABLES,
!       GO BACK AND FIND ANOTHER ONE.
!
      If (list(chr:chr)/='/') Go To 100
!
!       NOW LOOK FOR THE CONTROL VARIABLES, USING THE SAME LOGIC.
!
!       LOOK FOR "CONTROL=" OR "CONT=".  IF THERE IS SOMETHING
!       WRONG WITH THE SYNTAX HERE, WRITE AN ERROR MESSAGE.
!
      chr = index(list, 'CONTROL=') + 8
      If (chr==8) Then
        chr = index(list, 'CONT=') + 5
        If (chr==5) Then
          If (index(list,'CONT')>0) Then
            Call error(311)
          Else
            Call error(313)
          End If
        End If
      End If
!
!       THE LIST OF CONTROL VARIABLES SHOULD BEGIN AT THIS COLUMN.
!
      beg = chr
!
!       LOOK FOR A COMMA OR A SLASH, WHICH WILL INDICATE THE END
!       OF THE VARIABLE NAME.
!
110   chr = chr + 1
      If (list(chr:chr)/=',' .And. list(chr:chr)/='/' .And. chr<len+1) &
        Go To 110
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER
!
      stdnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) stdnum = i
      End Do
!
!       IF THE VARIABLE WAS NOT FOUND, WRITE AN ERROR MESSAGE.
!
      If (stdnum==0) Then
        Call errnam(2, thisv)
      End If
!
!       IF THE VARIABLE WAS PREVIOUSLY SPECIFIED AS THE WEIGHT
!       VARIABLE, WRITE AN ERROR MESSAGE.
!
      If (vardes(stdnum,5)==7) Then
        Call error(301)
      End If
!
!       INCREMENT NUMBER OF CONTROL VARIABLES.
!       IF MORE THAN FIVE SPECIFIED, WRITE ERROR MESSAGE.
!
      numstd = numstd + 1
      If (numstd>5) Then
        Call error(304)
      End If
!
!       IF FIVE CONTROL VARIABLES WERE SPECIFIED, THE BREAKDOWN
!       VARIABLES MUST BE IDENTICAL SO THAT ONLY FOUR
!       CONTROL VARIABLES ARE USED SIMULTANEOUSLY.  CHECK TO SEE IF
!       THIS IS THE CASE.
!
      If (numstd==5 .And. numbrk/=5) Call error(304)
!
!       CHECK TO ENSURE THAT THE VARIABLE WAS MENTIONED IN
!       THE PRECEDING MAKETAB COMMAND.
!
      mloc = vardes(stdnum, 5) - 1
      If (mloc<1 .Or. mloc>5) Then
        Call error(303)
      End If
!
!       STICK THE VARIABLE NUMBER IN THE APPROPRIATE POSITION
!       OF ARRAY STDLOC
!
      stdloc(numstd) = stdnum
!
!       IF THIS IS NOT THE END OF THE CONTROL VARIABLES,
!       GO BACK AND FIND ANOTHER ONE.
!
      beg = chr + 1
      If (list(chr:chr)/='/' .And. chr<len+1) Go To 110
!
!       NOW INTERPRET THE "STANDARD=" SUBCOMMAND, IF PRESENT.
!
      stcat = 0
      stpop = 0
      form = 1
      chr = index(list, 'STANDARD=AV')
      If (chr>0) stpop = 1
      If (stpop==0) Then
        chr = index(list, 'STANDARD=TO')
        If (chr>0) stpop = 2
        If (stpop==0) Then
          chr = index(list, 'STANDARD=CA')
          If (chr>0) Then
!
!       IF THE STANDARD IS A CATEGORY OF THE BREAKDOWN VARIABLE,
!       FIND THE CATEGORY, WHICH SHOULD BE IN PARENTHESES.
!       ONLY ONE BREAKDOWN VARIABLE IS ALLOWED IF THIS OPTION
!       IS CHOSEN; IF THERE IS MORE THAN ONE, WRITE ERROR MESSAGE.
!
            If (numbrk>1) Then
              Call error(314)
            End If
            stpop = 3
            beg = index(list, '(') + 1
            end = index(list, ')') - 1
            If (end<beg .Or. beg<=chr) Then
              Call error(305)
            End If
            stcat = select(list, beg, end)
          End If
        End If
      End If
!
!       IF A STANDARD SUBCOMMAND IS THERE BUT HAS NOT BEEN
!       INTERPRETED, WRITE ERROR. OTHERWISE, ASSIGN DEFAULT.
!
      If (stpop==0) Then
        chr = index(list, 'STANDARD=')
        If (chr>0) Then
          Call error(306)
        Else
          stpop = 2
        End If
      End If
!
!       NOW INTERPRET "FORMAT=" SUBCOMMAND.
!
      chr = index(list, 'FORMAT=PE')
      If (chr>0) Then
        form = 2
      Else
        chr = index(list, 'FORMAT=DE')
        If (chr>0) form = 3
      End If
!
!       NOW CHECK TO SEE IF A "/WRITE EXCLUDED" SUBCOMMAND IS PRESENT.
!
      wrex = index(list, '/WRI')
120   Continue
!
!       CHECK TO MAKE SURE A STANDARDIZE VARIABLE WAS SPECIFIED.
!
      If (numstd<1 .And. numbrk>=1) Then
        Call error(315)
      End If
!
!       CHECK TO MAKE SURE A BREAKDOWN VARIABLE WAS SPECIFIED.
!
      If (numbrk<1) Then
        Call error(307)
      End If
      Return
    End Subroutine


    Subroutine weight(list, len, varnam, vardes, varnum, wtdig)
!
!       THIS ROUTINE INTERPRETS THE WEIGHT COMMAND.
!
      Implicit Integer (A-Z)
      Character *10 varnam, thisv
      Character *500 list
      Dimension vardes(30, 5)
      Dimension varnam(30)

      len = 80
      wtnum = 0
!
!       CHECK TO SEE IF " BY " WAS SPECIFIED.
!       IF NOT, ASSUME THAT THE VARIABLE NAME BEGINS RIGHT AFTER
!       THE "WEIGHT" KEYWORD.
!
      chr = index(list, ' BY ') + 3
      If (chr==3) Then
        chr = index(list, 'WEIGHT ') + 6
        If (chr==6) Then
          Call error(702)
        End If
      End If
!
!       SKIP ANY LEADING SPACES.
!
100   chr = chr + 1
      If (chr>len+1) Go To 120
      If (list(chr:chr)==' ') Go To 100
      beg = chr
!
!       LOOK FOR A TRAILING SPACE, PERIOD OR COMMA TO DEFINE END
!       OF VARIABLE NAME.
!
110   chr = chr + 1
      If (chr>len+1) Go To 120
      If (list(chr:chr)/=' ' .And. list(chr:chr)/='.' .And. &
        list(chr:chr)/='(') Go To 110
      end = chr - 1
      thisv = list(beg:end)
!
!       LOOK FOR VARIABLE NAME
!
      Do i = 1, varnum
        If (thisv==varnam(i)) wtnum = i
      End Do
!
!       IF NO NAME FOUND, WRITE ERROR MESSAGE.
!
120   If (wtnum==0) Then
        Call errnam(6, thisv)
      Else
        vardes(wtnum, 5) = 7
      End If
!
!       NOW LOOK FOR DECIMAL PLACES, IF THEY EXIST.
!
      wtdig = 0
      chr = index(list, '(')
      If (chr>0) Then
        beg = chr + 1
        end = index(list, ')') - 1
!
!       CHECK SYNTAX
!
        If (beg>end .Or. end<1) Then
          Call error(702)
        End If
!
!       PUT NUMBER OF DIGITS IN WTDIG
!
        wtdig = select(list, beg, end)
        If (wtdig>5) wtdig = 5
      End If
      Return
    End Subroutine


    Block Data
!
!       THIS INITIALIZES THE COMMON BLOCK
!
      Common /tabs/table, maxvars
      Dimension table(640000, 2)
      Dimension maxvars(5)
      Data table/1280000*0.0/, maxvars/5*0/
    End Block Data



    Subroutine make(list, len, varnam, vardes, varnum, depdig, wrtab)
!
!       THIS ROUTINE INTERPRETS THE MAKETAB COMMAND.
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *10 thisv
      Character *500 list
      Dimension vardes(30, 5)
      Dimension varnam(30)
!
!       INITIALIZE ALL VARIABLE TYPES TO "OTHER," IF THEY HAVE NOT
!       ALREADY BEEN SET TO BE A WEIGHT VARIABLE.
!
      Do i = 1, varnum
        If (vardes(i,5)/=7) vardes(i, 5) = 8
      End Do
      numind = 0
!
!       LOOK FOR DEPENDENT VARIABLE, WHICH SHOULD BEGIN IMMEDIATELY
!       AFTER THE EQUALS SIGN.
!
      chr = index(list, '=') + 1
      beg = chr
!
!       THE LAST CHARACHTER OF THE DEPENDENT VARIABLE SHOULD BE
!       FOLLOWED BY A "/" OR A "(".  LOCATE IT, AND GRAB VARIABLE
!       NAME.
!
100   chr = chr + 1
      If (chr>len+1) Go To 160
      If (list(chr:chr)/='/' .And. list(chr:chr)/='(') Go To 100
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER.  IF NOT FOUND, WRITE ERROR MESSAGE.
!
      depnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) depnum = i
      End Do
      If (depnum==0) Then
        Call errnam(1, thisv)
      End If
!
!       DEPENDENT VARIABLE CANNOT BE WEIGHT VARIABLE.  IF IT IS,
!       WRITE ERROR MESSAGE.
!
      If (vardes(depnum,5)==7) Then
        Call error(201)
      End If
!
!       NOW IDENTIFY DEPENDENT VARIABLE IN VARDES.
!
      vardes(depnum, 5) = 1
!
!       LOOK FOR NEXT EQUALS SIGN.  THE FIRST INDEPENDENT VARIABLE
!       SHOULD FOLLOW IT IMMEDIATELY.
!
110   chr = chr + 1
      If (chr>len+1) Go To 160
      If (list(chr:chr)/='=') Go To 110
!
!       BEFORE PROCEEDING, CHECK TO SEE IF THERE ARE DIGITS TO
!       THE RIGHT OF THE DECIMAL PLACE ON THE DEPENDENT VARIABLE
!       SPECIFIED.  IF SO, SEE HOW MANY AND PUT THAT NUMBER IN
!       DEPDIG.
!
      depdig = 0
      nchr = index(list, '(')
      If (nchr<chr .And. nchr>0) Then
        beg = nchr + 1
        end = index(list, ')') - 1
        If (end>=beg) depdig = select(list, beg, end)
      End If
!
!       NOW WE ARE READY TO LOCATE THE DEPENDENT VARIABLES.
!
120   beg = chr + 1
130   chr = chr + 1
!
!       IF WE ENCOUNTER A SLASH OR THE END OF THE LINE, WE ARE DONE.
!
      If (chr>len+1 .Or. list(chr:chr)=='/') Go To 160
!
!       A ")" OR A "," AT THIS POINT INDICATES A SYNTAX ERROR.
!
      If (list(chr:chr)==')' .Or. list(chr:chr)==',') Then
        Call error(202)
      End If
!
!       THE LAST COLUMN OF THE INDEPENDENT VARIABLE NAME WILL
!       BE PRECEDED BY A "("
!
      If (list(chr:chr)/='(') Go To 130
      end = chr - 1
!
!       SEARCH FOR INDEPENDENT VARIABLE.  IF IT HAS NOT BEEN
!       PREVIOUSLY DEFINED, WRITE AN ERROR MESSAGE.
!
      thisv = list(beg:end)
      indnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) indnum = i
      End Do
      If (indnum==0) Then
        Call errnam(1, thisv)
      End If
!
!       INCREMENT NUMBER OF INDEPENDENT VARIABLES.  IF THERE ARE
!       MORE THAN FIVE, WRITE AN ERROR MESSAGE.
!
      numind = numind + 1
      If (numind>5) Then
        Call error(203)
      End If
!
!       IF THE VARIABLE HAS BEEN PREVIOUSLY SPECIFIED AS THE
!       DEPENDENT VARIABLE OR THE WEIGHT VARIABLE, WRITE AN
!       ERROR MESSAGE.
!
      If (vardes(indnum,5)==7 .Or. vardes(indnum,5)==1) Then
        Call error(201)
      End If
!
!       IDENTIFY THIS VARIABLE AS AN INDEPENDENT VARIABLE IN
!       ARRAY VARDES.
!
      vardes(indnum, 5) = numind + 1
!
!         NOW GET MINIMUM AND MAXIMUM.  THE NEXT COLUMN SHOULD BE
!         THE START OF THE MINIMUM.
!
      beg = chr + 1
140   chr = chr + 1
      If (chr>len+1) Go To 160
!
!         AN ASTERISK IS ALLOWED IN LIEU OF THE MINIMUM AND THE
!         MAXIMUM ONLY IF THE MIN AND MAX HAVE BEEN PREVIOUSLY
!         DEFINED BY A COMBINE COMMAND. OTHERWISE, WRITE AN ERROR
!         MESSAGE.
!
      If (list(chr:chr)=='*') Then
        If (vardes(indnum,3)/=0 .Or. vardes(indnum,4)/=0) Then
          chr = chr + 1
          Go To 120
        Else
          Call errnam(251, varnam(indnum))
        End If
      End If
!
!         CHECK SYNTAX.
!
      If (list(chr:chr)=='(' .Or. list(chr:chr)==')') Then
        Call error(202)
      End If
!
!         A COMMA INDICATES THE END OF THE MINIMUM
!
      If (list(chr:chr)/=',') Go To 140
      end = chr - 1
!
!         STICK MINIMUM IN APPROPRIATE SLOT OF VARDES
!
      vardes(indnum, 3) = select(list, beg, end)
!
!         NOW FIND MAXIMUM
!
      beg = chr + 1
150   chr = chr + 1
      If (chr>len+1) Go To 160
!
!         CHECK FOR SYNTAX ERRORS
!
      If (list(chr:chr)=='(' .Or. list(chr:chr)==',') Then
        Call error(202)
      End If
!
!         ")" INDICATES THE END OF THE MAXIMUM
!
      If (list(chr:chr)/=')') Go To 150
!
!         STICK MAXIMUM IN APPROPRIATE SLOT OF VARDES
!
      end = chr - 1
      vardes(indnum, 4) = select(list, beg, end)
!
!         NO GO LOOK FOR MORE INDEPENDENT VARIABLES
!
      Go To 120
160   Continue
!
!         CHECK FOR A "WRITE TABLE" SUBCOMMAND
!
      wrtab = index(list, '/WRI')
!
!         CHECK FOR ERRORS:
!           MISSING DEPENDENT SUBCOMMAND
!           MISSING VARIABLES SUBCOMMAND
!           MISSING EQUALS SIGN
!           TOO FEW INDEPENDENTS
!           WRONG END OF STATEMENT
!
      test1 = index(list, 'DEPENDENT=')
      If (test1==0) test1 = index(list, 'DEP=')
      test2 = index(list, 'DEP')
      test3 = index(list, 'VARIABLES=')
      If (test3==0) test3 = index(list, 'VARIABLE=')
      If (test3==0) test3 = index(list, 'VARS=')
      If (test3==0) test3 = index(list, 'VAR=')
      test4 = index(list, 'VAR')
      If (test2==0) Then
        Call error(205)
      Else If (test4==0) Then
        Call error(206)
      Else If (test4<=test2) Then
        Call error(207)
      Else If (test1==0 .Or. test3==0) Then
        Call error(208)
      Else If (numind<1) Then
        Call error(204)
      Else If (list(len:len)/=')' .And. wrtab==0) Then
        Call error(209)
      End If
!
!         CHECK TO SEE IF MINIMUM IS LESS THAN MAXIMUM
!
      Do i = 1, varnum
        If (vardes(i,5)>=2 .And. vardes(i,5)<=5) Then
          If (vardes(i,3)>vardes(i,4)) Then
            Call error(210)
          End If
        End If
      End Do
      Return
    End Subroutine



    Function ind(i, j, k, l, m)
!
!      THIS FUNCTION CALCULATES THE INDEX FOR TABULATING THE
!      VARIABLES
!
      Common /tabs/table(640000, 2), maxvars(5)

      ind = i + j*maxvars(1) + k*maxvars(1)*maxvars(2) + &
        l*maxvars(1)*maxvars(2)*maxvars(3) + m*maxvars(1)*maxvars(2)*maxvars(3 &
        )*maxvars(4) - (maxvars(1)+maxvars(1)*maxvars(2)+maxvars(1)*maxvars(2) &
        *maxvars(3)+maxvars(1)*maxvars(2)*maxvars(3)*maxvars(4))
      If (ind>640000) Call error(703)
      Return
    End Function



    Subroutine seldat(list, len, seldes, varnam, selnum, tratab, tranum)
!
!      THIS ROUTINE INTERPRETS THE SELECT IF COMMAND.  THE RESULTS
!      ARE STORED IN ARRAY SELDES (SELECT IF DESCRIPTIONS)
!
!      SELDES: SELECT IF DESCRIPTIONS
!              ROWS  - SELECT IF NUMBER
!              COL 1 - VARIABLE NUMBER
!              COL 2 - VALUE
!              COL 3 - OPERATOR
!                      1 - EQ
!                      2 - GT
!                      3 - LT
!                      4 - GE
!                      5 - LE
!              COL 4 - CONSTANT FOR COMPARISON
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *500 list
      Character *10 thisv

      Dimension seldes(50, 3)
      Dimension varnam(30)
      Dimension tratab(300, 2)
!
!       FIND FIRST VARIABLE NAME; IT SHOULD APPEAR AFTER THE "("
!
      chr = index(list, '(') + 1
      beg = chr
!
!       FIND THE "~" DELIMITER, WHICH INDICATES AN OPERATOR
!       BEGINS AND THE VARIABLE NAME ENDS.
!
! ****** WARNING: THE CHARACTER ~ MAY NOT BE AVAILABLE ON SOME
!                 SYSTEMS; SUBSTITUTE SOMETHING ELSE.
!
100   chr = chr + 1
      If (chr>len+1) Go To 120
      If (list(chr:chr)/='~') Go To 100
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER;
!
      varnum = 0
      Do i = 1, 30
        If (thisv==varnam(i)) varnum = i
      End Do
!
!       FIND THE OPERATOR
!
      op = 0
      chr = index(list, '~EQ~') + 4
      If (chr>4) Then
        op = 1
      Else
        chr = index(list, '~GT~') + 4
        If (chr>4) Then
          op = 2
        Else
          chr = index(list, '~LT~') + 4
          If (chr>4) Then
            op = 3
          Else
            chr = index(list, '~GE~') + 4
            If (chr>4) Then
              op = 4
            Else
              chr = index(list, '~LE~') + 4
              If (chr>4) Then
                op = 5
              Else
                chr = index(list, '~NE~') + 4
                If (chr>4) Then
                  op = 6
                Else
                  Call error(601)
                End If
              End If
            End If
          End If
        End If
      End If
!
!     INCREMENT SELECT IF NUMBER AND TRANSFORMATION NUMBER
!
      selnum = selnum + 1
      tranum = tranum + 1
!
!     A MAXIMUM OF 50 SELECT IF'S ARE ALLOWED.
!
      If (selnum>50) Then
        Call error(602)
      End If
!
!     STICK THE INFORMATION GATHERED SO FAR IN THE APPROPRIATE
!     SLOTS OF SELDES.
!
      seldes(selnum, 1) = varnum
      seldes(selnum, 2) = op
      tratab(tranum, 1) = 1
      tratab(tranum, 2) = selnum
!
!       NOW FIND THE CONSTANT; IT SHOULD BEGIN RIGHT AFTER THE
!       OPERATOR
!
      beg = chr
!
!       THE CONSTANT ENDS RIGHT BEFORE THE ")"
!
110   chr = chr + 1
      If (chr>len+1) Go To 120
      If (list(chr:chr)/=')') Go To 110
!
!       CONVERT IT TO AN INTEGER WITH SELECT AND STICK IT IN THE
!       APPROPRIATE SLOT OF SELDES
!
      end = chr - 1
      seldes(selnum, 3) = select(list, beg, end)
!
!       CHECK FOR SELECT ERRORS
!
120   Continue
!
!       IF THE VARIABLE WAS NOT FOUND,
!       WRITE AN ERROR MESSAGE.
!
      If (seldes(selnum,1)==0) Then
        Call errnam(5, thisv)
      End If
!
!      MAKE SURE STATEMENT ENDS WITH A ")".
!
      If (list(len:len)/=')') Then
        Call error(603)
      End If
      Return
    End Subroutine


    Subroutine decomp(list, len, varnam, vardes, varnum, poploc, pop1d, pop2d, &
      numfac, facloc, wrex)
!
!       THIS SUBROUTINE INTERPRETS THE DECOMPOSE COMMAND
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *10 thisv
      Character *500 list
      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension facloc(5)
!
!       INITIALIZE NUMBER OF FACTORS
!
      numfac = 0
!
!       LOCATE "COMPARE" SUBCOMMAND; WRITE ERROR MESSAGE IF
!       IT IS MISSING OR INCOMPLETE
!
      chr = index(list, 'COMPARE=') + 8
      If (chr==8) Then
        If (index(list,'COMPARE')>0) Then
          Call error(408)
        Else
          Call error(410)
        End If
      End If
!
!       FIND COMPARISON VARIABLE; IT SHOULD BEGIN RIGHT AFTER
!       THE EQUALS SIGN.
!
      beg = chr
100   chr = chr + 1
      If (chr>len+1) Go To 140
!
!       THE VARIABLE WILL END RIGHT BEFORE THE "("
!
      If (list(chr:chr)/='(') Go To 100
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER; WRITE ERROR IF VARIABLE NOT
!       FOUND
!
      popnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) popnum = i
      End Do
      If (popnum==0) Then
        Call errnam(3, thisv)
      End If
!
!       WRITE ERROR MESSAGE IF THE COMPARISON VARIABLE IS ALSO
!       THE WEIGHT VARIABLE
!
      If (vardes(popnum,5)==7) Then
        Call error(401)
      End If
!
!       WRITE ERROR MESSAGE IF COMPARISON VARIABLE WAS NOT DECLARED
!       AS A VARIABLE IN THE PRECEDING MAKETAB COMMAND.
!
      mloc = vardes(popnum, 5) - 1
      If (mloc<1 .Or. mloc>5) Then
        Call error(402)
      End If
!
!       STORE COMPARISON VARIABLE IN POPLOC AD SEARCH FOR COMPARISON
!       CATEGORIES, WHICH SHOULD BEGIN IN THE NEXT COLUMN
!
      poploc = popnum
      pop1d = 9999
      pop2d = 9999
      beg = chr + 1
110   chr = chr + 1
      If (chr>len+1) Go To 140
!
!       A COLON SHOULD SEPARATE THE TWO COMPARISON CATEGORIES
!
      If (list(chr:chr)/=':') Go To 110
      end = chr - 1
!
!       CONVERT CATEGORY 1 TO AN INTEGER USING FUNCTION SELECT
!
      pop1d = select(list, beg, end)
!
!       LOCATE SECOND COMPARISON CATEGORY, WHICH SHOULD BE FOLLOWED
!       BY A ")"
!
      beg = chr + 1
120   chr = chr + 1
      If (chr>len+1) Go To 140
      If (list(chr:chr)/=')') Go To 120
      end = chr - 1
!
!       CONVERT CATEGORY 2 TO AN INTEGER USING FUNCTION SELECT
!
      pop2d = select(list, beg, end)
!
!       NOW LOCATE "FACTORS=" SUBCOMMAND; WRITE ERROR MESSAGE IF IT
!       IS MISSING OR INCOMPLETE
!

      chr = index(list, 'FACTORS=') + 8
      If (chr==8) Then
        chr = index(list, 'FACTOR=') + 7
        If (chr==7) Then
          If (index(list,'FACTOR')>0) Then
            Call error(408)
          Else
            Call error(409)
          End If
        End If
      End If
!
!       THE FIRST FACTOR SHOULD BEGIN IMMEDIATELY FOLLOWING THE
!       EQUALS SIGN, AND END RIGHT BEFORE A COMMA, SLASH, OR THE
!       END OF THE COMMAND
!
      beg = chr

130   chr = chr + 1
      If (list(chr:chr)/=',' .And. chr<len+1 .And. list(chr:chr)/='/') &
        Go To 130
      end = chr - 1
      thisv = list(beg:end)
!
!       SEARCH FOR VARIABLE NUMBER; IF IT IS NOT FOUND, WRITE AN
!       ERROR MESSAGE.
!
      facnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) facnum = i
      End Do
      If (facnum==0) Then
        Call errnam(3, thisv)
      End If
!
!       THE WEIGHT VARIABLE IS NOT ALLOWED TO BE A FACTOR
!
      If (vardes(facnum,5)==7) Then
        Call error(401)
      End If
!
!       NO MORE THAN FOUR FACTORS ARE ALLOWED
!
      numfac = numfac + 1
      If (numfac>4) Then
        Call error(403)
      End If
!
!       ALL FACTORS MUST HAVE BEEN DECLARED AS VARIABLES IN THE
!       PRECEDING MAKETAB COMMAND
!
      facloc(numfac) = facnum
      mloc = vardes(facnum, 5) - 1
      If (mloc<1 .Or. mloc>5) Then
        Call error(402)
      End If
!
!       IF WE HAVE NOT ENCOUNTERED THE END OF THE LIST OF FACTORS,
!       AS INDICATED BY A SLASH OR THE END OF THE COMMAND, GO BACK
!       AND LOCATE MORE FACTORS.
!
      beg = chr + 1
      If (chr<len+1 .And. list(chr:chr)/='/') Go To 130
!
!       SEE IF THERE IS A "/WRITE EXCLUDED" SUBCOMMAND.
!
      wrex = index(list, '/WRI')
140   Continue
!
!       NOW CHECK TO MAKE SURE THAT THE COMPARISON CATEGORIES WERE
!       DEFINED
!
      If (pop1d==9999 .Or. pop2d==9999) Then
        Call error(404)
      End If
      Return
    End Subroutine




    Subroutine guts(n1t, n2t, num1, num2, tot1, tot2, comb, rate, effect)
!
!      THIS ROUTINE IS THE GUTS OF THE DECOMPOSITION
!
      Implicit Integer (A-Z)
      Real a1(4), a2(6), a3(4), b1(4), b2(6), b3(4)
      Real effect(4), tot1, tot2, comb, rate
      Real num1, num2
      Real table, n1, n2, av, n1t, n2t, t1, t2
      Dimension table(640000, 2)
      Dimension maxvars(5)
      Character *500 order(24)
      Common /tabs/table, maxvars
!cc       INQUIRE(UNIT=6,NAME=SCRNAM)

!      THERE ARE 24 EQUATIONS NECESSARY FOR DAS GUPTA'S METHODS,
!      BUT THEY CAN BE EXPRESSED AS A SINGLE EQUATION BY SUBSTITUTING
!      TERMS.  THE ARRAY ORDER CONTAINS THE INFORMATION NECESSARY TO
!      SUBSTITUTE THE TERMS IN THE CORRECT ORDER.
!
      order(1) = '125314'
      order(2) = '136314'
      order(3) = '144224'
      order(4) = '125224'
      order(5) = '136134'
      order(6) = '144134'
      order(7) = '215613'
      order(8) = '233613'
      order(9) = '242423'
      order(10) = '215423'
      order(11) = '233143'
      order(12) = '242143'
      order(13) = '316512'
      order(14) = '323512'
      order(15) = '341432'
      order(16) = '316432'
      order(17) = '323242'
      order(18) = '341242'
      order(19) = '414521'
      order(20) = '422521'
      order(21) = '431631'
      order(22) = '414631'
      order(23) = '422341'
      order(24) = '431341'
!
!      THE FOLLOWING LOOP CALCULATES THE RATE EFFECT AND THE
!      COMBINED EFFECT OF FACTORS.  THESE ARE PRETTY SIMPLE:
!      THE COMBINED EFFECT OF FACTORS IS THE SUM OF THE AVERAGE
!      OF THE RATES FOR THE TWO POPULATIONS TIMES THE DIFFERENCE
!      OF THE PROPORTIONS OF THE POPULATION IN EACH CELL.  THE
!      RATE EFFECT IS THE SUM OF THE AVERAGE OF THE PROPORTIONS
!      TIMES THE DIFFERENCE IN THE RATES OVER ALL CELLS.  THE
!      SUM OF THE RATE EFFECT AND THE EFFECT OF FACTORS SHOULD
!      BE EQUAL TO THE TOTAL POPULATION DIFFERENCE.  THIS LOOP
!      ALSO SUMS THE CASES IN EACH POPULATION AND CALCULATES THE
!      OVERALL RATE IN EACH POPULATION.
!
      Do i = 2, maxvars(2)
        Do j = 2, maxvars(3)
          Do k = 2, maxvars(4)
            Do l = 2, maxvars(5)
              n1 = table(ind(1,i,j,k,l), 2)
              n2 = table(ind(2,i,j,k,l), 2)
              If (n1>0.0 .And. n2>0.0) Then
                t1 = table(ind(1,i,j,k,l), 1)
                t2 = table(ind(2,i,j,k,l), 1)
                av = table(ind(3,i,j,k,l), 1)
                num1 = num1 + n1
                num2 = num2 + n2
                tot1 = tot1 + (t1*n1/n1t)
                tot2 = tot2 + (t2*n2/n2t)
                comb = comb + av*(n1/n1t-n2/n2t)
                rate = rate + (n1/n1t+n2/n2t)/2.0*(t1-t2)
              End If
            End Do
          End Do
        End Do
      End Do
!
!       NOW WE ARE READY FOR THE GUTS OF THE GUTS.
!       FOR EACH FACTOR, THE EQUATION MUST BE CARRIED OUT IN
!       EACH OF SIX POSSIBLE ORDERS, MAKING A TOTAL OF 24
!       PASSES THROUGH THE MAIN TABLE.
!
      ordfac = 0
      Do fact = 1, 4
        Do ord = 1, 6
!
!       SEND A MESSAGE TO THE SCREEN INFORMING THE ANALYST OF
!       OUR PROGRESS (WITH A COMPLICATED PROBLEM AND A SLOW
!       COMPUTER, THIS LOOP CAN TAKE QUITE A LONG TIME, SO THIS
!       MESSAGE JUST REASSURES THE USER THAT SOMETHING IS HAPPENING.
!
          ordfac = ordfac + 1
!
!       EXTRACT THE POINTERS FROM ARRAY ORDER.
!
          p1 = select(order(ordfac), 1, 1)
          p2 = select(order(ordfac), 2, 2)
          p3 = select(order(ordfac), 3, 3)
          p4 = select(order(ordfac), 4, 4)
          p5 = select(order(ordfac), 5, 5)
          p6 = select(order(ordfac), 6, 6)
!
!       THIS LOOP RUNS THROUGH ALL THE CELLS
!
          Do i = 2, maxvars(2)
            Do j = 2, maxvars(3)
              Do k = 2, maxvars(4)
                Do l = 2, maxvars(5)
!
!       GET THE NUMBER OF CASES IN THE CELL AND THE AVERAGE OF
!       THE RATES
!
                  n1 = table(ind(1,i,j,k,l), 2)
                  n2 = table(ind(2,i,j,k,l), 2)
                  av = table(ind(3,i,j,k,l), 1)
!
!       SKIP THE CELL IF IT IS EMPTY IN EITHER POPULATION
!
                  If (n1<=0 .Or. n2<=0) Go To 100
!
!       GET ALL THE TERMS OF THE EQUATION, WHICH ARE ALL OF THE
!       POSSIBLE MARGINALS
!
                  a3(1) = table(ind(1,i,j,k,1), 2)
                  a3(2) = table(ind(1,i,j,1,l), 2)
                  a3(3) = table(ind(1,i,1,k,l), 2)
                  a3(4) = table(ind(1,1,j,k,l), 2)
                  a2(1) = table(ind(1,1,1,k,l), 2)
                  a2(2) = table(ind(1,1,j,1,l), 2)
                  a2(3) = table(ind(1,1,j,k,1), 2)
                  a2(4) = table(ind(1,i,1,1,l), 2)
                  a2(5) = table(ind(1,i,j,1,1), 2)
                  a2(6) = table(ind(1,i,1,k,1), 2)
                  a1(1) = table(ind(1,i,1,1,1), 2)
                  a1(2) = table(ind(1,1,j,1,1), 2)
                  a1(3) = table(ind(1,1,1,k,1), 2)
                  a1(4) = table(ind(1,1,1,1,l), 2)
                  b3(1) = table(ind(2,i,j,k,1), 2)
                  b3(2) = table(ind(2,i,j,1,l), 2)
                  b3(3) = table(ind(2,i,1,k,l), 2)
                  b3(4) = table(ind(2,1,j,k,l), 2)
                  b2(1) = table(ind(2,1,1,k,l), 2)
                  b2(2) = table(ind(2,1,j,1,l), 2)
                  b2(3) = table(ind(2,1,j,k,1), 2)
                  b2(4) = table(ind(2,i,1,1,l), 2)
                  b2(5) = table(ind(2,i,j,1,1), 2)
                  b2(6) = table(ind(2,i,1,k,1), 2)
                  b1(1) = table(ind(2,i,1,1,1), 2)
                  b1(2) = table(ind(2,1,j,1,1), 2)
                  b1(3) = table(ind(2,1,1,k,1), 2)
                  b1(4) = table(ind(2,1,1,1,l), 2)
!
!       THIS IS IT! (SEE DAS GUPTA, 1978 FOR AN EXPLANATION)
!
                  effect(fact) = effect(fact) + (av*(n1/a1(p1)+n2/b1(p1))/2.0* &
                    (a1(p1)/n1t-b1(p1)/n2t)) + (av*(n1/a2(p3)+n2/b2(p3))/2.0*( &
                    a1(p2)/n1t+b1(p2)/n2t)/2.0*(a2(p3)/a1(p2)-b2(p3)/b1(p2))) &
                    + (av*(n1/a3(p5)+n2/b3(p5))/2.0*(a1(p2)/n1t+b1(p2)/n2t)/ &
                    2.0*(a2(p4)/a1(p2)+b2(p4)/b1(p2))/2.0*(a3(p5)/a2(p4)-b3( &
                    p5)/b2(p4))) + (av*(a1(p2)/n1t+b1(p2)/n2t)/2.0*(a2(p4)/a1( &
                    p2)+b2(p4)/b1(p2))/2.0*(a3(p6)/a2(p4)+b3(p6)/b2( &
                    p4))/2.0*(n1/a3(p6)-n2/b3(p6)))
100             End Do
              End Do
            End Do
          End Do
        End Do
      End Do
      Return
    End Subroutine






    Subroutine combdat(list, len, comdes, varnam, vardes, varnum, comnum, &
      tratab, tranum, pg, ln)
!
!       THIS ROUTINE INTERPRETS THE COMBINE COMMAND.
!
      Implicit Integer (A-Z)
      Character *10 varnam
      Character *500 list
      Character *10 thisv

      Dimension vardes(30, 5)
      Dimension varnam(30)
      Dimension comdes(10, 5)
      Dimension tratab(300, 2)
!
!       COMDES: COMBINE DESCRIPTIONS
!
!                 ROWS: COMBINATION VARIABLE NUMBER
!                COL 1: VARIABLE NUMBER OF FIRST VARIABLE TO COMBINE
!                COL 2: VARIABLE NUMBER OF SECOND VARIABLE TO COMBINE
!                COL 3: VARIABLE NUMBER OF THIRD VARIABLE TO COMBINE
!                COL 4: VARIABLE NUMBER OF FOURTH VARIABLE TO COMBINE
!                COL 5: VARIABLE NUMBER OF NEW COMBINED VARIABLE
!
!       FIND FIRST VARIABLE NAME.  THIS WILL BE THE NAME OF THE NEW
!       VARIABLE TO BE CREATED.
!
      beg = index(list, 'COMBINE') + 7
      end = index(list, '=') - 1
      thisv = list(beg:end)
!
!       INCREMENT VARIABLE NUMBER, COMBINATION NUMBER, AND
!       TRANSFORMATION NUMBER.  INSERT FIRST VARIABLE NUMBER
!       IN THE ARRAY COMBDES, TRANSFORMATION TYPE IN ARRAY TRATAB,
!       ETC.
!
      numind = 0
      varnum = varnum + 1
      comnum = comnum + 1
      tranum = tranum + 1
      comdes(comnum, 5) = varnum
      tratab(tranum, 1) = 3
      tratab(tranum, 2) = comnum
      vardes(varnum, 1) = 9999
      Do i = 1, 4
        comdes(comnum, i) = 99
      End Do
      varnam(varnum) = thisv
!
!       NOW FIND THE OTHER VARIABLES, WHICH SHOULD BEGIN RIGHT
!       AFTER THE EQUALS SIGN.
!
      chr = index(list, '=')
100   beg = chr + 1
110   chr = chr + 1
      If (chr>len+1) Go To 140
!
!       CHECK FOR SYNTAX ERRORS.
!
      If (list(chr:chr)==')' .Or. list(chr:chr)==',') Then
        Call error(801)
      End If
!
!       THE END OF THE VARIABLE NAME SHOULD APPEAR RIGHT BEFORE
!       THE "("
!
      If (list(chr:chr)/='(') Go To 110
      end = chr - 1
!
!       SEARCH FOR INDEPENDENT VARIABLE.  CHECK TO MAKE SURE IT
!       HAS BEEN DEFINED IN THE PRECEDING DATA LIST COMMAND.
!
      thisv = list(beg:end)
      indnum = 0
      Do i = 1, varnum
        If (thisv==varnam(i)) indnum = i
      End Do
      If (indnum==0) Then
        Call errnam(7, thisv)
      End If
!
!       INCREMENT NUMBER OF VARIABLES FOUND.
!       ONLY FOUR COMBINE VARIABLES ARE ALLOWED.
!
      numind = numind + 1
      If (numind>4) Then
        Call error(802)
      End If
!
!       NOW GET MINIMUM AND MAXIMUM; THEY ARE IN PARENTHESES,
!       SEPARATED BY A COMMA.
!
      beg = chr + 1
120   chr = chr + 1
      If (chr>len+1) Go To 140
!
!       CHECK PUNCTUATION.
!
      If (list(chr:chr)=='(' .Or. list(chr:chr)==')') Then
        Call error(801)
      End If
      If (list(chr:chr)/=',') Go To 120
      end = chr - 1
!
!       STICK MINIMUM IN APPROPRIATE POSITION OF VARDES.
!
      vardes(indnum, 3) = select(list, beg, end)
      beg = chr + 1
130   chr = chr + 1
      If (chr>len+1) Go To 140
!
!       CHECK PUNCTUATION
!
      If (list(chr:chr)=='(' .Or. list(chr:chr)==',') Then
        Call error(801)
      End If
      If (list(chr:chr)/=')') Go To 130
      end = chr - 1
!
!       STICK MAXIMUM IN APPROPRIATE POSITION OF VARDES, AND STICK
!       VARIABLE NUMBER IN APPROPRIATE POSITION OF COMBDES.
!
      vardes(indnum, 4) = select(list, beg, end)
      comdes(comnum, numind) = indnum
!
!       ENSURE THAT MINIMUM IS LESS THAN MAXIMUM.
!
      If (vardes(indnum,3)>vardes(indnum,4)) Then
        Call error(803)
      End If
!
!       NOW GO BACK AND LOOK FOR MORE COMBINE VARIABLES.
!
      Go To 100
140   Continue
!
!       CHECK FOR APPROPRIATE COMMAND TERMINATION AND ENSURE THAT
!       AT LEAST TWO VARIABLES HAVE BEEN SPECIFIED TO COMBINE.
!
      If (list(len:len)/=')') Then
        Call error(804)
      Else If (numind<2) Then
        Call error(805)
      End If
!
!         CALCULATE MINIMUM AND MAXIMUM OF NEW COMBINED VARIABLE
!
      i = 1
      j = 1
      k = 1
      l = 1
      loc1 = 30
      loc2 = 30
      loc3 = 30
      loc4 = 30
      If (comdes(comnum,1)/=99) Then
        loc1 = comdes(comnum, 1)
        i = vardes(loc1, 4) - vardes(loc1, 3) + 1
      End If
      If (comdes(comnum,2)/=99) Then
        loc2 = comdes(comnum, 2)
        j = vardes(loc2, 4) - vardes(loc2, 3) + 1
      End If
      If (comdes(comnum,3)/=99) Then
        loc3 = comdes(comnum, 3)
        k = vardes(loc3, 4) - vardes(loc3, 3) + 1
      End If
      If (comdes(comnum,4)/=99) Then
        loc4 = comdes(comnum, 4)
        l = vardes(loc4, 4) - vardes(loc4, 3) + 1
      End If
      min = 1
      max = ind2(l, k, j, i, l, k, j)
!
!       STICK MIN AND MAX IN THE APPROPRIATE SLOTS OF VARDES.
!
      vardes(varnum, 3) = min
      vardes(varnum, 4) = max
!
!       WRITE TRANSLATION TABLE TO OUTPUT FILE.  START NEW PAGE IF
!       NECESSARY.  THE FOLLOWING CODE WILL WRITE THE HEADING.
!
      ln = ln + 6
      Call page(pg, ln)
      Write (23, 150) varnam(varnum)
150   Format (/, 1X, 'COMBINE has constructed the new variable ', A10, /, 1X, &
        'The following codes were assigned:', /)
      If (numind==2) Then
        Write (23, 160) varnam(varnum), varnam(loc1), varnam(loc2)
      Else If (numind==3) Then
        Write (23, 170) varnam(varnum), varnam(loc1), varnam(loc2), &
          varnam(loc3)
      Else If (numind==4) Then
        Write (23, 180) varnam(varnum), varnam(loc1), varnam(loc2), &
          varnam(loc3), varnam(loc4)
      End If
160   Format (21X, 3(A10,5X), /)
170   Format (13X, 4(A10,5X), /)
180   Format (6X, 5(A10,5X), /)
!
!       NOW WRITE THE TABLE.
!
      Do a = 1, i
        Do b = 1, j
          Do c = 1, k
            Do d = 1, l
              v1 = a + (vardes(loc1,3)-1)
              v2 = b + (vardes(loc2,3)-1)
              v3 = c + (vardes(loc3,3)-1)
              v4 = d + (vardes(loc4,3)-1)
!
!       FUNCTION IND2 CALCULATES THE COMBINED VALUES.
!
              vcom = ind2(d, c, b, a, l, k, j)
!
!       IF WE HAVE STARTED A NEW PAGE OF THE OUTPUT FILE,
!       WRITE A NEW HEADING.
!
              Call page(pg, ln)
              If (ln==1) Then
                ln = ln + 4
                Write (23, *)
                Call page(pg, ln)
                If (numind==2) Then
                  Write (23, 160) varnam(varnum), varnam(loc1), varnam(loc2)
                Else If (numind==3) Then
                  Write (23, 170) varnam(varnum), varnam(loc1), varnam(loc2), &
                    varnam(loc3)
                Else If (numind==4) Then
                  Write (23, 180) varnam(varnum), varnam(loc1), varnam(loc2), &
                    varnam(loc3), varnam(loc4)
                End If
              End If
!
!       WRITE VALUES
!
              If (numind==2) Then
                Write (23, 190) vcom, v1, v2
              Else If (numind==3) Then
                Write (23, 200) vcom, v1, v2, v3
              Else If (numind==4) Then
                Write (23, 210) vcom, v1, v2, v3, v4
              End If
            End Do
          End Do
        End Do
      End Do
      Write (23, *)
190   Format (15X, 3(I10,5X))
200   Format (7X, 4(I10,5X))
210   Format (5(I10,5X))
      Return
    End Subroutine



    Function ind2(i, j, k, l, a, b, c)
!
!       THIS FUNCTION IS USED TO CALCULATE COMBINED VALUES.
!       IT IS CALLED FROM SUBROUTINES COMBINE AND TRANS
!
      Implicit Integer (A-Z)

      ind2 = i + j*a + k*a*b + l*a*b*c - (a+a*b+a*b*c)
      If (ind2<1 .Or. ind2>640000) Then
        Call error(806)
      End If
      Return
    End Function
!
!
!
    Subroutine setup(optyp, setcmd)
!
!       THIS SUBROUTINE INTERPRETS THE SETUP.CMD FILE
!
      Implicit Integer (A-Z)
      Logical exists
      Character *80 setcmd
      Character *35 setnam
      Dimension setcmd(10)
!
!       INITIALIZE THE SETUP COMMAND ARRAY, "SETCMD",
!       AND THE SETUP.CMD LOCATION, "SETNAM"
!
      Do i = 1, 10
        setcmd(i)(1:50) = '                                                  '
      End Do
      setnam = 'UNKNOWN'
!
!       FIND THE SETUP.CMD FILE.
!       FIRST LOOK LOCALLY, THEN AT \DECOMP\ SUBDIRECTORIES
!       ON DRIVES C: D: AND E:
!

      Inquire (File='SETUP.CMD', Exist=exists)
      If (exists) Then
        setnam = 'SETUP.CMD'
      Else
        Inquire (File='C:\DECOMP\SETUP.CMD', Exist=exists)
        If (exists) Then
          setnam = 'C:\DECOMP\SETUP.CMD'
        Else
          Inquire (File='D:\DECOMP\SETUP.CMD', Exist=exists)
          If (exists) Then
            setnam = 'D:\DECOMP\SETUP.CMD'
          Else
            Inquire (File='E:\DECOMP\SETUP.CMD', Exist=exists)
            If (exists) setnam = 'E:\DECOMP\SETUP.CMD'
          End If
        End If
      End If
!
!       IF A SETUP.CMD FILE HAS BEEN FOUND, LOOK FOR SET PROMPT AND
!       SET PAGENUMS COMMANDS.
!
      If (setnam(1:7)/='UNKNOWN') Then
        Open (2, File=setnam, Status='OLD')
        Do i = 1, 10
          Read (2, 120, End=110) setcmd(i)
          If (setcmd(i)(1:3)/='set' .And. setcmd(i)(1:3)/='SET') Go To 100
          p = index(setcmd(i), 'PROMPT')
          If (p==0) p = index(setcmd(i), 'prompt')
          If (p/=0) Then
            p = index(setcmd(i), 'OFF')
            If (p==0) p = index(setcmd(i), 'off')
!
!          IF SET PROMPT = OFF, ASSIGN 2 TO OPTYP
!
            If (p/=0) optyp = 2
            setcmd(i)(1:50) = &
              '                                                  '
          Else
            p = index(setcmd(i), 'PAGE')
            If (p==0) p = index(setcmd(i), 'page')
            If (p/=0) Then
              p = index(setcmd(i), 'OFF')
              If (p==0) p = index(setcmd(i), 'off')
!
!             IF SET PAGENUMS = OFF, TELL THE "PAGE" SUBROUTINE
!
              If (p/=0) Call page(9999, 9999)
              setcmd(i)(1:50) = &
                '                                                  '
            Else
              p = index(setcmd(i), 'SCREEN')
              If (p==0) p = index(setcmd(i), 'screen')
              If (p/=0) Then
                p = index(setcmd(i), 'OFF')
                If (p==0) p = index(setcmd(i), 'off')
!
!               IF SET SCREEN = OFF, OPEN A SCRATCH FILE FOR SCREEN OUTPUT
!
!cc                IF (P.NE.0)OPEN(6,FILE='ERRORS.TMP')
                setcmd(i)(1:50) = &
                  '                                                  '
              End If
            End If
          End If

100     End Do
      End If
110   Continue
120   Format (A80)
      Return
    End Subroutine
