      subroutine acinp(ierr)
      parameter(maxnc=200)
C
C
C IF ACFILE IS SPECIFIED THEN THIS READS IT FREE FORMAT
C    FIELD 1 - ACVAL
C    FIELD 2 - THE CONTOUR INTERVAL APPLYING TO THE ACVAL
C              (THE MAXIMUM NUMBER OF CONTOUR LINES PER INCH
C               IS CALCULATED BASED ON THIS CONTOUR INTERVAL)
C    FIELD 3 - THE MAXIMUM GRADIENT FOR THE ACVAL IN UNITS OF
C              CONTOUR LINES PER INCH (CONTOUR INTERVAL FOR
C              CALCULATING THIS IS THE INTERVAL IN FIELD 2)
C    FIELD 4 - THE TYPE OF LINE TO BE DRAWN FOR THIS ACVAL
C    FIELD 5 - THE SIZE OF THE LABELS
C    FIELD 6 - THE MINIMUM CONTOUR LEVEL WHEN ACDEL IS NEGATIVE.
C    FIELD 7 - THE MAXIMUM CONOTUR LEVEL WHEN ACDEL IS NEGATIVE.
C
      common /concom/ ncol,nrow,bmin,bmax,
     1 grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     1 fltmax,lmult(0:3),idashs,linet,sigma,nsig
      common /labcom/ lchars(6),fmtc,nchar,size,idum(3),delb
      common /contrc/ cmin,cmax,dcval,ncval,nsec,gradi,
     1                acdel(maxnc),acgrad(maxnc),jcdash(maxnc),
     2                ncval2,acsize(maxnc),acmin(maxnc),acmax(maxnc)
      common /contur/ acval(maxnc)
      common/files/ifile,ifile2,ifile3,acfile,mxhach,access
c
      character*40 ifile,ifile2,ifile3,acfile,jblank
      character*16 fmtc
      character*1 access

      REAL*4 ACFARR(7,maxnc),JCFARR(7,maxnc)
      INTEGER*4 LNKACF(maxnc+10),IKEY(2)
      character*20 fields(40),iblank
      logical exists
      equivalence(acfarr(1,1),jcfarr(1,1))
      data dval/1.0e37/,iblank/' '/,jblank/' '/,icmd/7/

      ierr=0
      if(acfile.eq.jblank) go to 35
      open(unit=icmd,file=acfile,status='old',form='formatted',
     & err=985,iostat=icheck)
      go to 3331
985   inquire(file=acfile,exist=exists)
      if(.not.exists) then
      print *,'can''t find file ',acfile
      ierr=1
      return
      else
      write(*,'(a,i5,2a)') 'error no. ',icheck,
     & ' on opening file: ',acfile
      close(icmd)
      ierr=1
      return
      endif
3331  open(24,status='scratch')
      ncval=0
   31 ncval=ncval+1
   36 call freer(icmd,nfield,fields,0,*32,*32,ierr)
      if(ierr.ne.0) return
      if(nfield.lt.1) goto 36
      do 37 i=1,nfield
        k=ichar(fields(i)(1:1))
        if(k.lt.43.or.k.gt.57.or.k.eq.47) goto 38
   37 continue
   38 nfield=i-1
      if(nfield.lt.1) goto 36
      if(nfield.lt.7) then
        do 39 i=nfield+1,7
          fields(i)=iblank
   39   continue
      endif
      goto (71,72,73,74,75,76,77),nfield
   77 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval),acmin(ncval),acmax(ncval)
      go to 711
   76 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval),acmin(ncval)
      go to 711
   75 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval)
      go to 711
   74 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval)
      go to 711
   73 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval)
      go to 711
   72 read(24,*) acval(ncval),acdel(ncval)
      go to 711
   71 read(24,*) acval(ncval)
  711 if(ncval.eq.1) then
        goto (51,52,53,54,55,56,57),nfield+1
   51   goto 58
   52   acdel(ncval)=1.e+7
   53   acgrad(ncval)=abs(gradi)
   54   jcdash(ncval)=0.
   55   acsize(ncval)=size
   56   acmin(ncval)=0.
   57   acmax(ncval)=0.
   58   goto 31
      endif
      goto (61,62,63,64,65,66,67),nfield+1
   61 goto 68
   62 acdel(ncval)=1.e+7
   63 acgrad(ncval)=abs(gradi)
   64 jcdash(ncval)=0.
   65 acsize(ncval)=size
   66 acmin(ncval)=0.
   67 acmax(ncval)=0.
   68 goto 31
C
   32 NCVAL=NCVAL-1
      CLOSE(UNIT=ICMD)
      close(24)
C
C SORT THE ACVAL, ACDEL, ACGRAD, JCDASH ,ACSIZE, ACMIN,
C AND ACMAX ARRAYS
C
   35 if(ncval.gt.0) then
        do 41 j=1,maxnc
          if(acdel(j).lt.0) acval(j)=acdel(j)
   41   continue
        if(acdel(1).eq.0.) then
          do 42 j=1,maxnc
            acdel(j)=1.e+7
   42   continue
        end if
        if(acgrad(1).eq.0.) then
          do 43 j=1,maxnc
            acgrad(j)=abs(gradi)
   43   continue
        end if
        do 44 j=1,maxnc
          acgrad(j)=acdel(j)*abs(acgrad(j))
          if(acgrad(j).gt.1.e+9) acgrad(j)=1.e+9
   44   continue
        if(jcdash(1).eq.999) then
          do 45 j=1,maxnc
            jcdash(j)=iabs(idashs)
   45   continue
        end if
        if(acsize(1).lt.0.) then
          do 46 j=1,maxnc
            acsize(j)=size
   46   continue
        end if
        if(acmin(1).gt.dval) then
          do 471 j=1,maxnc
            acmin(j)=0.
  471   continue
        end if
        if(acmax(1).lt.-dval) then
          do 48 j=1,maxnc
            acmax(j)=0.
   48   continue
        end if
        nkey=2
        ikey(1)=3
        ikey(2)=1
        do 49 j=1,ncval
          acfarr(1,j)=acval(j)
          acfarr(2,j)=acdel(j)
          acfarr(3,j)=acgrad(j)
          jcfarr(4,j)=iabs(jcdash(j))
          acfarr(5,j)=acsize(j)
          acfarr(6,j)=acmin(j)
          acfarr(7,j)=acmax(j)
   49   continue
        call sortm(acfarr,7,maxnc,ncval,ikey,nkey,lnkacf,210,
     1   ier2)
        ll=lnkacf(ncval+1)
        do 50 j=1,ncval
          acval(j)=acfarr(1,ll)
          acdel(j)=acfarr(2,ll)
          acgrad(j)=acfarr(3,ll)
          jcdash(j)=jcfarr(4,ll)
          acsize(j)=acfarr(5,ll)
          acmin(j)=acfarr(6,ll)
          acmax(j)=acfarr(7,ll)
          ll=lnkacf(ll)
   50   continue
        do 81 j=1,ncval
          if(acdel(j).ge.0) goto 33
          acfarr(1,j)=acval(j)
          acfarr(2,j)=acdel(j)
          acfarr(3,j)=acgrad(j)
          acfarr(4,j)=jcdash(j)
          acfarr(5,j)=acsize(j)
          acfarr(6,j)=acmin(j)
          acfarr(7,j)=acmax(j)
   81   continue
   33   ncval2=j-1
        if(ncval2.le.0) goto 34
        nkey=1
        ikey(1)=2
        call sortm(acfarr,7,maxnc,ncval2,ikey,nkey,lnkacf,210,
     1   ier2)
        ll=lnkacf(ncval2+1)
        do 82 j=1,ncval2
          acval(j)=acfarr(1,ll)
          acdel(j)=abs(acfarr(2,ll))
          acgrad(j)=abs(acfarr(3,ll))
          jcdash(j)=jcfarr(4,ll)
          acsize(j)=acfarr(5,ll)
          acmin(j)=acfarr(6,ll)
          acmax(j)=acfarr(7,ll)
          ll=lnkacf(ll)
   82   continue
   34   continue
      end if
      return
      end
C
C      SUBROUTINE FREERD
C________________________________________________________________
C SUBROUTINE FREE READ UTILIZES ASCIIP TO READ
C AN ASCII RECORD AND BREAK IT INTO FIELDS. FIELS MAY BE
C SEPARATED BY ANY COMBINATION OF BLANKS AND COMMMAS.  CARRIAGE
C RETURNS (END-OF-RECORD) MAY INDICATE FIELD DELIMITERS, EOR, OR
C EOF DEPENDING ON THE VALUE OF NREAD AND NFIELD.  DITTO MARKS
C (A SINGLE DECIMAL (.)) MAY BE USED TO ASSIGN THE VALUE FROM THE
C LAST (MOST RECENT) CALL TO FREEREAD3 IN THE CORROSPONDING FIELD
C TO THE FIELD OF THE DITTO MARK IN THE PRESENT CALL.
C
C THE MAXIMUM ACCUMULATIVE FIELD LENGTH FOR ANY ONE CALL IS
C 80 CHARACTERS.  THIS SUBROUTINE IS INTENDED FOR TERMINAL
C ENTRY OR SMALL ASCII FILES WHICH REQUIRE GREAT READING
C FLEXIBILITY.
C
C   IUNIT  - INTEGER*4.  UNIT NUMBER TO BE READ.  THIS MAY BE A
C            TERMINAL KEYBOARD OR ANY FILE (DISK, TAPE, ETC).
C   NFIELD - INTEGER*4.  TOTAL NUMBER OF FIELDS ACTUALLY READ.
C   FIELDS - CHARACTER*(*).  CHARACTER ARRAY CONTAINING THE
C            FIELDS WHICH ARE READ IN.  EACH ELEMENT OF THE ARRAY
C            CONTAINS ONE FIELD.
C   NREAD  - INTEGER*4.  DETERMINES WHETHER A CARRIAGE RETURN IS
C            A FIELD DELIMITER (TYPICALLY SPACE, COMMA),
C            AN END-OF-REC MARKER (TYPICALLY CARRIAGE RETURN),
C            OR AN END-OF-FILE MARKER (TYPICALLY CONTROL-Z,
C            E-O-F, EXIT).  THE FIELD DELIMITER INDICATES
C            BEGINNING OF A NEW FIELD AND CONTROL REMAINS IN
C            FREEREAD3.  THE END-OF-REC MARKER INDICATES THAT
C            ALL FIELDS FOR THIS RECORD HAVE BEEN ACCUMULATED AND
C            CONTROL IS TRANSFERED BACK TO THE CALLING ROUTINE.
C            THE END-OF-FILE MARKER INDICATES THAT THERE IS NO
C            DATA AVAILABLE ON THIS CALL AND CONTROL IS
C            TRANSFERED BACK TO THE CALLING ROUTINE VIA THE
C            STARRED (*) RETURN.  FOLLOWING IS A DESCRIPTION OF
C            HOW A CARRIAGE RETURN (OR END OF RECORD) BEHAVES AS
C            A FUNCTION OF THE VALUE OF NREAD AND THE NUMBER OF
C            FIELDS READ IN (NFIELD):
C
C              NREAD > 0:  0 =< NFIELD < NREAD:  FIELD DELIM.
C                          NFIELD >= NREAD:  END-OF-REC.
C              NREAD = 0:  FOR ALL NFIELD:  END-OF-REC.
C              NREAD < 0:  NFIELD = 0:  END-OF-FILE *2.
C                          0 < NFIELD < ABS(NREAD):  FIELD DELIM.
C                          NFIELD >= ABS(NREAD):  END-OF-REC.
C                          LAST FIELD WAS 'EXIT':  END-OF-FILE *1
C
C              FOR ALL NREAD AND NFIELD CONTROL-Z AND E-O-F ARE
C              END-OF-FILE *1 MARKERS.  WARNING:  CONTROL-Z DOES
C              SOME UNPREDICTABLE THINGS IN ADDITION TO WHAT IT'S
C              SUPPOSED TO DO.  IT IS RECOMMENDED THAT THE USER
C              AVOID IT ALL TOGETHER.
C
C   *1     - LINE NUMBER TO WHICH CONTROL IS RETURNED IN THE
C            CASE OF AN END-OF-FILE *1
C   *2     - LINE NUMBER TO WHICH CONTROL IS RETURNED IN THE
C            CASE OF AN END-OF-FILE *2
C
C SUBROUTINE FREEREAD3 WRITTEN BY ROB BRACKEN, USGS, 20JUL87.
C
C
      subroutine freer(iunit,nfield,fields,nread,*,*,ier)
      save
C
      character*(*) fields(1)
      character*80 alphnu,alphn2
      character*5 ftest
      character*1 ditto
      data alphnu /' '/
      ier=0
C
C REMEMBER LAST READ
C
      alphn2=alphnu
C
C READ THE FILE
C
      i=1
      nfield=0
  102 read(iunit,801,end=998,err=900,iostat=icheck) alphnu
  801 format(a80)
C
C PARSE THE LINE FROM THE FILE
C
      ditto=' '
      call asciir(80,alphnu,ditto,nfield,fields(i))
C
C DETERMINE IF THE PROPER NUMBER OF FIELDS HAVE BEEN READ
C
      i=i+nfield
      nfield=i-1
      if(nread.lt.0) then
        if(i.eq.1) goto 997
        ftest=fields(i-1)
        call upper(ftest,5)
        if(ftest.eq.'EXIT ') goto 998
      endif
      if(nread.ne.0.and.i.le.abs(nread)) goto 102
C
C FILL IN THE DITTO MARKS
C
      ditto='.'
      call asciir(80,alphn2,ditto,nfld2,fields)
C
C SET UP ALPHNU WITH PRESENT SET OF FIELDS
C
      lent=0
      do 20 j=1,nfield
        do 10 len=1,80
          if(fields(j)(len:len).eq.' ') goto 101
   10   continue
  101   alphnu(lent+1:lent+len)=fields(j)(1:len)
        lent=lent+len
   20 continue
      if(lent.lt.80) alphnu(lent+1:80)=' '
C
C CUT OUT ANY FIELDS BEYOND THE LAST FIELD IN THIS READ
C
      if(nfield.lt.nfld2) then
        do 30 i=nfield+1,nfld2
          fields(i)=' '
   30   continue
      endif
C
C NORMAL RETURN
C
      rewind(24)
      write(24,*) (fields(i),i=1,nfield)
      rewind(24)
  999 return
C
C END OF FILE REACHED
C
  900 print *,'error no. ',icheck,' processing free form input'
      ier=1
      return
  998 alphnu=' '
      return 1
  997 alphnu=' '
      return 2
      end
C________________________________________________________________
C
C     SUBROUTINE  A S C I I P
C________________________________________________________________
C
C SUBROUTINE ASCII PARSE TAKES A CHARACTER VARIABLE OF ANY LENGTH
C AND DETERMINES WHERE THE DATA FIELDS ARE.  THE FIELDS MUST BE
C SEPARATED BY A BLANK SPACE OR A COMMA (AS IN THE VMS STANDARD).
C
C   MXLEN  - INTEGER*4.  MAXIMUM ALLOWABLE LENGTH OF THE INPUT
C            CHARACTER VARIABLE.
C   ALPHNU - CHARACTER*(*).  CHARACTER (ALPHA-NUMERIC) VARIABLE
C            TO BE PARSED.
C   DITTO  - CHARACTER*1.  THE SYMBOL WHICH IS TO BE USED AS A
C            DITTO MARK.  IF DITTO IS BLANK THEN NO DITTO
C            FUNCTION WILL BE PERFORMED.  SEE "FIELDS".
C   NFIELD - INTEGER*4.  THE NUMBER OF FIELDS WHICH WERE READ IN.
C   FIELDS - CHARACTER*(*).  CHARACTER (ALPHA-NUMERIC) VARIABLE
C            ARRAY TO WHICH THE FIELDS WILL BE READ FROM ALPHNU.
C            EACH ARRAY POSITION IN "FIELDS" WILL CONTAIN A
C            SINGLE FIELD.  IF DITTO IS BLANK ALL FIELDS WILL BE
C            READ FROM ALPHNU TO "FIELDS".  IF DITTO IS NOT BLANK
C            THEN ONLY THOSE ARRAY POSITIONS IN "FIELDS" WHICH
C            CONTAIN A SINGLE DITTO CHARACTER WILL BE FILLED.
C
C SUBROUTINE ASCIIPARSE WRITTEN BY ROB BRACKEN, USGS, 20JUL87.
C
C
      subroutine asciir(mxlen,alphnu,ditto,nfield,fields)
C
      character*(*) alphnu,fields(1)
      character*1 ditto,apos
      logical*4 blank,text,iapos
      data apos/''''/
C
C INIALIZE VARIABLES
C
      blank=.true.
      iapos=.false.
      text=.false.
      nfield=0
C
C SET UP ARRAY CONTAINING FIELDS.
C
      do 10 i=1,mxlen
      if(alphnu(i:i).eq.apos) then
       if(.not.text) then
       if(.not.iapos) then
C      FIRST APOSTROPHY READ
       text=.true.
       ibeg=i
       else
C      TWO APOSTROPHYS IN A ROW SO IT IS PART OF TEXT
       text=.true.
       nfield=nfield-1
       endif
      else
C      SECOND APOSTROPHY READ-TEXT FIELD COMPLETE
       text=.false.
       nfield=nfield+1
       if(ditto.eq.' '.or.fields(nfield).eq.ditto)
     1 fields(nfield)=alphnu(ibeg:i)
       blank=.true.
       iapos=.true.
      endif
        else if(alphnu(i:i).eq.' '.or.alphnu(i:i).eq.','
     1          .and.(.not.text)) then
          iapos=.false.
          if(.not.blank) then
C ---       PREVIOUS FIELD NON-BLANK/PRESENT FIELD BLANK ---
            nfield=nfield+1
            if(ditto.eq.' '.or.fields(nfield).eq.ditto)
     1       fields(nfield)=alphnu(ibeg:i-1)
            blank=.true.
          end if
        else
          if(blank.and.(.not.text)) then
C ---       PREVIOUS FIELD BLANK/PRESENT FIELD NON-BLANK ---
            ibeg=i
            blank=.false.
          end if
        end if
   10 continue
      if(.not.blank) then
C ---   FINAL FIELD NON-BLANK ---
        nfield=nfield+1
        if(ditto.eq.' '.or.fields(nfield).eq.ditto)
     1   fields(nfield)=alphnu(ibeg:i-1)
        blank=.true.
      end if
C
C EXIT PROCEDURE
C
  199 return
      end
C________________________________________________________________
C
C     SUBROUTINE  U P P E R C 
C________________________________________________________________
C
C SUBROUTINE UPPER CASE CONVERTS ALPHA-NUMERIC VARIABLES INTO
C ALL UPPER CASE LETTERS.
C
C   ALPHA  - CHARACTER*(*).  ALPHA-NUMERIC TO BE CONVERTED.
C   LENGTH - INTEGER*4.  LENGTH OF ALPHA.
C
C SUBROUTINE UPPER WRITTEN BY ROB BRACKEN, USGS, 17JUL87.
C
C
      subroutine upper(alpha,length)
C
      character*(*) alpha
C
      do 10 i=1,length
        k=ichar(alpha(i:i))
        if(k.ge.97.and.k.le.122) then
          alpha(i:i)=char(k-32)
        end if
   10 continue
      return
      end
C________________________________________________________________
C
C     SUBROUTINE   S O R T I L
C________________________________________________________________
C
c  Hart modification to Wooodrum sort from Creative Computing
c  1978,vol 4,no. 1,p. 96-101.
c  Translated to FORTRAN by R.R. Wahl UTR Branch
c  US Geol. Survey, Denver, CO
c
c     modified by r.godson,usgs
c     to handle multi-dimensioned input arrays
c     and return sorted data in another array
c     equivalent in size to the input array    
c
c     return array option is deleted in this version
c
C
C  DESCRIPTION OF PASSED VARIABLES:
C
C     AIN    - REAL*4.  2 DIMENSIONAL INPUT UNSORTED ARRAY.
C     IDIM   - INTEGER*4.  FIRST DIMENSION OF AIN.
C     JDIM   - INTEGER*4.  SECOND DIMENSION OF AIN.
C     NSORT  - INTEGER*4.  NUMBER OF COLUMNS TO SORT.
C     IKEY   - INTEGER*4.  1 DIMENSIONAL ARRAY CONTAINING A LIST
C              OF KEY FIELDS (ROW NUMBERS, FIRST INDEX) TO SORT
C              FROM MAJOR TO MINOR.
C     NKEY   - INTEGER*4.  THE NUMBER OF KEY FIELDS.
C     LNK    - INTEGER*4.  1 DIMENSIONAL OUTPUT LINK ARRAY.
C              AIN(ISORT,LNK(NSORT+1)) = LEAST VALUE.
C              AIN(ISORT,LNK(LNK(NSORT+1))) = SECOND LEAST VALUE.
C     LDIM   - INTEGER*4.  DIMENSION OF LINK ARRAY.
C              MINIMUM LDIM = NSORT + LN( NSORT-2 ) + 2
C     IER    - INTEGER*4.  ERROR INDICATOR.  0 = SUCCESSFUL RUN.
C              1 = LESS THAN 2 ELEMENTS TO SORT.  2 = NOT ENOUGH
C              WORK SPACE IN LINK ARRAY.
C
C
      subroutine sortm(ain,idim,jdim,
     1                  nsort,ikey,nkey,lnk,ldim,ier)
C
      dimension ain(idim,jdim),lnk(ldim),ikey(nkey)
c
      ier=0
      if (nsort.le.1)  go to 1000
      mt=nsort+ifix(alog(float(nsort))/.69314718+.5)+2
      if (ldim.lt.mt)  go to 1010
c
c  Initialize variables
      k1=0
      i=0
      m1=0
      it2=0
      it4=0
c     Head of sequence 1
      j=nsort+1
      lnk(1)=1
      lnk(j)=1
      k2=1
c     Number of leaves
      nleav=nsort
c
c  Climb the tree
  250 continue
c     Low-order twig value
      if (nleav.lt.4)  go to 320
      ib2=nleav/2
      it4=it4+k2*(nleav-2*ib2)
      nleav=ib2
c     Total number of twigs at an altitude
      k2=k2*2
      go to 250
c
c  Initial calculations
  320 continue
c     Number of low-order twigs
      it4=k2-it4
c     High bit value of binary counter
      ib2=k2/2
c
c  Next twig
  350 continue
c     Exit -- sort complete
      if (k1.eq.k2)  go to 940
c     Twig number
      k1=k1+1
      it1=k1
c     High bit value
      ib1=ib2
c     Previous reflected twig number
      it3=it2
c
c  Add 1 to reflected binary counter and carry
  400 continue
c     Any more carries
      it1s=it1/2
      if (2*it1s.lt.it1)  go to 470
      it1=it1s
c     Number of merges
      m1=m1+1
      it2=it2-ib1
c     Next bit value
      ib1=ib1/2
c     Carry one
      go to 400
c
c  Twig calculations
  470 continue
c     Reflected twig number
      it2=it2+ib1
c     2-twigs and 3-twigs
      if (nleav.eq.2)  go to 550
c
c  3-twigs and 4-twigs
c     low-order twig (3-twig)
      if (it3.lt.it4)  go to 560
c
c  4-twig
c     Dis-engage number of merges
      m1=-m1
      go to 630
  550 continue
c     Low-order twig (2-twig)
      if (it3.lt.it4)  go to 610
c
c  3-twig
  560 continue
c     Number of merges
      m1=m1+1
c     Next leaf
      i=i+1
c     Generate a leaf
      lnk(i)=i
      lnk(j)=i
c     Generate next sequence head
      j=j+1
c
c  2-twig
  610 continue
c     Number of merges
      m1=m1+1
  630 continue
c     Next leaf
      i=i+1
c     Generate a new leaf
      l1=i
      lnk(i)=i
      lnk(j)=i
c     Head of older leaf (last line)
      l0=j
c     Head of latest leaf (next two lines)
      j=j+1
c     Next leaf
      i=i+1
c     Generate a leaf
      l2=i
      lnk(i)=i
      lnk(j)=i
c     Merge leaves
      go to 750
c
c  Merge twigs and branches
  700 continue
c     Head of latest branch or twig
      j=j-1
c     Head of older branch or twig
      l0=j-1
c     Head of sequence 1
      l1=lnk(l0)
c     Head of sequence 2
      l2=lnk(j)
  750 continue
c     Stay in sequence 1
C
      kkey=1
  763 if(ain(ikey(kkey),l1)-ain(ikey(kkey),l2)) 820,762,761
  762 kkey=kkey+1
      if(kkey.le.nkey) goto 763
      goto 820
  761 continue
C
C      if (ain(isort,l1).le.ain(isort,l2))  go to 820
c     Switch to sequence 2
      lnk(l0)=l2
  770 continue
c     Top leaf in sequence 2
      l0=l2
c     Next leaf in sequence 2
      l2=lnk(l0)
c     End of sequence 2
      if (l2.eq.l0)  go to 870
c     Stay in sequence 2
C
      kkey=1
  783 if(ain(ikey(kkey),l1)-ain(ikey(kkey),l2)) 781,782,770
  782 kkey=kkey+1
      if(kkey.le.nkey) goto 783
  781 continue
C
C      if (ain(isort,l1).gt.ain(isort,l2))  go to 770
c     Switch to sequence 1
      lnk(l0)=l1
  820 continue
c     Top leaf in sequence 1
      l0=l1
c     Next leaf in sequence 2
      l1=lnk(l0)
c     Not end of sequence 1
      if (l1.ne.l0)  go to 750
c     Switch to sequence 2
      lnk(l0)=l2
      go to 880
  870 continue
c     Switch to sequence 1
      lnk(l0)=l1
  880 continue
c     Number of merges
      m1=m1-1
      if (m1.gt.0)  go to 700
      if (m1.eq.0)  go to 350
c
c  Generate second half of 4-twig
c     Re-engage number of merges
      m1=1-m1
      go to 630
c
c  All done
c
c     sort output array
c
  940 ll=lnk(nsort+1)
c      do 960 i=1,nsort
c      do 950 j=1,jdim
c      aout(i,j)=ain(ll,j)
c  950 continue
c      ll=lnk(ll)
c  960 continue
  970 return
c     Less than two elements to sort
 1000 ier=1
      lnk(nsort+1)=nsort
      lnk(nsort)=nsort+1
      go to 970
c     Not enough work space in link array
 1010 ier=2
      go to 970
      end
